|
| 1 | +/* vim:set noexpandtab tabstop=4 wrap */ |
| 2 | +#include "WCSimDemo.h" |
| 3 | + |
| 4 | +WCSimDemo::WCSimDemo():Tool(){} |
| 5 | + |
| 6 | + |
| 7 | +bool WCSimDemo::Initialise(std::string configfile, DataModel &data){ |
| 8 | + |
| 9 | + /////////////////// Usefull header /////////////////////// |
| 10 | + //loading config file |
| 11 | + if(configfile!="") m_variables.Initialise(configfile); |
| 12 | + //m_variables.Print(); |
| 13 | + |
| 14 | + //assign transient data pointer |
| 15 | + m_data= &data; |
| 16 | + |
| 17 | + // Get the Tool configuration variables |
| 18 | + m_variables.Get("verbosity",verbosity); |
| 19 | + |
| 20 | + Log("Initializing Tool WCSimDemo",v_message,verbosity); |
| 21 | + return true; |
| 22 | +} |
| 23 | + |
| 24 | + |
| 25 | +bool WCSimDemo::Execute(){ |
| 26 | + Log("WCSimDemo Tool: Executing",v_debug,verbosity); |
| 27 | + get_ok = m_data->Stores.count("ANNIEEvent"); |
| 28 | + if(!get_ok){ |
| 29 | + Log("WCSimDemo Tool: No ANNIEEvent store!",v_error,verbosity); |
| 30 | + return false; |
| 31 | + }; |
| 32 | + |
| 33 | + // First, see if this is a delayed trigger in the event |
| 34 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCTriggernum",MCTriggernum); |
| 35 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCTriggernum from ANNIEEvent!",v_error,verbosity); return false; } |
| 36 | + // if so, truth analysis is probably not interested in this trigger. Primary muon will not be in the listed tracks. |
| 37 | + if(MCTriggernum>0){ Log("WCSimDemo Tool: Skipping delayed trigger",v_debug,verbosity); return true;} |
| 38 | + |
| 39 | + // Retrieve the info from the ANNIEEvent |
| 40 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCFile",MCFile); |
| 41 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCFile from ANNIEEvent!",v_error,verbosity); return false; } |
| 42 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCEventNum",MCEventNum); |
| 43 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCEventNum from ANNIEEvent!",v_error,verbosity); return false; } |
| 44 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("EventNumber",EventNumber); |
| 45 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving EventNumber from ANNIEEvent!",v_error,verbosity); return false; } |
| 46 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCParticles",MCParticles); |
| 47 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCParticles,true from ANNIEEvent!",v_error,verbosity); return false; } |
| 48 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCHits",MCHits); |
| 49 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCHits,true from ANNIEEvent!",v_error,verbosity); return false; } |
| 50 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("MCLAPPDHits",MCLAPPDHits); |
| 51 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving MCLAPPDHits,true from ANNIEEvent!",v_error,verbosity); return false; } |
| 52 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("TDCData",TDCData); |
| 53 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving TDCData,true from ANNIEEvent!",v_error,verbosity); return false; } |
| 54 | + get_ok = m_data->Stores.at("ANNIEEvent")->Get("EventTime",EventTime); |
| 55 | + if(not get_ok){ Log("WCSimDemo Tool: Error retrieving EventTime,true from ANNIEEvent!",v_error,verbosity); return false; } |
| 56 | + |
| 57 | + // Print the information. Jingbo, do what you do here. :) |
| 58 | + logmessage = "WCSimDemo Tool: Processing truth tracks and digits for "+MCFile |
| 59 | + +", MCEvent "+to_string(MCEventNum)+", MCTrigger "+to_string(MCTriggernum); |
| 60 | + Log(logmessage,v_debug,verbosity); |
| 61 | + |
| 62 | + // loop over the MCParticles to find the highest enery primary muon |
| 63 | + // MCParticles is a std::vector<MCParticle> |
| 64 | + MCParticle primarymuon; // primary muon |
| 65 | + bool mufound=false; |
| 66 | + if(MCParticles){ |
| 67 | + Log("WCSimDemo Tool: Num MCParticles = "+to_string(MCParticles->size()),v_message,verbosity); |
| 68 | + for(int particlei=0; particlei<MCParticles->size(); particlei++){ |
| 69 | + MCParticle aparticle = MCParticles->at(particlei); |
| 70 | + //if(v_debug<verbosity) aparticle.Print(); // print if we're being *really* verbose |
| 71 | + if(aparticle.GetParentPdg()!=0) continue; // not a primary particle |
| 72 | + if(aparticle.GetPdgCode()!=13) continue; // not a muon |
| 73 | + primarymuon = aparticle; // note the particle |
| 74 | + mufound=true; // note that we found it |
| 75 | + break; // won't have more than one primary muon |
| 76 | + } |
| 77 | + } else { |
| 78 | + Log("WCSimDemo Tool: No MCParticles in the event!",v_error,verbosity); |
| 79 | + } |
| 80 | + if(not mufound){ |
| 81 | + Log("WCSimDemo Tool: No muon in this event",v_warning,verbosity); |
| 82 | + return true; |
| 83 | + } |
| 84 | + |
| 85 | + // retrieve desired information from the particle |
| 86 | + const Position neutrinovtx = primarymuon.GetStartVertex(); // only true if the muon is primary |
| 87 | + const Direction muondirection = primarymuon.GetStartDirection(); |
| 88 | + double muonenergy = primarymuon.GetStartEnergy(); |
| 89 | + logmessage = "WCSimDemo Tool: Interaction Vertex is at ("+to_string(neutrinovtx.X()) |
| 90 | + +", "+to_string(neutrinovtx.Y())+", "+to_string(neutrinovtx.Z())+")\n" |
| 91 | + +"Primary muon has energy "+to_string(muonenergy)+"GeV and direction (" |
| 92 | + +to_string(muondirection.X())+", "+to_string(muondirection.Y())+", "+to_string(muondirection.Z())+")"; |
| 93 | + Log(logmessage,v_debug,verbosity); |
| 94 | + // see Particle.h for other information in the MCParticle class |
| 95 | + |
| 96 | + // now move to digit retrieval |
| 97 | + // MCHits is a std::map<ChannelKey,std::vector<Hit>> |
| 98 | + if(MCHits){ |
| 99 | + Log("WCSimDemo Tool: Num PMT Digits = "+to_string(MCHits->size()),v_message,verbosity); |
| 100 | + // iterate over the map of sensors with a measurement |
| 101 | + for(std::pair<ChannelKey,std::vector<Hit>>&& apair : *MCHits){ |
| 102 | + ChannelKey chankey = apair.first; |
| 103 | + // a ChannelKey is a detector descriptor, containing 2 elements: |
| 104 | + // a 'subdetector' (enum class), with types ADC, LAPPD, TDC |
| 105 | + // and a DetectorElementIndex, i.e. the ID of the detector of that type |
| 106 | + |
| 107 | + if(chankey.GetSubDetectorType()==subdetector::ADC){ |
| 108 | + std::vector<Hit>& hits = apair.second; |
| 109 | + for(Hit& ahit : hits){ |
| 110 | + //if(v_message<verbosity) ahit.Print(); // << VERY verbose |
| 111 | + // a Hit has tubeid, charge and time |
| 112 | + } |
| 113 | + } |
| 114 | + } // end loop over MCHits |
| 115 | + } else { |
| 116 | + cout<<"No MCHits"<<endl; |
| 117 | + } |
| 118 | + |
| 119 | + // repeat for LAPPD hits |
| 120 | + // MCLAPPDHits is a std::map<ChannelKey,std::vector<LAPPDHit>> |
| 121 | + if(MCLAPPDHits){ |
| 122 | + Log("WCSimDemo Tool: Num LAPPD Digits = "+to_string(MCLAPPDHits->size()),v_message,verbosity); |
| 123 | + // iterate over the map of sensors with a measurement |
| 124 | + for(std::pair<ChannelKey,std::vector<LAPPDHit>>&& apair : *MCLAPPDHits){ |
| 125 | + ChannelKey chankey = apair.first; |
| 126 | + if(chankey.GetSubDetectorType()==subdetector::LAPPD){ // redundant |
| 127 | + std::vector<LAPPDHit>& hits = apair.second; |
| 128 | + for(LAPPDHit& ahit : hits){ |
| 129 | + //if(v_message<verbosity) ahit.Print(); // << VERY verbose |
| 130 | + // an LAPPDHit has adds (global x-y-z) position, (in-tile x-y) local position |
| 131 | + // and time psecs |
| 132 | + } |
| 133 | + } |
| 134 | + } // end loop over MCLAPPDHits |
| 135 | + } else { |
| 136 | + cout<<"No MCLAPPDHits"<<endl; |
| 137 | + } |
| 138 | + |
| 139 | + return true; |
| 140 | +} |
| 141 | + |
| 142 | + |
| 143 | +bool WCSimDemo::Finalise(){ |
| 144 | + |
| 145 | + return true; |
| 146 | +} |
0 commit comments