diff --git a/UserTools/LoadWCSim/LoadWCSim.cpp b/UserTools/LoadWCSim/LoadWCSim.cpp index 5d9600108..16ff3636a 100644 --- a/UserTools/LoadWCSim/LoadWCSim.cpp +++ b/UserTools/LoadWCSim/LoadWCSim.cpp @@ -1,1587 +1,1582 @@ -/* vim:set noexpandtab tabstop=4 wrap */ - #include "LoadWCSim.h" LoadWCSim::LoadWCSim():Tool(){} -bool LoadWCSim::Initialise(std::string configfile, DataModel &data){ - - /////////////////// Useful header /////////////////////// - - if(verbosity) cout<<"Initializing Tool LoadWCSim"<CStore.Set("SplitSubTriggers",splitSubtriggers); - get_ok = m_variables.Get("TriggerType",Triggertype); - if (not get_ok){ - Log("LoadWCSim Tool: No Triggertype specified. Assuming TriggerType = Beam",v_warning,verbosity); - Triggertype = "Beam"; //other options: Cosmic / No Loopback - } + if (!m_variables.Get("MaxEntries", MaxEntries)) MaxEntries = -1; + + if (!m_variables.Get("InputFile", MCFile)) { + logmessage = "LoadWCSim::Initialise: NO InputFile set in the config!"; + Log(logmessage, v_error, verbosity); + return false; + } + + if (!m_variables.Get("HistoricTriggeroffset", HistoricTriggeroffset)) { + logmessage = "LoadWCSim::Initialise: NO HistoricTriggeroffset set in the config! Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + + if (!m_variables.Get("WCSimVersion", WCSimVersion)) { + logmessage = "LoadWCSim::Initialise: WCSimVersion not set in the config! Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + m_data->CStore.Set("WCSimVersion", WCSimVersion); + + if (!m_variables.Get("UseDigitSmearedTime", use_smeared_digit_time)) { + use_smeared_digit_time = 1; + logmessage = "LoadWCSim::Initialise: UseDigitSmearedTime flag not set in the config."; + logmessage += "Using default value: " + std::to_string(use_smeared_digit_time); + Log(logmessage, v_warning, verbosity); + + } + + if (!m_variables.Get("LappdNumStrips", LappdNumStrips)) { + LappdNumStrips = 56; + logmessage = "LoadWCSim::Initialise: LappdNumStrips not set in config. "; + logmessage += "Using default value: " + std::to_string(LappdNumStrips); + Log(logmessage, v_warning, verbosity); + } + + if(!m_variables.Get("LappdStripLength", LappdStripLength)){ + LappdStripLength = 200; + logmessage = "LoadWCSim::Initialise: LappdStripLength not set in config. "; + logmessage += "Using default value: " + std::to_string(LappdStripLength); + Log(logmessage, v_warning, verbosity); + } + + if (!m_variables.Get("LappdStripSeparation", LappdStripSeparation)) { + LappdStripSeparation = 7.14; + logmessage = "LoadWCSim::Initialise: LappdStripSeparation not set in config. "; + logmessage += "Using default value: " + std::to_string(LappdStripSeparation); + Log(logmessage, v_warning, verbosity); + } + + if (!m_variables.Get("RunStartDate", RunStartUser)) { + RunStartUser = 0; + logmessage = "LoadWCSim::Initialise: RunStartUser not set in config. "; + logmessage += "Using default value: " + std::to_string(LappdStripSeparation) + " ns since unix epoch"; + Log(logmessage, v_warning, verbosity); + } + + if (!m_variables.Get("SplitSubTriggers", splitSubtriggers)) { + splitSubtriggers = false; + logmessage = "LoadWCSim::Initialise: splitSubtriggers flag not set in config. "; + logmessage += "Using default value: " + std::to_string(splitSubtriggers); + Log(logmessage, v_warning, verbosity); + } + m_data->CStore.Set("SplitSubTriggers", splitSubtriggers); + + if (!m_variables.Get("TriggerType", TriggerType)) { + TriggerType = "Beam"; + logmessage = "LoadWCSim::Initialise: TriggerType not set in config. "; + logmessage += "Using default value: " + TriggerType; + Log(logmessage, v_warning, verbosity); + } - get_ok = m_variables.Get("TriggerWord",TriggerWord); - if (not get_ok){ - Log("LoadWCSim Tool: No Triggerword specified. Assuming TriggerWord = 5 (Beam)",v_warning,verbosity); - TriggerWord = 5; - } - path_chankeymap = "./configfiles/LoadWCSim/Chankey_WCSimID_v7.txt"; - get_ok = m_variables.Get("ChankeyToPMTIDMap",path_chankeymap); - if (not get_ok){ - Log("LoadWCSim Tool: No Channelkey map provided. Use the standard one at ./configfiles/LoadWCSim/Chankey_WCSimID_v7.txt",v_warning,verbosity); - } - ifstream file_pmtid(path_chankeymap.c_str()); - if (file_pmtid.is_open()){ - // watch out: comment or empty lines not supported here - while (!file_pmtid.eof()){ - unsigned long chankey; - int pmtid; - file_pmtid >> chankey >> pmtid; - channelkey_to_pmtid.emplace(chankey,pmtid); - pmtid_to_channelkey.emplace(pmtid,chankey); - } + if (!m_variables.Get("TriggerWord", TriggerWord)) { + TriggerWord = 5; + logmessage = "LoadWCSim::Initialise: TriggerWord not set in config. "; + logmessage += "Using default value: " + std::to_string(TriggerWord) + " (ie. Beam)"; + Log(logmessage, v_warning, verbosity); + } + + if (!m_variables.Get("RunType", RunType)) { + RunType = 3; + logmessage = "LoadWCSim::Initialise: RunType not set in config. "; + logmessage += "Using default value: " + std::to_string(RunType) + " (ie. Beam)"; + Log(logmessage, v_warning, verbosity); + } + + if (!m_variables.Get("PMTMask", PMTMask)) { + PMTMask = "None"; + logmessage = "LoadWCSim::Initialise: RunType not set in config. "; + logmessage += "Using default value: " + PMTMask; + Log(logmessage, v_warning, verbosity); + } + if (PMTMask != "None") masked_ids = this->LoadPMTMask(PMTMask); + + if (!m_variables.Get("FileStartOffset", WCSimEntryNum)) { + WCSimEntryNum = 0; + logmessage = "LoadWCSim::Initialise: FileStartOffset not set in config. "; + logmessage += "Using default value: " + std::to_string(WCSimEntryNum); + Log(logmessage, v_warning, verbosity); + } + + // TODO: should not use relative path. + // There should be an env var that points to the top directory + if (!m_variables.Get("ChankeyToPMTIDMap", path_chankeymap)) { + path_chankeymap = "./configfiles/LoadWCSim/Chankey_WCSimID_v7.txt"; + logmessage = "LoadWCSim::Initialise: ChankeyToPMTIDMap not set in config. "; + logmessage += "Using default path: " + path_chankeymap; + Log(logmessage, v_warning, verbosity); + } + ifstream file_pmtid(path_chankeymap.c_str()); + if (file_pmtid.is_open()){ + // watch out: comment or empty lines not supported here + while (!file_pmtid.eof()) { + unsigned long chankey; + int pmtid; + file_pmtid >> chankey >> pmtid; + channelkey_to_pmtid.emplace(chankey,pmtid); + pmtid_to_channelkey.emplace(pmtid,chankey); + } - file_pmtid.close(); - m_data->CStore.Set("pmt_tubeid_to_channelkey_data",pmtid_to_channelkey); - m_data->CStore.Set("channelkey_to_pmtid_data",channelkey_to_pmtid); - } else { - Log("LoadWCSim Tool: PMT ID Configuration file "+path_chankeymap+" could not be opened! Is the path valid? Abort",v_warning,verbosity); - return false; - } - path_mrd_chankeymap = "./configfiles/LoadWCSim/MRD_Chankey_WCSimID.dat"; - get_ok = m_variables.Get("ChankeyToMRDIDMap",path_mrd_chankeymap); - ifstream file_mrdid(path_mrd_chankeymap.c_str()); - if (file_mrdid.is_open()){ - // watch out: comment or empty lines not supported here - while (!file_mrdid.eof()){ - unsigned long chankey; - int mrdid; - file_mrdid >> chankey >> mrdid; - mrdid_to_channelkey.emplace(mrdid,chankey); - } - file_mrdid.close(); - } else { - Log("LoadWCSim Tool: MRD ID Configuration file "+path_mrd_chankeymap+" could not be opened! Is the path valid? Abort",v_warning,verbosity); - return false; - } - path_fmv_chankeymap = "./configfiles/LoadWCSim/FMV_Chankey_WCSimID.dat"; - get_ok = m_variables.Get("ChankeyToFMVIDMap",path_fmv_chankeymap); - ifstream file_fmvid(path_fmv_chankeymap.c_str()); - if (file_fmvid.is_open()){ - // watch out: comment or empty lines not supported here - while (!file_fmvid.eof()){ - unsigned long chankey; - int fmvid; - file_fmvid >> chankey >> fmvid; - fmvid_to_channelkey.emplace(fmvid,chankey); - } - file_fmvid.close(); - } else { - Log("LoadWCSim Tool: FMV ID Configuration file "+path_fmv_chankeymap+" could not be opened! Is the path valid? Abort",v_warning,verbosity); - return false; - } - get_ok = m_variables.Get("RunType",RunType); - if (not get_ok){ - Log("LoadWCSim Tool: No RunType specified. Assuming RunType = 3 (Beam)",v_warning,verbosity); - RunType = 3; - } - - get_ok = m_variables.Get("PMTMask",PMTMask); - if (not get_ok){ - Log("LoadWCSim Tool: Assuming to use no PMTMask",v_warning,verbosity); - PMTMask = "None"; - } - - MCEventNum=0; - get_ok = m_variables.Get("FileStartOffset",MCEventNum); - - // put version in the CStore for downstream tools - m_data->CStore.Set("WCSimVersion", WCSimVersion); - - // Short Stores README - ////////////////////// - // n.b. m_data->vars is a Store (of ben's Store type) that is not saved to disk? - // m_data->CStore is a single entry binary BoostStore that is not saved to disk. - // m_data->Stores["StoreName"] is a map of binary BoostStores that are saved to disk. - // If using Stores->BoostStore->Set("MyVariable") it will always be saved to disk - // Using Stores->BoostStore.Set("MyVariable",&myvar) if myvar is a pointer (to an object on - // the heap) puts myvar in the Store and it's deletion will be handled by the Store. - // (provided your class has a suitable destructor.) - // Is 'BoostStore::Save' needed for single-entry stores? - // ---------------- - // create a new BoostStore with key "ANNIEEvent" in the Stores std::map - // BoostStore constructor args: typechecking (bool), m_format (0=binary, 1=ASCII, 2=multievent) - // A BoostStore has a header where useful constants may be saved. The header is a BoostStore itself, - // and can be accessed via: Store.Header->Get() and Store.Header->Set(). - // The method 'Store::Save()' writes everything 'Set' since the last 'Save' to the current entry. - // Store::Clear clears the map of the current entry, to start building a new one. - // Use 'Store::GetEntry(int entrynum)' to load an entry to then be able to 'Get' it's contents. - // 'Store->Header->Get("TotalEntries",NumEvents)' will load the num entries into NumEvents - // ------------------ - // When adding a BoostStore (or class object in general) to a BoostStore (such as ANNIEEvent) - // you call BoostStore::Set("key",ObjectPointer) - BUT be aware that serialization happens when the - // 'Set' method is called - although you pass it a pointer, any subsequent changes to the object - // will NOT get saved! You must call 'Set' AFTER making ALL changes to your object! - ///////////////////////////////////////////////////////////////// - - // Make class private members; e.g. the WCSimT and WCSimRootGeom - // ============================================================= -// file= new TFile(MCFile.c_str(),"READ"); -// wcsimtree= (TTree*) file->Get("wcsimT"); -// WCSimEntry= new wcsimT(wcsimtree); - WCSimEntry= new wcsimT(MCFile.c_str(),verbosity); - - gROOT->cd(); - wcsimrootgeom = WCSimEntry->wcsimrootgeom; - wcsimrootopts = WCSimEntry->wcsimrootopts; - int pretriggerwindow=wcsimrootopts->GetNDigitsPreTriggerWindow(); - int posttriggerwindow=wcsimrootopts->GetNDigitsPostTriggerWindow(); - - // put useful constants into the CStore - // ==================================== - //m_data->CStore.Set("WCSimEntry",WCSimEntry,false); // pass on the WCSim entry - not possible - // pass on root options. TODO should these be saved somewhere? - m_data->CStore.Set("WCSimPreTriggerWindow",pretriggerwindow); - m_data->CStore.Set("WCSimPostTriggerWindow",posttriggerwindow); - // store the WCSimRootGeom for LAPPD reader calculation of global coords - intptr_t geomptr = reinterpret_cast(wcsimrootgeom); - if(verbosity>1) cout<<"wcsimrootgeom at "<Stores.count("WCSimRootGeomStore"); -// if(wcsimgeomexists==0){ -// m_data->Stores["WCSimRootGeomStore"] = new BoostStore(false,0); -// //m_data->Stores.at("WCSimRootGeomStore")->Header->Set("WCSimRootGeom",wcsimrootgeom); -// m_data->Stores.at("WCSimRootGeomStore")->Set("WCSimRootGeom",&wcsimrootgeom); -// } -// -// // Make a WCSimStore to store additional WCSim info passed between tools -// // ===================================================================== -// int wcsimstoreexists = m_data->Stores.count("WCSimStore"); -// cout<<"wcsimstoreexists="<Stores["WCSimStore"] = new BoostStore(false,0); -// } -// m_data->Stores.at("WCSimStore")->Set("WCSimRootGeom",geomptr); -// m_data->Stores.at("WCSimStore")->Set("WCSimPreTriggerWindow",pretriggerwindow); -// m_data->Stores.at("WCSimStore")->Set("WCSimPostTriggerWindow",posttriggerwindow); - - // Make the ANNIEEvent Store if it doesn't exist - // ============================================= - int annieeventexists = m_data->Stores.count("ANNIEEvent"); - if(annieeventexists==0) m_data->Stores["ANNIEEvent"] = new BoostStore(false,2); - - // Convert WCSimRootGeom into ToolChain Geometry class - Geometry* anniegeom = ConstructToolChainGeometry(); - - // Set run-level information in the ANNIEEvent - // =========================================== - /* - At time of writing ANNIEEvent has the following that need to be Set before each Save: - RunNumber - SubrunNumber - EventNumber - MCParticles - RecoParticles - MCHits - TDCData - RawADCData - CalibratedADCData - RawLAPPDData - CalibratedLAPPDData - TriggerData - MCFlag - EventTime - MCEventNum - MCFile - BeamStatus - */ - - EventNumber=0; - MCTriggernum=0; - // pull the first entry to get the MCFile - int nbytesread = WCSimEntry->GetEntry(MCEventNum); // <0 if out of file - if(nbytesread<=0){ - logmessage = "LoadWCSim Tool had no entry "+to_string(MCEventNum); - if(nbytesread==-4){ - logmessage+=": Overran end of TChain! Have you specified more iterations than are available in ToolChainConfig?"; - } else if(nbytesread==0){ - logmessage+=": No TChain loaded! Is your filepath correct?"; - } - Log(logmessage,v_error,verbosity); - cerr<<"############################"<vars.Set("StopLoop",1); - return false; - } - - MCFile = WCSimEntry->GetCurrentFile()->GetName(); - m_data->Stores.at("ANNIEEvent")->Set("MCFile",MCFile); - - if (PMTMask != "None"){ - masked_ids = this->LoadPMTMask(PMTMask); - } - - // use nominal beam values TODO - double beaminten=4.777e+12; - double beampow=3.2545e+16; - uint64_t beamtimestamp=0; - RunStartTime.SetNs(RunStartUser); - beamstat.set_time(TimeClass(beamtimestamp)); - beamstat.set_pot(beampow); - BeamCondition bc = BeamCondition::Ok; - beamstat.set_condition(bc); - - // Construct the other objects we'll be setting at event level, - // pass managed pointers to the ANNIEEvent Store - MCParticles = new std::vector; - MCHits = new std::map>; - TDCData = new std::map>; - EventTime = new TimeClass(); - TriggerClass beamtrigger("beam",5,0,true,0); - TriggerData = new std::vector{beamtrigger}; // FIXME ? one trigger and resetting time is ok? - - // we'll put these in the CStore: so don't delete them in Finalise! It'll get handled by the Store - ParticleId_to_TankTubeIds = new std::map>; - ParticleId_to_MrdTubeIds = new std::map>; - ParticleId_to_VetoTubeIds = new std::map>; - ParticleId_to_TankCharge = new std::map; - ParticleId_to_MrdCharge = new std::map; - ParticleId_to_VetoCharge = new std::map; - trackid_to_mcparticleindex = new std::map; - - //anniegeom->GetChannel(0); // trigger InitChannelMap + file_pmtid.close(); + m_data->CStore.Set("pmt_tubeid_to_channelkey_data", pmtid_to_channelkey); + m_data->CStore.Set("channelkey_to_pmtid_data", channelkey_to_pmtid); + } else { + logmessage = "LoadWCSim::Initialise: PMT ID Configuration file: " + path_chankeymap; + logmessage += " could not be opened! Is the path valid? Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + + // TODO: should not use relative path. + // There should be an env var that points to the top directory + if (!m_variables.Get("ChankeyToMRDIDMap", path_mrd_chankeymap)) { + path_mrd_chankeymap = "./configfiles/LoadWCSim/MRD_Chankey_WCSimID.dat"; + logmessage = "LoadWCSim::Initialise: ChankeyToMRDIDMap not set in config. "; + logmessage += "Using default path: " + path_mrd_chankeymap; + Log(logmessage, v_warning, verbosity); + } + ifstream file_mrdid(path_mrd_chankeymap.c_str()); + if (file_mrdid.is_open()){ + // watch out: comment or empty lines not supported here + while (!file_mrdid.eof()){ + unsigned long chankey; + int mrdid; + file_mrdid >> chankey >> mrdid; + mrdid_to_channelkey.emplace(mrdid,chankey); + } + file_mrdid.close(); + } else { + logmessage = "LoadWCSim::Initialise: MRD ID Configuration file: " + path_mrd_chankeymap; + logmessage += " could not be opened! Is the path valid? Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } + + // TODO: should not use relative path. + // There should be an env var that points to the top directory + if (!m_variables.Get("ChankeyToFMVIDMap", path_fmv_chankeymap)) { + path_fmv_chankeymap = "./configfiles/LoadWCSim/FMV_Chankey_WCSimID.dat"; + logmessage = "LoadWCSim::Initialise: ChankeyToFMVIDMap not set in config. "; + logmessage += "Using default path: " + path_fmv_chankeymap; + Log(logmessage, v_warning, verbosity); + } + ifstream file_fmvid(path_fmv_chankeymap.c_str()); + if (file_fmvid.is_open()){ + // watch out: comment or empty lines not supported here + while (!file_fmvid.eof()){ + unsigned long chankey; + int fmvid; + file_fmvid >> chankey >> fmvid; + fmvid_to_channelkey.emplace(fmvid,chankey); + } + file_fmvid.close(); + } else { + logmessage = "LoadWCSim::Initialise: FMV ID Configuration file: " + path_fmv_chankeymap; + logmessage += " could not be opened! Is the path valid? Aborting!"; + Log(logmessage, v_error, verbosity); + return false; + } - m_data->CStore.Set("UserEvent",false); //enables the ability for other tools to select a specific event number - triggers_event = 0; + // Short Stores README + // ====================================================== + // n.b. m_data->vars is a Store (of ben's Store type) that is not saved to disk. + // m_data->CStore is a single entry binary BoostStore that is not saved to disk. + // m_data->Stores["StoreName"] is a map of binary BoostStores that are saved to disk. + // With BoostStore::Set("MyVariable", myvar), if myvar is not a pointer it will always + // be saved to disk, but if myvar is a pointer the persist flag (default true) controls + // whether it will be saved to disk. In either case the BoostStore becomes the owner + // of the object and will handle its deletion. + // Is 'BoostStore::Save' needed for single-entry stores? + // ---------------- + // BoostStore constructor args: typechecking (bool), m_format (0=binary, 1=ASCII, 2=multievent) + // A BoostStore has a header where useful constants may be saved. The header is a BoostStore itself, + // and can be accessed via: Store.Header->Get() and Store.Header->Set(). + // The method 'BoostStore::Save()' writes everything 'Set' since the last 'Save' to the current entry. + // BoostStore::Clear clears the map of the current entry, to start building a new one. + // Use 'BoostStore::GetEntry(int entrynum)' to load an entry to then be able to 'Get' it's contents. + // 'BoostStore->Header->Get("TotalEntries",NumEvents)' will load the num entries into NumEvents + // ------------------ + // When adding a BoostStore (or class object in general) to a BoostStore (such as ANNIEEvent) + // you call BoostStore::Set("key",ObjectPointer) - BUT be aware that serialization happens when the + // 'Set' method is called - although you pass it a pointer, any subsequent changes to the object + // will NOT get saved! You must call 'Set' AFTER making ALL changes to your object! + ///////////////////////////////////////////////////////// + + // Make class private members; e.g. the WCSimT and WCSimRootGeom + WCSimEntry = new wcsimT(MCFile.c_str(),verbosity); + gROOT->cd(); + + // Make the ANNIEEvent Store if it doesn't exist + // ============================================= + if (m_data->Stores.count("ANNIEEvent") == 0) + m_data->Stores["ANNIEEvent"] = new BoostStore(false, 2); + + MCFile = WCSimEntry->GetCurrentFile()->GetName(); + m_data->Stores.at("ANNIEEvent")->Set("MCFile", MCFile); + + // Grab the geom and record it for LAPPD reader to use later + wcsimrootgeom = WCSimEntry->wcsimrootgeom; + intptr_t geomptr = reinterpret_cast(wcsimrootgeom); + m_data->CStore.Set("WCSimRootGeom", geomptr); + std::ostringstream ss_wcsimrootgeom; + std::ostringstream ss_geomptr; + ss_wcsimrootgeom << wcsimrootgeom; + ss_geomptr << geomptr; + logmessage = "LoadWCSim::Initialise: WCSimRootGeom at: " + ss_wcsimrootgeom.str(); + logmessage += ", and geomptr at: " + ss_geomptr.str(); + Log(logmessage, v_message, verbosity); + + // Convert WCSimRootGeom into ToolChain Geometry class + Geometry* anniegeom = ConstructToolChainGeometry(); + + // Grab root options and record them + wcsimrootopts = WCSimEntry->wcsimrootopts; + int pretriggerwindow = wcsimrootopts->GetNDigitsPreTriggerWindow(); + int posttriggerwindow = wcsimrootopts->GetNDigitsPostTriggerWindow(); + m_data->CStore.Set("WCSimPreTriggerWindow", pretriggerwindow); + m_data->CStore.Set("WCSimPostTriggerWindow", posttriggerwindow); + + // Set run-level information in the ANNIEEvent + // =========================================== + /* + At time of writing ANNIEEvent has the following that need to be Set before each Save: + RunNumber + SubrunNumber + EventNumber + MCParticles + RecoParticles + MCHits + TDCData + RawADCData + CalibratedADCData + RawLAPPDData + CalibratedLAPPDData + TriggerData + MCFlag + EventTime + WCSimEntryNum + MCFile + BeamStatus + */ + + EventNumber = 0; + MCTriggerNum = 0; + + // pull the first entry so it's ready to look at during the first execute loop + int nbytesread = WCSimEntry->GetEntry(WCSimEntryNum); + if (nbytesread <= 0) { + logmessage = "LoadWCSim::Initialise: Issue loading event number: "; + logmessage += std::to_string(WCSimEntryNum) + ". Aborting!"; + Log(logmessage, v_error, verbosity); + + if (nbytesread == -4) { + logmessage = "LoadWCSim::Initialise: Overran the end of TChain! "; + logmessage += "Have you specified more iterations than are available in ToolChainConfig?"; + Log(logmessage, v_error, verbosity); + } + + if (nbytesread==0) { + logmessage = "LoadWCSim::Initialise: No TChain loaded! Is your filepath correct?"; + Log(logmessage, v_error, verbosity); + } + + m_data->vars.Set("StopLoop", 1); + return false; + } + + // TODO: use nominal beam values + double beaminten = 4.777e+12; + double beampow = 3.2545e+16; + uint64_t beamtimestamp = 0; + RunStartTime.SetNs(RunStartUser); + beamstat.set_time(TimeClass(beamtimestamp)); + beamstat.set_pot(beampow); + BeamCondition bc = BeamCondition::Ok; + beamstat.set_condition(bc); + + + // Construct the other objects we'll be setting at event level, + // pass managed pointers to the ANNIEEvent Store + // ============================================= + MCParticles = new std::vector; + MCHits = new std::map>; + TDCData = new std::map>; + EventTime = new TimeClass(); + + // FIXME ? one trigger and resetting time is ok? + // TriggerClass should have a lower case type + std::string triggertype = TriggerType; + std::transform(triggertype.begin(), triggertype.end(), triggertype.begin(), + [](unsigned char c){ return std::tolower(c); }); + + TriggerClass beamtrigger(triggertype, TriggerWord, 0, true, 0); + TriggerData = new std::vector{beamtrigger}; - return true; + // we'll put these in the CStore: so don't delete them in Finalise! It'll get handled by the Store + ParticleId_to_TankTubeIds = new std::map>; + ParticleId_to_MrdTubeIds = new std::map>; + ParticleId_to_VetoTubeIds = new std::map>; + ParticleId_to_TankCharge = new std::map; + ParticleId_to_MrdCharge = new std::map; + ParticleId_to_VetoCharge = new std::map; + trackid_to_mcparticleindex = new std::map; + mapNeutronIsPrim = new std::map; + + // If true then other tools can select a specific event number + m_data->CStore.Set("UserEvent", false); + trigsInEntry = 0; + + return true; } +// Overview of process to extract WCSim files: +// The wcsimT object is a tree that we load into WCSimEntry. +// It contains WCSimRootEvent branches for the tank, MRD, and veto +// (along with geom. and root options). +// Each G4 "event" is obtained by using wcsimT::GetEntry(), +// which loads the specified entry from the WCSimRootEvent branches. +// We loaded the first entry in the initialize stage so we're ready +// at the first execute. We then load the next entry at the end of +// Execute (depending on the splitSubtriggers flag (see below)) +// WCSimRootEvent objects contain "triggers" +// The first trigger contains all the particles initially passed G4. +// Subsequent triggers contain delayed particles (eg. from decays). +// If splitSubtriggers is set to True then each Execute loop will +// add one subtrigger as a single event. +// If splitSubtriggers is set to False then each Execute loop will +// add all subtriggers as a single event. +// In this way the WCSimEntryNum may not match the EventNumber assigned by LoadWCSim. +bool LoadWCSim::Execute() +{ + // Check if another tool has specified an entry and trigger number to load. + // Currently e.g. the EventDisplay has the ability to do that + // ============================================= + bool userEvent; + m_data->CStore.Get("UserEvent", userEvent); + if (userEvent) { + // Grab the user-specified entry + if (m_data->CStore.Get("LoadEvNr", WCSimEntryNum)) { + logmessage = "LoadWCSim::Execute: Going to the user-defined event number: "; + logmessage += std::to_string(WCSimEntryNum); + Log(logmessage, v_message, verbosity); + } else { + logmessage = "LoadWCSim::Execute: Requested to a user-defined event number, "; + logmessage += "but LoadEvNr is not set in the CStore."; + Log(logmessage, v_error, verbosity); + return false; + } + + // Reset the flag and go to the first trigger in that event + m_data->CStore.Set("UserEvent", false); + MCTriggerNum = 0; -bool LoadWCSim::Execute(){ - - // probably not necessary, clears the map for this entry. We're going to re-Set the event entry anyway... - //m_data->Stores.at("ANNIEEvent")->Clear(); - - //check if another tool has specified a specific evnumber to load (currently e.g. the EventDisplay has the ability to do that) - bool user_event; - m_data->CStore.Get("UserEvent",user_event); - if (user_event){ - m_data->CStore.Set("UserEvent",false); - MCTriggernum = 0; //look at first trigger for user-specified event numbers - int user_evnum; - uint16_t currentTriggernum; - bool check_further_triggers=false; - m_data->CStore.Get("LoadEvNr",user_evnum); - m_data->CStore.Get("CheckFurtherTriggers",check_further_triggers); - m_data->CStore.Get("CurrentTriggernum",currentTriggernum); - MCEventNum = user_evnum; - currentTriggernum++; - if (verbosity > 3){ - std::cout <<"check_further_triggers = "<=MaxEntries && MaxEntries>0){ - std::cout<<"LoadWCSim Tool: Reached max entries specified in config file, terminating ToolChain"<vars.Set("StopLoop",1); - } else { - int nbytesread = WCSimEntry->GetEntry(MCEventNum); // <0 if out of file - if (verbosity > v_debug) std::cout <<"LoadWCSim tool: Trying to get next event, MCEventNum: "<vars.Set("StopLoop",1); - return true; - } - } - } - if(verbosity) cout<<"Executing tool LoadWCSim with MC entry "<vars.Get("StopLoop",loopstopped); - if(get_ok && loopstopped){ - // setting StopLoop doesn't terminate the ToolChain if the number of iterations - // is specified manually in the ToolChainConfig. - // This is almost certainly going to result in a segfault somewhere, - // (e.g. if this tool set it in the last loop iteration because it ran out of entries) - // but let's do what we can - Log("WARNING: STOPLOOP HAS BEEN SET. RETURNING",v_error,verbosity); - return 0; - } - MCFile = WCSimEntry->GetCurrentFile()->GetName(); - - MCHits->clear(); - TDCData->clear(); - MCNeutCap.clear(); - MCNeutCapGammas.clear(); - mrd_firstlayer=false; - mrd_lastlayer=false; + // Check if an additional trigger is available and which trigger was requested + bool checkFurtherTriggers = false; + uint16_t currentTriggerNum; + m_data->CStore.Get("CheckFurtherTriggers", checkFurtherTriggers); + m_data->CStore.Get("CurrentTriggernum", currentTriggerNum); + ++currentTriggerNum; - std::map neutcap_is_primary; //map to store whether a neutron capture was from primary neutron or secondary, key = ncapture time, value = was the capture primary? + logmessage = "LoadWCSim::Execute: Check further triggers: " + std::to_string(checkFurtherTriggers); + logmessage += ", current trigger num: " + std::to_string(currentTriggerNum); + logmessage += "\n total number triggers: " + std::to_string(trigsInEntry); + Log(logmessage, v_debug, verbosity); - triggers_event = WCSimEntry->wcsimrootevent->GetNumberOfEvents(); - - int MaxEventNr = MCTriggernum+1; - if (!splitSubtriggers){ - // we'll merge all subtriggers, so loop from 0 to triggers_event; - MCTriggernum = 0; - MaxEventNr = WCSimEntry->wcsimrootevent->GetNumberOfEvents(); - } // else if we're splitting subtriggers, loop from [current MCTriggernum] to [current MCTriggernum+1] - - while (MCTriggernum < MaxEventNr){ - if(verbosity>1) cout<<"getting triggers"<wcsimrootevent->GetTrigger(0); - firsttrigm=WCSimEntry->wcsimrootevent_mrd->GetTrigger(0); - firsttrigv=WCSimEntry->wcsimrootevent_facc->GetTrigger(0); - atrigt = WCSimEntry->wcsimrootevent->GetTrigger(MCTriggernum); - if(MCTriggernum<(WCSimEntry->wcsimrootevent_mrd->GetNumberOfEvents())){ - atrigm = WCSimEntry->wcsimrootevent_mrd->GetTrigger(MCTriggernum); - } else { atrigm=nullptr; } - if(MCTriggernum<(WCSimEntry->wcsimrootevent_facc->GetNumberOfEvents())){ - atrigv = WCSimEntry->wcsimrootevent_facc->GetTrigger(MCTriggernum); - } else { atrigv=nullptr; } - if(verbosity>2) cout<<"wcsimrootevent="<wcsimrootevent<2) cout<<"wcsimrootevent_mrd="<wcsimrootevent_mrd<2) cout<<"wcsimrootevent_facc="<wcsimrootevent_facc<2) cout<<"atrigt="<GetHeader()->GetRun(); - SubrunNumber = 0; - EventTimeNs = atrigt->GetHeader()->GetDate(); - EventTime->SetNs(EventTimeNs); - if(verbosity>2) cout<<"EventTime is "<clear(); - trackid_to_mcparticleindex->clear(); - primarymuonindex=-1; - - std::string geniefilename = firsttrigt->GetHeader()->GetGenieFileName().Data(); - int genieentry = firsttrigt->GetHeader()->GetGenieEntryNum(); - if(verbosity>1) cout<<"Genie file is "<CStore.Set("GenieFile",geniefilename); - m_data->CStore.Set("GenieEntry",genieentry); - - for(int trigi=0; trigiwcsimrootevent->GetNumberOfEvents(); trigi++){ - - WCSimRootTrigger* atrigtt = WCSimEntry->wcsimrootevent->GetTrigger(trigi); - if(verbosity>1) cout<<"getting "<GetNtrack()<<" tracks from trigger "<GetNtrack(); tracki++){ - if(verbosity>2) cout<<"getting track "<GetTracks()->At(tracki); - /* a WCSimRootTrack has methods: - Int_t GetIpnu() pdg - Int_t GetFlag() -1: neutrino primary, -2: neutrino target, 0: other - Float_t GetM() mass - Float_t GetP() momentum magnitude - Float_t GetE() energy (inc rest mass^2) - Float_t GetEndE() energy on stopping of particle tracking - Float_t GetEndP() momentum on stopping of particle tracking - Int_t GetStartvol() starting volume: 10 is tank, 20 is facc, 30 is mrd - Int_t GetStopvol() stopping volume: but these may not be set. - Float_t GetDir(Int_t i=0) momentum unit vector - Float_t GetPdir(Int_t i=0) momentum vector - Float_t GetPdirEnd(Int_t i=0) direction vector on stop tracking - Float_t GetStop(Int_t i=0) stopping vertex x,y,z for i=0-2, in cm - Float_t GetStart(Int_t i=0) starting vertex x,y,z for i=0-2, in cm - Int_t GetParenttype() parent pdg, 0 for primary. - Float_t GetTime() trj->GetGlobalTime(); starting time of particle - Float_t GetStopTime() - Int_t GetId() wcsim trackid - */ - - tracktype startstoptype = tracktype::UNDEFINED; - //MC particle times are relative to the trigger time - if(nextrack->GetFlag()!=0) { - if (nextrack->GetFlag()==-1){ - double starttime, stoptime = -1; - if (splitSubtriggers){ - //MC particle times now stored relative to the trigger time - starttime = (static_cast(nextrack->GetTime()-EventTimeNs)); - stoptime = (static_cast(nextrack->GetStopTime()-EventTimeNs)); - } else { - starttime = (static_cast(nextrack->GetTime())); - stoptime = (static_cast(nextrack->GetStopTime())); - } - MCParticle neutrino( - nextrack->GetIpnu(), nextrack->GetE(), nextrack->GetEndE(), - Position(nextrack->GetStart(0) / 100., - nextrack->GetStart(1) / 100., - nextrack->GetStart(2) / 100.), - Position(nextrack->GetStop(0) / 100., - nextrack->GetStop(1) / 100., - nextrack->GetStop(2) / 100.), - starttime, - stoptime, - Direction(nextrack->GetDir(0), nextrack->GetDir(1), nextrack->GetDir(2)), - (sqrt(pow(nextrack->GetStop(0)-nextrack->GetStart(0),2.)+ - pow(nextrack->GetStop(1)-nextrack->GetStart(1),2.)+ - pow(nextrack->GetStop(2)-nextrack->GetStart(2),2.))) / 100., - startstoptype, - nextrack->GetId(), - nextrack->GetParenttype(), - nextrack->GetFlag(), - trigi); - - //Set the neutrino as its own particle - m_data->Stores["ANNIEEvent"]->Set("NeutrinoParticle",neutrino); - - continue; // flag 0 only is normal particles: excludes neutrino - } - } - double starttime, stoptime = -1; - if (splitSubtriggers){ - //MC particle times now stored relative to the trigger time - starttime = (static_cast(nextrack->GetTime()-EventTimeNs)); - stoptime = (static_cast(nextrack->GetStopTime()-EventTimeNs)); - } else { - starttime = (static_cast(nextrack->GetTime())); - stoptime = (static_cast(nextrack->GetStopTime())); - } - if (verbosity > 2) std::cout <<"LoadWCSim tool, loaded particle with PDG: "<GetIpnu()<<", Time: "<GetParenttype() << ", EndProcess: "<GetEndProcess()<GetIpnu() == 2112 && nextrack->GetEndProcess() == "nCapture"){ - if (verbosity > 2) std::cout <<"LoadWCSim tool: Neutron capture! Parent: "<GetParenttype()<<", Time: "<GetParenttype()==0)); //store whether neutron was primary - } - MCParticle thisparticle( - nextrack->GetIpnu(), nextrack->GetE(), nextrack->GetEndE(), - Position(nextrack->GetStart(0) / 100., - nextrack->GetStart(1) / 100., - nextrack->GetStart(2) / 100.), - Position(nextrack->GetStop(0) / 100., - nextrack->GetStop(1) / 100., - nextrack->GetStop(2) / 100.), - starttime, - stoptime, - Direction(nextrack->GetDir(0), nextrack->GetDir(1), nextrack->GetDir(2)), - (sqrt(pow(nextrack->GetStop(0)-nextrack->GetStart(0),2.)+ - pow(nextrack->GetStop(1)-nextrack->GetStart(1),2.)+ - pow(nextrack->GetStop(2)-nextrack->GetStart(2),2.))) / 100., - startstoptype, - nextrack->GetId(), - nextrack->GetParenttype(), - nextrack->GetFlag(), - trigi); - // not currently in constructor call, but we now have it in latest WCSim files - // XXX this will fall over with older WCSim files, whose WCSimLib doesn't have this method! - thisparticle.SetTankExitPoint(Position(nextrack->GetTankExitPoint(0)/ 100., - nextrack->GetTankExitPoint(1)/ 100., - nextrack->GetTankExitPoint(2)/ 100.)); - if( (nextrack->GetIpnu()==13) && - (nextrack->GetParenttype()==0) && - (nextrack->GetFlag()==0) && - (primarymuonindex<0) ){ - // call this the primary muon. If we have more than one, use the first - primarymuonindex = MCParticles->size(); - } - if((abs(nextrack->GetIpnu())==13)|| - (abs(nextrack->GetIpnu())==211)|| - (nextrack->GetIpnu()==111)){ - if (verbosity > 0) { - std::cout<<"Found "<GetIpnu()<<" with parent pdg " - <GetParenttype()<<", flag "<GetFlag() - <<" track id "<GetId() - << ", start vertex (" + to_string(nextrack->GetStart(0)/100.) - << ", " + to_string(nextrack->GetStart(1)/100.) - << ", " + to_string(nextrack->GetStart(2)/100.) - << "), and end vertex (" + to_string(nextrack->GetStop(0)/100.) - << ", " + to_string(nextrack->GetStop(1)/100.) - << ", " + to_string(nextrack->GetStop(2)/100.) - << ")" - <GetIpnu()==13) && (nextrack->GetParenttype()==0))|| - (MCParticles->size()==0)){ - if (verbosity) std::cout<<"Found "<GetIpnu()<<" with parent pdg " - <GetParenttype()<<", flag "<GetFlag() - <<" track id "<GetId() - <<" at position "<size()<GetIpnu())==211)|| - (nextrack->GetIpnu()==111)){ - if (verbosity) std::cout<<"Found "<GetIpnu()<<" with parent pdg " - <GetParenttype()<<", flag "<GetFlag() - <<" track id "<GetId() - << ", start vertex (" + to_string(nextrack->GetStart(0)/100.) - << ", " + to_string(nextrack->GetStart(1)/100.) - << ", " + to_string(nextrack->GetStart(2)/100.) - << "), and end vertex (" + to_string(nextrack->GetStop(0)/100.) - << ", " + to_string(nextrack->GetStop(1)/100.) - << ", " + to_string(nextrack->GetStop(2)/100.) - << ")" - <GetIpnu())==13)||(nextrack->GetIpnu()==22)) && - ((abs(nextrack->GetParenttype())==211)||(nextrack->GetParenttype()==111))){ - if (verbosity) std::cout<<"Found "<GetIpnu()<<" with parent pdg " - <GetParenttype()<<", flag "<GetFlag() - <<" track id "<GetId() - <<" at position "<size()<GetIpnu()==13){ - logmessage = "Muon found with flag: "+to_string(nextrack->GetFlag()) - + ", parent type " + to_string(nextrack->GetParenttype()) - + ", Id " + to_string(nextrack->GetId()) - + ", start vertex (" + to_string(nextrack->GetStart(0)/100.) - + ", " + to_string(nextrack->GetStart(1)/100.) - + ", " + to_string(nextrack->GetStart(2)/100.) - + "), and end vertex (" + to_string(nextrack->GetStop(0)/100.) - + ", " + to_string(nextrack->GetStop(1)/100.) - + ", " + to_string(nextrack->GetStop(2)/100.) - + ")"; - Log(logmessage,v_debug,verbosity); - } - - trackid_to_mcparticleindex->emplace(nextrack->GetId(),MCParticles->size()); - MCParticles->push_back(thisparticle); - } - if(verbosity>2) cout<<"MCParticles has "<size()<<" entries"<0, since particle times are relative to the trigger time, - // we need to update all the particle times - double timediff = EventTimeNs - firsttrigt->GetHeader()->GetDate(); - for(MCParticle& aparticle : *MCParticles){ - if (splitSubtriggers){ - aparticle.SetStartTime(aparticle.GetStartTime()-timediff); - aparticle.SetStopTime (aparticle.GetStopTime() -timediff); - } - else { - aparticle.SetStartTime(aparticle.GetStartTime()); - aparticle.SetStopTime (aparticle.GetStopTime()); - } - } - } // end updating particle times - - int numtankdigits = atrigt ? atrigt->GetCherenkovDigiHits()->GetEntries() : 0; - if(verbosity>1) cout<<"looping over "<2) cout<<"getting digit "<GetCherenkovDigiHits()->At(digiti); - //WCSimRootChernkovDigiHit has methods GetTubeId(), GetT(), GetQ(), GetPhotonIds() - if(verbosity>2) cout<<"LoadWCSim tool: next digihit at "<GetTubeId(); // geometry TubeID->channelkey map is made INCLUDING offset of 1 - if(pmt_tubeid_to_channelkey.count(tubeid)==0){ - cerr<<"LoadWCSim ERROR: tank PMT with no associated ChannelKey!"<2) cout<<"LoadWCSim tool: tubeid = "<(digihit->GetT()-HistoricTriggeroffset); // relative to trigger - } else { - // instead take the true time of the first photon - std::vector photonids = digihit->GetPhotonIds(); // indices of the digit's photons - double earliestphotontruetime=999999999999; - for(int& aphotonindex : photonids){ - WCSimRootCherenkovHitTime* thehittimeobject = - (WCSimRootCherenkovHitTime*)firsttrigt->GetCherenkovHitTimes()->At(aphotonindex); - if(thehittimeobject==nullptr){ - cerr<<"LoadWCSim Tool: ERROR! Retrieval of photon from digit returned nullptr!"<(thehittimeobject->GetTruetime()); - if(aphotontime2){ cout<<"digittime is "<GetQ(); - if(verbosity>2) cout<<"digit Q is "< parents = GetHitParentIds(digihit, firsttrigt); - if(verbosity>2){ std::cout <<"digittime before adding event time: "< therefore we shouldn't add the event / trigger time if using unsmeared - if(use_smeared_digit_time){ - if (!splitSubtriggers) digittime += EventTimeNs; - } - if(verbosity>2){ std::cout <<"digittime after adding event time: "<count(key)==0) MCHits->emplace(key, std::vector{nexthit}); - else MCHits->at(key).push_back(nexthit); - if(verbosity>2) cout<<"digit added"<2) cout<<"done with tank digits"<GetCherenkovDigiHits()->GetEntries() : 0; - if(verbosity>1) cout<<"adding "<2) cout<<"getting digit "<GetCherenkovDigiHits()->At(digiti); - if(verbosity>2) cout<<"next digihit at "<GetTubeId(); - if(verbosity>2) cout<<"tubeid="<(digihit->GetT()-HistoricTriggeroffset); // relative to trigger - } else { - // instead take the true time of the first photon - std::vector photonids = digihit->GetPhotonIds(); // indices of the digit's photons - double earliestphotontruetime=999999999999; - for(int& aphotonindex : photonids){ - WCSimRootCherenkovHitTime* thehittimeobject = - (WCSimRootCherenkovHitTime*)firsttrigm->GetCherenkovHitTimes()->At(aphotonindex); - if(thehittimeobject==nullptr){ - cerr<<"LoadWCSim Tool: ERROR! Retrieval of photon from digit returned nullptr!"<(thehittimeobject->GetTruetime()); - if(aphotontime2){ cout<<"digittime is "<GetQ(); - if(verbosity>2) cout<<"digit Q is "< parents = GetHitParentIds(digihit, firsttrigm); - if (!splitSubtriggers) digittime += EventTimeNs; - - MCHit nexthit(key, digittime, digiq, parents); - if(TDCData->count(key)==0) {TDCData->emplace(key, std::vector{nexthit}); if (Mrd_Chankey_Layer.at(key)==0) mrd_firstlayer=true; if (Mrd_Chankey_Layer.at(key)==10) mrd_lastlayer=true;} - else TDCData->at(key).push_back(nexthit); - if(verbosity>2) cout<<"digit added"<2) cout<<"done with mrd digits"<GetCherenkovDigiHits()->GetEntries() : 0; - if(verbosity>1) cout<<"adding "<2) cout<<"getting digit "<GetCherenkovDigiHits()->At(digiti); - if(verbosity>2) cout<<"next digihit at "<GetTubeId(); - if(verbosity>2) cout<<"tubeid="<(digihit->GetT()-HistoricTriggeroffset); // relative to trigger - } else { - // instead take the true time of the first photon - std::vector photonids = digihit->GetPhotonIds(); // indices of the digit's photons - double earliestphotontruetime=999999999999; - for(int& aphotonindex : photonids){ - WCSimRootCherenkovHitTime* thehittimeobject = - (WCSimRootCherenkovHitTime*)firsttrigv->GetCherenkovHitTimes()->At(aphotonindex); - if(thehittimeobject==nullptr){ - cerr<<"LoadWCSim Tool: ERROR! Retrieval of photon from digit returned nullptr!"<(thehittimeobject->GetTruetime()); - if(aphotontime2){ cout<<"digittime is "<GetQ(); - if(verbosity>2) cout<<"digit Q is "< parents = GetHitParentIds(digihit, firsttrigv); - digittime += EventTimeNs; - - MCHit nexthit(key, digittime, digiq, parents); - if(TDCData->count(key)==0) TDCData->emplace(key, std::vector{nexthit}); - else TDCData->at(key).push_back(nexthit); - if(verbosity>2) cout<<"digit added"<2) cout<<"done with veto digits"<2) cout<<"setting triggerdata time to "<front().SetTime(EventTimeNs); + // If there is a further trigger in the previous event then we'll load that entry + if (checkFurtherTriggers && currentTriggerNum != trigsInEntry) { + --WCSimEntryNum; + MCTriggerNum = currentTriggerNum; + } + + // Load the desired entry + // If this is the last one in the chain then stop the loop + if ((int)WCSimEntryNum >= MaxEntries && MaxEntries > 0) { + logmessage = "LoadWCSim::Execute: Reached max entries specified in the config file. Terminating ToolChain"; + Log(logmessage, v_error, verbosity); + m_data->vars.Set("StopLoop", 1); + } else { + int nBytesRead = WCSimEntry->GetEntry(WCSimEntryNum); // <0 if out of file + + logmessage = "LoadWCSim::Execute: Trying to get next event, WCSimEntryNum: " + std::to_string(WCSimEntryNum); + logmessage += ", N bytes read: " + std::to_string(nBytesRead); + Log(logmessage, v_debug, verbosity); + + if (nBytesRead <= 0) { + logmessage = "LoadWCSim Tool: Reached last entry of WCSim input file, terminating ToolChain"; + Log(logmessage, v_warning, verbosity); + m_data->vars.Set("StopLoop", 1); + return true; + } + } + }// endif userEvent + + logmessage = "LoadWCSim::Execute: Executing with MC entry: " + std::to_string(WCSimEntryNum); + logmessage += " and trigger: " + std::to_string(MCTriggerNum); + Log(logmessage, v_warning, verbosity); + + // setting StopLoop doesn't terminate the ToolChain if the number of iterations + // is specified manually in the ToolChainConfig. + // This is almost certainly going to result in a segfault somewhere, + // (e.g. if this tool set it in the last loop iteration because it ran out of entries) + // but let's do what we can + int loopStopped = 0; + get_ok = m_data->vars.Get("StopLoop",loopStopped); + if (get_ok && loopStopped){ + logmessage = "LoadWCSim::Execute: WARNING: STOPLOOP HAS BEEN SET. RETURNING"; + Log(logmessage, v_warning, verbosity); + return false; + } + + // Start extracting the things from the WCSim file + // ============================================= + MCFile = WCSimEntry->GetCurrentFile()->GetName(); + + // Clean slate + TDCData->clear(); + MCHits->clear(); + MCNeutCap.clear(); + MCNeutCapGammas.clear(); + //MCHitsToParticles.clear(); + mrd_firstlayer = false; + mrd_lastlayer = false; + + // CherenkovHit(Times) are all in first trigger so just grab them once per execute + WCSimRootTrigger* firstTrigTank = WCSimEntry->wcsimrootevent ->GetTrigger(0); + WCSimRootTrigger* firstTrigMRD = WCSimEntry->wcsimrootevent_mrd ->GetTrigger(0); + WCSimRootTrigger* firstTrigVeto = WCSimEntry->wcsimrootevent_facc->GetTrigger(0); + + // Grab GENIE info and record it + std::string genieFilename = firstTrigTank->GetHeader()->GetGenieFileName().Data(); + int genieEntry = firstTrigTank->GetHeader()->GetGenieEntryNum(); + m_data->CStore.Set("GenieFile",genieFilename); + m_data->CStore.Set("GenieEntry",genieEntry); + + logmessage = "LoadWCSim::Execute: GENIE file is " + genieFilename; + logmessage += ", GENIE evt num is: " + std::to_string(genieEntry); + Log(logmessage, v_message, verbosity); + + + // How many triggers to loop over in this Execute? + // If splitSubtriggers is True then we will "loop" over a single trigger + // defined by previously set MCTriggerNum + // If splitSubtriggers is False then we will loop from 0 up to the total + // number of triggers + trigsInEntry = WCSimEntry->wcsimrootevent->GetNumberOfEvents(); + MCTriggerNum = splitSubtriggers ? MCTriggerNum : 0; + int MaxEventNr = splitSubtriggers ? MCTriggerNum + 1 : trigsInEntry; + + logmessage = "LoadWCSim::Execute: There are " + std::to_string(trigsInEntry); + logmessage += " triggers in this entry"; + Log(logmessage, v_warning, verbosity); + + int nMRDTriggers = WCSimEntry->wcsimrootevent_mrd->GetNumberOfEvents(); + int nVetoTriggers = WCSimEntry->wcsimrootevent_facc->GetNumberOfEvents(); + + // Loop over over the triggers + // ============================================= + while (MCTriggerNum < MaxEventNr) { + logmessage = "LoadWCSim::Execute: Getting trigger " + std::to_string(MCTriggerNum); + logmessage += " of " + std::to_string(trigsInEntry); + Log(logmessage, v_message, verbosity); + + WCSimRootTrigger* aTrigTank = WCSimEntry->wcsimrootevent->GetTrigger(MCTriggerNum); + WCSimRootTrigger* aTrigMRD = ( (MCTriggerNum < nMRDTriggers) + ? WCSimEntry->wcsimrootevent_mrd->GetTrigger(MCTriggerNum) + : nullptr ); + WCSimRootTrigger* aTrigVeto = ( (MCTriggerNum < nVetoTriggers) + ? WCSimEntry->wcsimrootevent_facc->GetTrigger(MCTriggerNum) + : nullptr ); + + std::ostringstream ssEventTank, ssEventMRD, ssEventVeto; + std::ostringstream ssTrigTank, ssTrigMRD, ssTrigVeto; + ssEventTank << WCSimEntry->wcsimrootevent; + ssEventMRD << WCSimEntry->wcsimrootevent_mrd; + ssEventVeto << WCSimEntry->wcsimrootevent_facc; + ssTrigTank << aTrigTank; ssTrigMRD << aTrigMRD; ssTrigVeto << aTrigVeto; + logmessage = "LoadWCSim::Execute: memory locations\n"; + logmessage += " wcsimrootevent: " + ssEventTank.str(); + logmessage += ", wcsimrootevent_mrd: " + ssEventMRD.str(); + logmessage += ", wcsimrootevent_facc: " + ssEventVeto.str() + "\n"; + logmessage += " aTrigTank: " + ssTrigTank.str(); + logmessage += ", aTrigMRD: " + ssTrigMRD.str(); + logmessage += ", aTrigVeto: " + ssTrigVeto.str(); + Log(logmessage, v_debug, verbosity); + + + logmessage = "LoadWCSim::Execute: Getting event date"; + Log(logmessage, v_message, verbosity); + + RunNumber = aTrigTank->GetHeader()->GetRun(); + SubrunNumber = 0; + EventTimeNs = aTrigTank->GetHeader()->GetDate(); + EventTime->SetNs(EventTimeNs); + + logmessage = "LoadWCSim::Execute: EventTime is: " + std::to_string(EventTimeNs); + Log(logmessage, v_debug, verbosity); + + // Load everything + // LoadHits returns a bool that we could use, but don't + // ============================================= + LoadMCParticles(firstTrigTank); + LoadNeutronCaptures(aTrigTank); + LoadHits(aTrigTank, firstTrigTank, "Tank"); + LoadHits(aTrigMRD, firstTrigMRD, "MRD" ); + LoadHits(aTrigVeto, firstTrigVeto, "Veto"); + + // Set the event time + logmessage = "LoadWCSim::Execute: Setting triggerdata time to " + std::to_string(EventTimeNs) + " ns"; + Log(logmessage, v_message, verbosity); + TriggerData->front().SetTime(EventTimeNs); - //Load neutron capture information - int numcaptures = atrigt ? atrigt->GetNcaptures() : 0; - if(verbosity>1) cout<<"LoadWCSim tool: looping over "<2) cout<<"getting capture # "<GetCaptures()->At(capi); - - int capt_parent = capt->GetCaptureParent(); - double capt_vtxx = capt->GetCaptureVtx(0); - double capt_vtxy = capt->GetCaptureVtx(1); - double capt_vtxz = capt->GetCaptureVtx(2); - int capt_ngamma = capt->GetNGamma(); - double capt_totalE = capt->GetTotalGammaE(); - double capt_t = capt->GetCaptureT(); - if (splitSubtriggers){ - capt_t -= EventTimeNs; - } - int capt_nucleus = capt->GetCaptureNucleus(); - double double_primary = -9999; - if (neutcap_is_primary.count(capt_t) > 0){ - bool is_primary = neutcap_is_primary.at(capt_t); - if (verbosity > 2) std::cout <<"LoadWCSim tool: NeutCap object at "< check WCSim options - Log("LoadWCSim tool: Capture parent: "+std::to_string(capt_parent)+", capt_nucleus: "+std::to_string(capt_nucleus)+ ", time " + std::to_string(capt_t) + " was not found in neutcap_is_primary map. Was the EndProcess saved in the WCSim output file?",v_warning,verbosity); - } - std::vector gamma_energies; - for (int i_gamma=0; i_gamma < capt_ngamma; i_gamma++){ - WCSimRootCaptureGamma* captgamma = (WCSimRootCaptureGamma*) capt->GetGammas()->At(i_gamma); - gamma_energies.push_back(captgamma->GetE()); - } - if (MCNeutCap.size()==0){ - MCNeutCap.emplace("CaptParent",std::vector{double(capt_parent)}); - MCNeutCap.emplace("CaptVtxX",std::vector{capt_vtxx}); - MCNeutCap.emplace("CaptVtxY",std::vector{capt_vtxy}); - MCNeutCap.emplace("CaptVtxZ",std::vector{capt_vtxz}); - MCNeutCap.emplace("CaptNGamma",std::vector{double(capt_ngamma)}); - MCNeutCap.emplace("CaptTotalE",std::vector{capt_totalE}); - MCNeutCap.emplace("CaptTime",std::vector{capt_t}); - MCNeutCap.emplace("CaptNucleus",std::vector{double(capt_nucleus)}); - MCNeutCapGammas.emplace("CaptGammas",std::vector>{gamma_energies}); - MCNeutCap.emplace("CaptPrimary",std::vector{double_primary}); - } else { - MCNeutCap.at("CaptParent").push_back(capt_parent); - MCNeutCap.at("CaptVtxX").push_back(capt_vtxx); - MCNeutCap.at("CaptVtxY").push_back(capt_vtxy); - MCNeutCap.at("CaptVtxZ").push_back(capt_vtxz); - MCNeutCap.at("CaptNGamma").push_back(capt_ngamma); - MCNeutCap.at("CaptTotalE").push_back(capt_totalE); - MCNeutCap.at("CaptTime").push_back(capt_t); - MCNeutCap.at("CaptNucleus").push_back(capt_nucleus); - MCNeutCapGammas.at("CaptGammas").push_back(gamma_energies); - MCNeutCap.at("CaptPrimary").push_back(double_primary); - } - } - // update the information about tracks and which tank/mrd/veto PMTs they hit - // this needs updating with each MC trigger, as digits are grouped into MC trigger - // so these maps will then only contain the digits respective particles create - // in the active trigger - // - // ParticleId_to_TankTubeIds is a std::map> - // where TotalCharge is the total charge from that particle on that tube - // (in the event that the particle generated several hits on the tube) - // ParticleId_to_TankCharge is a std::map - // where TotalCharge is summed over all digits, on all pmts, which contained - // light from that particle - ParticleId_to_TankTubeIds->clear(); - ParticleId_to_MrdTubeIds->clear(); - ParticleId_to_VetoTubeIds->clear(); - ParticleId_to_TankCharge->clear(); - ParticleId_to_MrdCharge->clear(); - ParticleId_to_VetoCharge->clear(); - MakeParticleToPmtMap(atrigt, firsttrigt, ParticleId_to_TankTubeIds, ParticleId_to_TankCharge, pmt_tubeid_to_channelkey); - MakeParticleToPmtMap(atrigm, firsttrigm, ParticleId_to_MrdTubeIds, ParticleId_to_MrdCharge, mrd_tubeid_to_channelkey); - MakeParticleToPmtMap(atrigv, firsttrigv, ParticleId_to_VetoTubeIds, ParticleId_to_VetoCharge, facc_tubeid_to_channelkey); - MCTriggernum++; - } //End of MCTriggerNum loop - - //int mrdentries; - //m_data->Stores.at("TDCData")->Get("TotalEntries",mrdentries); // ?? -// m_data->Stores("WCSimEntries")->Set("wcsimrootevent",WCSimEntry->wcsimrootevent); -// m_data->Stores("WCSimEntries")->Set("wcsimrootevent_mrd",WCSimEntry->wcsimrootevent_mrd); -// m_data->Stores("WCSimEntries")->Set("wcsimrootevent_facc",*(WCSimEntry->wcsimrootevent_facc)); - - - // set event level variables - if(verbosity>1) cout<<"setting the store variables"<Stores.at("ANNIEEvent")->Set("RunNumber",RunNumber); - m_data->Stores.at("ANNIEEvent")->Set("SubrunNumber",SubrunNumber); - m_data->Stores.at("ANNIEEvent")->Set("RunType",RunType); - m_data->Stores.at("ANNIEEvent")->Set("EventNumber",EventNumber); - if(verbosity>2) cout<<"particles"<Stores.at("ANNIEEvent")->Set("MCParticles",MCParticles,true); - //Set up Particles object for reconstructed particles - std::vector Particles_Reco; - m_data->Stores["ANNIEEvent"]->Set("Particles",Particles_Reco); - if(verbosity>2) cout<<"hits"<Stores.at("ANNIEEvent")->Set("MCHits",MCHits,true); - if(verbosity>2) cout<<"tdcdata"<Stores.at("ANNIEEvent")->Set("TDCData",TDCData,true); - // TODO? - // right now we have three Time variables: - // 1. "RunStartTime" which stores just a user passed start time. - // 2. "TriggerData", a one-element vector of TriggerClass objects, each of which has a time member. - // The one used entry has its time member set to the time between MC event start (always 0) - // and the MC trigger time. - // 3. "EventTime" which stores this same time difference between MC Trigger and MC event start - // This means we simulate many events all happening ~ the unix epoch (time 0). - // If we want to simulate a 'run' of events, we need to throw some cumulative running time and - // add this to RunStartTime, and store the Event and Trigger times separately. - // This is done by PulseSimulation tool (timefileout_Time), but maybe should be here... - if(verbosity>2) cout<<"triggerdata"<Stores.at("ANNIEEvent")->Set("TriggerData",TriggerData,true); - if(verbosity>2) cout<<"eventtime"<Stores.at("ANNIEEvent")->Set("RunStartTime",runstarttime); - m_data->Stores.at("ANNIEEvent")->Set("EventTimeTank",runstarttime); - m_data->Stores.at("ANNIEEvent")->Set("EventTimeMRD",RunStartTime); - m_data->Stores.at("ANNIEEvent")->Set("EventTime",EventTime,true); - m_data->Stores.at("ANNIEEvent")->Set("MCEventNum",MCEventNum); - // if we merged the subtriggers, report an MCTriggernum of 0 - // so downstream tools know it's a new event - int reported_triggernum = splitSubtriggers ? MCTriggernum-1 : 0; - m_data->Stores.at("ANNIEEvent")->Set("MCTriggernum",reported_triggernum); - m_data->Stores.at("ANNIEEvent")->Set("MCFile",MCFile); - m_data->Stores.at("ANNIEEvent")->Set("MCFlag",true); - //m_data->Stores.at("ANNIEEvent")->Set("BeamStatus",BeamStatus,true); - m_data->Stores.at("ANNIEEvent")->Set("BeamStatus",beamstat); - m_data->Stores.at("ANNIEEvent")->Set("MCNeutCap",MCNeutCap); - m_data->Stores.at("ANNIEEvent")->Set("MCNeutCapGammas",MCNeutCapGammas); - m_data->CStore.Set("NumTriggersThisMCEvt",WCSimEntry->wcsimrootevent->GetNumberOfEvents()); - // auxilliary information about MC Truth particles - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_TankTubeIds", ParticleId_to_TankTubeIds, false); - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_TankCharge", ParticleId_to_TankCharge, false); - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_MrdTubeIds", ParticleId_to_MrdTubeIds, false); - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_MrdCharge", ParticleId_to_MrdCharge, false); - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_VetoTubeIds", ParticleId_to_VetoTubeIds, false); - m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_VetoCharge", ParticleId_to_VetoCharge, false); - m_data->Stores.at("ANNIEEvent")->Set("TrackId_to_MCParticleIndex",trackid_to_mcparticleindex,false); - //Change MRD Triggertype if first and last layer saw a hit (Hardware Cosmic Trigger) - if (mrd_lastlayer && mrd_firstlayer) Triggertype = "Cosmic"; - m_data->Stores.at("ANNIEEvent")->Set("MRDTriggerType",Triggertype); - m_data->Stores.at("ANNIEEvent")->Set("PrimaryMuonIndex",primarymuonindex); - m_data->Stores.at("ANNIEEvent")->Set("TriggerWord",TriggerWord); - - std::map DataStreams; - DataStreams.emplace(std::make_pair("Tank",true)); - DataStreams.emplace(std::make_pair("MRD",true)); - DataStreams.emplace(std::make_pair("Trigger",true)); - DataStreams.emplace(std::make_pair("LAPPD",true)); - m_data->Stores.at("ANNIEEvent")->Set("DataStreams",DataStreams); + // update the information about tracks and which tank/mrd/veto PMTs they hit + // this needs updating with each MC trigger, as digits are grouped into MC trigger + // so these maps will then only contain the digits respective particles create + // in the active trigger + // + // ParticleId_to_TankTubeIds is a std::map> + // where TotalCharge is the total charge from that particle on that tube + // (in the event that the particle generated several hits on the tube) + // ParticleId_to_TankCharge is a std::map + // where TotalCharge is summed over all digits, on all pmts, which contained + // light from that particle + ParticleId_to_TankTubeIds->clear(); + ParticleId_to_MrdTubeIds->clear(); + ParticleId_to_VetoTubeIds->clear(); + ParticleId_to_TankCharge->clear(); + ParticleId_to_MrdCharge->clear(); + ParticleId_to_VetoCharge->clear(); + MakeParticleToPmtMap(aTrigTank, firstTrigTank, ParticleId_to_TankTubeIds, + ParticleId_to_TankCharge, pmt_tubeid_to_channelkey); + MakeParticleToPmtMap(aTrigMRD, firstTrigMRD, ParticleId_to_MrdTubeIds, + ParticleId_to_MrdCharge, mrd_tubeid_to_channelkey); + MakeParticleToPmtMap(aTrigVeto, firstTrigVeto, ParticleId_to_VetoTubeIds, + ParticleId_to_VetoCharge, facc_tubeid_to_channelkey); + MCTriggerNum++; + } //End of MCTriggerNum loop - int TriggerExtended = 1; //1: We have an extended readout for all MC events - m_data->Stores.at("ANNIEEvent")->Set("TriggerExtended",TriggerExtended); - //Things that need to be set by later tools: - //RawADCData - //CalibratedADCData - //RawLAPPDData - //CalibratedLAPPDData - //RecoParticles - if(verbosity>1) cout<<"done loading event"<Stores.at("ANNIEEvent")->Set("RunNumber", RunNumber); + m_data->Stores.at("ANNIEEvent")->Set("SubrunNumber", SubrunNumber); + m_data->Stores.at("ANNIEEvent")->Set("RunType", RunType); + m_data->Stores.at("ANNIEEvent")->Set("EventNumber", EventNumber); + + logmessage = "LoadWCSim::Execute: Setting the MCParticles"; + Log(logmessage, v_debug, verbosity); + m_data->Stores.at("ANNIEEvent")->Set("MCParticles", MCParticles, true); + + // Set up Particles object for reconstructed particles + std::vector Particles_Reco; + m_data->Stores.at("ANNIEEvent")->Set("Particles", Particles_Reco); + + logmessage = "LoadWCSim::Execute: Setting the MCHits"; + Log(logmessage, v_debug, verbosity); + m_data->Stores.at("ANNIEEvent")->Set("MCHits", MCHits, true); + + logmessage = "LoadWCSim::Execute: Setting the TDCData"; + Log(logmessage, v_debug, verbosity); + m_data->Stores.at("ANNIEEvent")->Set("TDCData", TDCData, true); + + // TODO? + // right now we have three Time variables: + // 1. "RunStartTime" which stores just a user passed start time. + // 2. "TriggerData", a one-element vector of TriggerClass objects, each of which has a time member. + // The one used entry has its time member set to the time between MC event start (always 0) + // and the MC trigger time. + // 3. "EventTime" which stores this same time difference between MC Trigger and MC event start + // This means we simulate many events all happening ~ the unix epoch (time 0). + // If we want to simulate a 'run' of events, we need to throw some cumulative running time and + // add this to RunStartTime, and store the Event and Trigger times separately. + // This is done by PulseSimulation tool (timefileout_Time), but maybe should be here... + logmessage = "LoadWCSim::Execute: Setting the TriggerData"; + Log(logmessage, v_debug, verbosity); + m_data->Stores.at("ANNIEEvent")->Set("TriggerData", TriggerData, true); + + logmessage = "LoadWCSim::Execute: Setting the EventTime"; + Log(logmessage, v_debug, verbosity); + long runstarttime = RunStartTime.GetNs(); + m_data->Stores.at("ANNIEEvent")->Set("RunStartTime", runstarttime); + m_data->Stores.at("ANNIEEvent")->Set("EventTimeTank", runstarttime); + m_data->Stores.at("ANNIEEvent")->Set("EventTimeMRD", RunStartTime); + m_data->Stores.at("ANNIEEvent")->Set("EventTime", EventTime, true); + m_data->Stores.at("ANNIEEvent")->Set("MCEventNum", WCSimEntryNum); + + // If we merged the subtriggers, report an MCTriggerNum of 0 + // so downstream tools know it's a new event + int reported_triggernum = splitSubtriggers ? MCTriggerNum-1 : 0; + m_data->Stores.at("ANNIEEvent")->Set("MCTriggerNum", reported_triggernum); + m_data->Stores.at("ANNIEEvent")->Set("MCFile", MCFile); + m_data->Stores.at("ANNIEEvent")->Set("MCFlag", true); + m_data->Stores.at("ANNIEEvent")->Set("BeamStatus", beamstat); + m_data->Stores.at("ANNIEEvent")->Set("MCNeutCap", MCNeutCap); + m_data->Stores.at("ANNIEEvent")->Set("MCNeutCapGammas", MCNeutCapGammas); + m_data->CStore.Set("NumTriggersThisMCEvt", trigsInEntry); + + // auxilliary information about MC Truth particles + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_TankTubeIds", ParticleId_to_TankTubeIds, false); + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_TankCharge", ParticleId_to_TankCharge, false); + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_MrdTubeIds", ParticleId_to_MrdTubeIds, false); + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_MrdCharge", ParticleId_to_MrdCharge, false); + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_VetoTubeIds", ParticleId_to_VetoTubeIds, false); + m_data->Stores.at("ANNIEEvent")->Set("ParticleId_to_VetoCharge", ParticleId_to_VetoCharge, false); + m_data->Stores.at("ANNIEEvent")->Set("TrackId_to_MCParticleIndex",trackid_to_mcparticleindex,false); + + //Change MRD TriggerType if first and last layer saw a hit (Hardware Cosmic Trigger) + if (mrd_lastlayer && mrd_firstlayer) TriggerType = "Cosmic"; + m_data->Stores.at("ANNIEEvent")->Set("MRDTriggerType", TriggerType); + m_data->Stores.at("ANNIEEvent")->Set("PrimaryMuonIndex", primaryMuonIndex); + m_data->Stores.at("ANNIEEvent")->Set("TriggerWord", TriggerWord); - EventNumber++; - if(verbosity>2) cout<<"checking if we're done on trigs in this event"<wcsimrootevent->GetNumberOfEvents()){ - MCTriggernum=0; - MCEventNum++; - newentry=true; - if(verbosity>2) cout<<"this is the last trigger in the event: next loop will process a new event"<2) cout<<"there are further triggers in this event: next loop will process the trigger "<wcsimrootevent->GetNumberOfEvents()<=MaxEntries && MaxEntries>0){ - cout<<"LoadWCSim Tool: Reached max entries specified in config file, terminating ToolChain"<vars.Set("StopLoop",1); - } else { - int nbytesread = WCSimEntry->GetEntry(MCEventNum); // <0 if out of file - if (verbosity > 2) std::cout <<"LoadWCSim tool: Trying to get next event, MCEventNum: "<vars.Set("StopLoop",1); - } - } - } + std::map DataStreams; + DataStreams.emplace(std::make_pair("Tank", true)); + DataStreams.emplace(std::make_pair("MRD", true)); + DataStreams.emplace(std::make_pair("Trigger", true)); + DataStreams.emplace(std::make_pair("LAPPD", true)); + m_data->Stores.at("ANNIEEvent")->Set("DataStreams", DataStreams); + + // We have an extended readout for all MC events, set to 1 + m_data->Stores.at("ANNIEEvent")->Set("TriggerExtended", 1); + + logmessage = "LoadWCSim::Execute: Done loading event"; + Log(logmessage, v_message, verbosity); + + ++EventNumber; + + // Things that will be set by later tools: + // RawADCData + // CalibratedADCData + // RawLAPPDData + // CalibratedLAPPDData + // RecoParticles + + logmessage = "LoadWCSim::Execute: Checking if we're done with trigs from this event"; + Log(logmessage, v_debug, verbosity); + + bool newentry = false; + if (MCTriggerNum == trigsInEntry) { + MCTriggerNum = 0; + ++WCSimEntryNum; + newentry = true; + + logmessage = "LoadWCSim::Execute: This is the last trigger in the event. "; + logmessage += "Next loop with process a new event"; + Log(logmessage, v_debug, verbosity); + } else { + logmessage = "LoadWCSim::Execute: There are additional trigger in this event. "; + logmessage += "Next loop with process trigger: " + std::to_string(MCTriggerNum); + logmessage += "/" + std::to_string(trigsInEntry); + + } + + // Load the next entry to get it ready for the next Execute loop + // If this is the last one in the chain then stop the loop + if (newentry) { + // if next loop is processing the next trigger in the same entry, no need to re-load it + if ((int)WCSimEntryNum >= MaxEntries && MaxEntries > 0) { + logmessage = "LoadWCSim::Execute: Reached the max number of entries specified in the config. "; + logmessage += "Terminating the ToolChain"; + Log(logmessage, v_warning, verbosity); + + m_data->vars.Set("StopLoop", 1); + } else { + int nBytesRead = WCSimEntry->GetEntry(WCSimEntryNum); // <0 if out of file + + logmessage = "LoadWCSim::Execute: Trying to get the next event. "; + logmessage += "WCSimEntryNum: " + std::to_string(WCSimEntryNum); + logmessage += ", nBytesRead: " + std::to_string(nBytesRead); + Log(logmessage, v_debug, verbosity); + + + // Check if we've reached the end of the file + if (nBytesRead <= 0){ + logmessage = "LoadWCSim::Execute: Reached the last entry of the WCSim input file. "; + logmessage += "Terminating the ToolChain"; + Log(logmessage, v_warning, verbosity); + + m_data->vars.Set("StopLoop", 1); + } + } + } - //gObjectTable->Print(); - return true; + return true; } -bool LoadWCSim::Finalise(){ - WCSimEntry->GetCurrentFile()->Close(); - delete WCSimEntry; - - // any pointers put in Stores to objects we do not want the Store to clean up - // must be nullified before in finalise to prevent double free - // can't just put 0 or nullptr directly as type must be recognisable as a pointer -// std::map>* ParticleId_to_TankTubeIds_nullptr = nullptr; -// std::map>* ParticleId_to_MrdTubeIds_nullptr = nullptr; -// std::map>* ParticleId_to_VetoTubeIds_nullptr = nullptr; -// std::map* ParticleId_to_TankCharge_nullptr = nullptr; -// std::map* ParticleId_to_MrdCharge_nullptr = nullptr; -// std::map* ParticleId_to_VetoCharge_nullptr = nullptr; +bool LoadWCSim::Finalise() +{ + WCSimEntry->GetCurrentFile()->Close(); + delete WCSimEntry; -// m_data->CStore.Set("ParticleId_to_TankTubeIds", ParticleId_to_TankTubeIds_nullptr, false); -// m_data->CStore.Set("ParticleId_to_TankCharge", ParticleId_to_TankCharge_nullptr, false); -// m_data->CStore.Set("ParticleId_to_MrdTubeIds", ParticleId_to_MrdTubeIds_nullptr, false); -// m_data->CStore.Set("ParticleId_to_MrdCharge", ParticleId_to_MrdCharge_nullptr, false); -// m_data->CStore.Set("ParticleId_to_VetoTubeIds", ParticleId_to_VetoTubeIds_nullptr, false); -// m_data->CStore.Set("ParticleId_to_VetoCharge", ParticleId_to_VetoCharge_nullptr, false); - - return true; + return true; } -Geometry* LoadWCSim::ConstructToolChainGeometry(){ - // Pull details from the WCSimRootGeom - // =================================== - double WCSimGeometryVer = 1; // TODO: pull this from some suitable variable - // some variables are available from wcsimrootgeom. TODO: put everything into wcsimrootgeom - numtankpmts = wcsimrootgeom->GetWCNumPMT(); - numlappds = wcsimrootgeom->GetWCNumLAPPD(); - nummrdpmts = wcsimrootgeom->GetWCNumMRDPMT(); - numvetopmts = wcsimrootgeom->GetWCNumFACCPMT(); - double tank_xcentre = (wcsimrootgeom->GetWCOffset(0)) / 100.; // convert [cm] to [m] - double tank_ycentre = (wcsimrootgeom->GetWCOffset(1)) / 100.; - double tank_zcentre = (wcsimrootgeom->GetWCOffset(2)) / 100.; - Position tank_centre(tank_xcentre, tank_ycentre, tank_zcentre); - double tank_radius = (wcsimrootgeom->GetWCCylRadius()) / 100.; //GetWCCylRadius() returns the black sheet radius not the tank radius - double tank_halfheight = (wcsimrootgeom->GetWCCylLength()) / 200.; //GetWCCylLength() returns the main annulus height not the tank height - //Currently hard-coded; estimated with a tape measure on the ANNIE frame :) - double pmt_enclosed_radius = 1.0; - double pmt_enclosed_halfheight = 1.45; - // geometry variables not yet in wcsimrootgeom are in MRDSpecs.hh - double mrd_width = (MRDSpecs::MRD_width) / 100.; - double mrd_height = (MRDSpecs::MRD_height) / 100.; - double mrd_depth = (MRDSpecs::MRD_depth) / 100.; - double mrd_start = (MRDSpecs::MRD_start) / 100.; - if(verbosity>1) cout<<"we have "<1){ - cout<<"constructed anniegeom at "<Stores.at("ANNIEEvent")->Header->Set("AnnieGeometry",anniegeom,true); +Geometry* LoadWCSim::ConstructToolChainGeometry() +{ + // Pull details from the WCSimRootGeom + // ====================================================== + double WCSimGeometryVer = 1; // TODO: pull this from some suitable variable + // some variables are available from wcsimrootgeom. TODO: put everything into wcsimrootgeom + numtankpmts = wcsimrootgeom->GetWCNumPMT(); + numlappds = wcsimrootgeom->GetWCNumLAPPD(); + nummrdpmts = wcsimrootgeom->GetWCNumMRDPMT(); + numvetopmts = wcsimrootgeom->GetWCNumFACCPMT(); + + Position tank_centre(wcsimrootgeom->GetWCOffset(0), + wcsimrootgeom->GetWCOffset(1), + wcsimrootgeom->GetWCOffset(2)); + tank_centre.UnitToMeter(); + + // GetWCCylRadius() returns the black sheet radius not the tank radius + double tank_radius = (wcsimrootgeom->GetWCCylRadius()) / 100.; + + // GetWCCylLength() returns the main annulus height not the tank height + double tank_halfheight = 0.5*(wcsimrootgeom->GetWCCylLength()) / 100.; + + //Currently hard-coded; estimated with a tape measure on the ANNIE frame :) + double pmt_enclosed_radius = 1.0; + double pmt_enclosed_halfheight = 1.45; + + // geometry variables not yet in wcsimrootgeom are in MRDSpecs.hh + double mrd_width = (MRDSpecs::MRD_width) / 100.; + double mrd_height = (MRDSpecs::MRD_height) / 100.; + double mrd_depth = (MRDSpecs::MRD_depth) / 100.; + double mrd_start = (MRDSpecs::MRD_start) / 100.; + + logmessage = "LoadWCSim::ConstructToolChainGeometry: We have "; + logmessage += std::to_string(numtankpmts) + " tank PMTs, "; + logmessage += std::to_string(nummrdpmts) + " MRD PMTs, and"; + logmessage += std::to_string(numlappds) + " LAPPDs"; + Log(logmessage, v_message, verbosity); - // Construct the Detectors and Channels - // ==================================== - // PMTs - unsigned int ADC_Crate_Num = 0; - unsigned int ADC_Card_Num = 0; - unsigned int ADC_Chan_Num = 0; - unsigned int MT_Crate_Num = 0; - unsigned int MT_Card_Num = 0; - unsigned int MT_Chan_Num = 0; - // LAPPDs - unsigned int ACDC_Crate_Num = 0; - unsigned int ACDC_Card_Num = 0; - unsigned int ACDC_Chan_Num = 0; - unsigned int ACC_Crate_Num = 0; - unsigned int ACC_Card_Num = 0; - unsigned int ACC_Chan_Num = 0; - // TDCs - unsigned int TDC_Crate_Num = 0; - unsigned int TDC_Card_Num = 0; - unsigned int TDC_Chan_Num = 0; - // HV - unsigned int CAEN_HV_Crate_Num = 0; - unsigned int CAEN_HV_Card_Num = 0; - unsigned int CAEN_HV_Chan_Num = 0; - unsigned int LeCroy_HV_Crate_Num = 0; - unsigned int LeCroy_HV_Card_Num = 0; - unsigned int LeCroy_HV_Chan_Num = 0; - unsigned int LAPPD_HV_Crate_Num = 0; - unsigned int LAPPD_HV_Card_Num = 0; - unsigned int LAPPD_HV_Chan_Num = 0; + // construct the ToolChain Goemetry + // ====================================================== + auto* anniegeom = new Geometry(WCSimGeometryVer, + tank_centre, + tank_radius, + tank_halfheight, + pmt_enclosed_radius, + pmt_enclosed_halfheight, + mrd_width, + mrd_height, + mrd_depth, + mrd_start, + numtankpmts, + nummrdpmts, + numvetopmts, + numlappds, + geostatus::FULLY_OPERATIONAL); + + logmessage = "LoadWCSim::ConstructToolChainGeometry: Constructed anniegeom "; + logmessage += "with tank origin (x, y, z): (" + std::to_string(tank_centre.X()); + logmessage += ", " + std::to_string(tank_centre.Y()); + logmessage += ", " + std::to_string(tank_centre.Z()) + ")"; + Log(logmessage, v_message, verbosity); + + m_data->Stores.at("ANNIEEvent")->Header->Set("AnnieGeometry", anniegeom, true); + // Construct the Detectors and Channels + // ====================================================== + ConstructDetectors(anniegeom, numtankpmts, "Tank" ); + ConstructDetectors(anniegeom, nummrdpmts, "MRD" ); + ConstructDetectors(anniegeom, numvetopmts, "Veto" ); + ConstructDetectors(anniegeom, numlappds, "LAPPD"); + + // for other WCSim tools that may need the WCSim Tube IDs + m_data->CStore.Set("lappd_tubeid_to_detectorkey", lappd_tubeid_to_detectorkey); + m_data->CStore.Set("pmt_tubeid_to_channelkey", pmt_tubeid_to_channelkey ); + m_data->CStore.Set("mrd_tubeid_to_channelkey", mrd_tubeid_to_channelkey ); + m_data->CStore.Set("facc_tubeid_to_channelkey", facc_tubeid_to_channelkey ); + // inverse + m_data->CStore.Set("detectorkey_to_lappdid", detectorkey_to_lappdid ); + m_data->CStore.Set("channelkey_to_pmtid", channelkey_to_pmtid ); + m_data->CStore.Set("channelkey_to_mrdpmtid", channelkey_to_mrdpmtid ); + m_data->CStore.Set("channelkey_to_faccpmtid", channelkey_to_faccpmtid); - // tank PMTs - for(int pmti=0; pmtiGetPMT(pmti); - - // Construct the detector associated with this PMT - //unsigned long uniquedetectorkey = anniegeom->ConsumeNextFreeDetectorKey(); - unsigned long uniquedetectorkey = pmtid_to_channelkey[pmti+1]; - std::string CylLocString; - int cylloc = apmt.GetCylLoc(); - if(apmt.GetName().find("EMI9954KB")!= std::string::npos){ cylloc=6; } // FIXME set OD pmts in WCSim - switch (cylloc){ - case 0: CylLocString = "TopCap"; break; - case 2: CylLocString = "BottomCap"; break; - case 1: CylLocString = "Barrel"; break; - case 4: CylLocString = "MRD"; break; // TODO set this as H or V paddle? And layer? - case 5: CylLocString = "Veto"; break; // TODO set layer? - case 6: CylLocString = "OD"; break; - default: CylLocString = "NA"; break; // unknown - } - Detector adet(uniquedetectorkey, - "Tank", - CylLocString, - Position( apmt.GetPosition(0)/100., - apmt.GetPosition(1)/100., - apmt.GetPosition(2)/100.), - Direction(apmt.GetOrientation(0), - apmt.GetOrientation(1), - apmt.GetOrientation(2)), - apmt.GetName(), - detectorstatus::ON, - 0.); + return anniegeom; +} + +//////////////////////////////////////////////////////////////////////////////// +void LoadWCSim::ConstructDetectors(Geometry* anniegeom, int numDets, std::string system) +{ + // Check the "system" string. Only Tank, MRD, and Veto are allowed + if (system != "Tank" && system != "MRD" && system != "Veto" && system != "LAPPD") { + logmessage = "LoadWCSim::ConstructDetectors: \"system\" must be Tank, MRD, Veto, or LAPPD; but you passed: "; + logmessage += system; + Log(logmessage, v_error, verbosity); + return; + } + + // To fill up crates, cards, and channels monotonically + // they're arbitrary for simulation + int lvl1_crt = 0; + int lvl1_crd = 0; + int lvl1_chn = 0; + int lvl2_crt = 0; + int lvl2_crd = 0; + int lvl2_chn = 0; + int hv_crt = 0; + int hv_crd = 0; + int hv_chn = 0; + + // Loop and build the detectors + // ====================================================== + for (int detIdx = 0; detIdx < numDets; ++detIdx) { + WCSimRootPMT apmt; + if (system == "Tank") apmt = wcsimrootgeom->GetPMT(detIdx); + if (system == "MRD") apmt = wcsimrootgeom->GetMRDPMT(detIdx); + if (system == "Veto") apmt = wcsimrootgeom->GetFACCPMT(detIdx); + if (system == "LAPPD") apmt = wcsimrootgeom->GetLAPPD(detIdx); - // construct the channel associated with this PMT - //unsigned long uniquechannelkey = anniegeom->ConsumeNextFreeChannelKey(); - unsigned long uniquechannelkey = uniquedetectorkey; + // Construct the detector associated with this PMT + // ====================================================== + unsigned long uniquedetectorkey; + if (system == "Tank") uniquedetectorkey = pmtid_to_channelkey[detIdx + 1]; + if (system == "MRD") uniquedetectorkey = mrdid_to_channelkey[detIdx]; + if (system == "Veto") uniquedetectorkey = fmvid_to_channelkey[detIdx]; + if (system == "LAPPD") uniquedetectorkey = anniegeom->ConsumeNextFreeDetectorKey(); + + std::string CylLocString; + int cylloc = apmt.GetCylLoc(); + + // FIXME set OD pmts in WCSim + if (apmt.GetName().find("EMI9954KB") != std::string::npos) cylloc = 6; + + switch (cylloc) { + case 0: CylLocString = "TopCap"; break; + case 2: CylLocString = "BottomCap"; break; + case 1: CylLocString = "Barrel"; break; + case 4: CylLocString = "MRD"; break; // TODO set this as H or V paddle? And layer? + case 5: CylLocString = "Veto"; break; // TODO set layer? + case 6: CylLocString = "OD"; break; + default: CylLocString = "NA"; break; // unknown + } + + Detector adet(uniquedetectorkey, system, CylLocString, + Position( apmt.GetPosition(0)/100., + apmt.GetPosition(1)/100., + apmt.GetPosition(2)/100.), + Direction(apmt.GetOrientation(0), + apmt.GetOrientation(1), + apmt.GetOrientation(2)), + apmt.GetName(), detectorstatus::ON, 0.); + + // Construct the channel associated with this PMT + // Tank, MRD, and Veto have one channel per detector + // LAPPDs have two channels per strip + // ====================================================== + unsigned long uniquechannelkey; + if (system != "LAPPD") { + // Tank, MRD, or Veto + uniquechannelkey = uniquedetectorkey; + ++lvl1_chn; + ++hv_chn; + + if (system == "Tank") { pmt_tubeid_to_channelkey.emplace(apmt.GetTubeNo(), uniquechannelkey); - if (verbosity > 3) std::cout <<"LoadWCSim tool: WCSim ID: "<=ADC_CHANNELS_PER_CARD) { ADC_Chan_Num=0; ADC_Card_Num++; MT_Chan_Num++; } - if(ADC_Card_Num>=ADC_CARDS_PER_CRATE) { ADC_Card_Num=0; ADC_Crate_Num++; } - if(MT_Chan_Num>=MT_CHANNELS_PER_CARD) { MT_Chan_Num=0; MT_Card_Num++; } - if(MT_Card_Num>=MT_CARDS_PER_CRATE) { MT_Card_Num=0; MT_Crate_Num++; } - // same for HV - CAEN_HV_Chan_Num++; - if(CAEN_HV_Chan_Num>=CAEN_HV_CHANNELS_PER_CARD) { CAEN_HV_Chan_Num=0; CAEN_HV_Card_Num++; } - if(CAEN_HV_Card_Num>=CAEN_HV_CARDS_PER_CRATE) { CAEN_HV_Card_Num=0; CAEN_HV_Crate_Num++; } - - Channel pmtchannel( uniquechannelkey, - Position(0,0,0.), - 0, // stripside - 0, // stripnum - ADC_Crate_Num, - ADC_Card_Num, - ADC_Chan_Num, - MT_Crate_Num, - MT_Card_Num, - MT_Chan_Num, - CAEN_HV_Crate_Num, - CAEN_HV_Card_Num, - CAEN_HV_Chan_Num, - channelstatus::ON); - - // Add this channel to the geometry - if(verbosity>4) cout<<"Adding channel "<4) cout<<"Adding detector "<AddDetector(adet); - if(verbosity>4) cout<<"printing geometry"<4) anniegeom->PrintChannels(); - } - - // mrd PMTs - for(int mrdpmti=0; mrdpmtiGetMRDPMT(mrdpmti); - - // Construct the detector associated with this PMT - //unsigned long uniquedetectorkey = anniegeom->ConsumeNextFreeDetectorKey(); - unsigned long uniquedetectorkey = mrdid_to_channelkey[mrdpmti]; - std::string CylLocString; - switch (apmt.GetCylLoc()){ - case 0: CylLocString = "TopCap"; break; - case 2: CylLocString = "BottomCap"; break; - case 1: CylLocString = "Barrel"; break; - case 4: CylLocString = "MRD"; break; // TODO set this as H or V paddle? And layer? - case 5: CylLocString = "Veto"; break; // TODO set layer? - default: CylLocString = "NA"; break; // unknown - } - Detector adet(uniquedetectorkey, - "MRD", - CylLocString, - Position( apmt.GetPosition(0)/100., - apmt.GetPosition(1)/100., - apmt.GetPosition(2)/100.), - Direction(apmt.GetOrientation(0), - apmt.GetOrientation(1), - apmt.GetOrientation(2)), - apmt.GetName(), - detectorstatus::ON, - 0.); - - // construct the channel associated with this PMT - //unsigned long uniquechannelkey = anniegeom->ConsumeNextFreeChannelKey(); - unsigned long uniquechannelkey = uniquedetectorkey; + channelkey_to_pmtid.emplace(uniquechannelkey, apmt.GetTubeNo()); + + if (lvl1_chn >= ADC_CHANNELS_PER_CARD) { lvl1_chn = 0; ++lvl1_crd; ++lvl2_chn; } + if (lvl1_crd >= ADC_CARDS_PER_CRATE) { lvl1_crd = 0; ++lvl1_crt; } + if (lvl2_chn >= MT_CHANNELS_PER_CARD) { lvl2_chn = 0; ++lvl2_crd; } + if (lvl2_crd >= MT_CARDS_PER_CRATE) { lvl2_crd = 0; ++lvl2_crt; } + if (hv_chn >= CAEN_HV_CHANNELS_PER_CARD) { hv_chn = 0; ++hv_crd; } + if (hv_crd >= CAEN_HV_CARDS_PER_CRATE) { hv_crd = 0; ++hv_crt; } + } + + // Technically the LeCroy is for both MRD and Veto. Thus the veto HV + // channel/card/crate should pickup from where the MRD left off + // But these things elements are never called in practice and these + // numbers are largely aribitrary for simulation anyway. + + if (system == "MRD") { mrd_tubeid_to_channelkey.emplace(apmt.GetTubeNo(), uniquechannelkey); channelkey_to_mrdpmtid.emplace(uniquechannelkey, apmt.GetTubeNo()); - - // fill up TDC cards and channels monotonically, they're arbitrary for simulation - TDC_Chan_Num++; - if(TDC_Chan_Num>=TDC_CHANNELS_PER_CARD) { TDC_Chan_Num=0; TDC_Card_Num++; } - if(TDC_Card_Num>=TDC_CARDS_PER_CRATE) { TDC_Card_Num=0; TDC_Crate_Num++; } - // same for HV - LeCroy_HV_Chan_Num++; - if(LeCroy_HV_Chan_Num>=LECROY_HV_CHANNELS_PER_CARD) { LeCroy_HV_Chan_Num=0; LeCroy_HV_Card_Num++; } - if(LeCroy_HV_Card_Num>=LECROY_HV_CARDS_PER_CRATE) { LeCroy_HV_Card_Num=0; LeCroy_HV_Crate_Num++; } - - Channel pmtchannel( uniquechannelkey, - Position(0,0,0.), - 0, // stripside - 0, // stripnum - TDC_Crate_Num, - TDC_Card_Num, - TDC_Chan_Num, - -1, // TDC has no level 2 signal handling - -1, - -1, - LeCroy_HV_Crate_Num, - LeCroy_HV_Card_Num, - LeCroy_HV_Chan_Num, - channelstatus::ON); - + + if (lvl1_chn >= TDC_CHANNELS_PER_CARD) { lvl1_chn = 0; ++lvl1_crd; } + if (lvl1_crd >= TDC_CARDS_PER_CRATE) { lvl1_crd = 0; ++lvl1_crt; } + if (hv_chn >= LECROY_HV_CHANNELS_PER_CARD) { hv_chn = 0; ++hv_crd; } + if (hv_crd >= LECROY_HV_CARDS_PER_CRATE) { hv_crd = 0; ++hv_crt; } + lvl2_crt = -1; lvl2_crd = -1; lvl2_chn = -1; + } + + if (system == "Veto") { + facc_tubeid_to_channelkey.emplace(apmt.GetTubeNo(), uniquechannelkey); + channelkey_to_faccpmtid.emplace(uniquechannelkey, apmt.GetTubeNo()); + + if (lvl1_chn >= TDC_CHANNELS_PER_CARD) { lvl1_chn = 0; ++lvl1_crd; } + if (lvl1_crd >= TDC_CARDS_PER_CRATE) { lvl1_crd = 0; ++lvl1_crt; } + if (hv_chn >= LECROY_HV_CHANNELS_PER_CARD) { hv_chn = 0; ++hv_crd; } + if (hv_crd >= LECROY_HV_CARDS_PER_CRATE) { hv_crd = 0; ++hv_crt; } + lvl2_crt = -1; lvl2_crd = -1; lvl2_chn = -1; + } + + Channel channel(uniquechannelkey, Position(0, 0, 0.), + 0, 0, + lvl1_crt, lvl1_crd, lvl1_chn, + lvl2_crt, lvl2_crd, lvl2_chn, + hv_crt, hv_crd, hv_chn, + channelstatus::ON); + + // Add this channel to the geometry + logmessage = "LoadWCSim::ConstructDetectors: " + system ; + logmessage += ", WCSim ID: " + std::to_string(apmt.GetTubeNo()); + logmessage += " Adding channel: " + std::to_string(uniquechannelkey); + logmessage += " to detector: " + std::to_string(uniquedetectorkey); + Log(logmessage, v_debuggier, verbosity); + + adet.AddChannel(channel); + } else { + // LAPPDs get one channel per strip + lappd_tubeid_to_detectorkey.emplace(apmt.GetTubeNo(), uniquedetectorkey); + detectorkey_to_lappdid.emplace(uniquedetectorkey, apmt.GetTubeNo()); + + for (int stripIdx = 0; stripIdx < LappdNumStrips; ++stripIdx) { + uniquechannelkey = anniegeom->ConsumeNextFreeChannelKey(); + + // Calc strip info + // StripSide == 0 for LHS, StripSide == 1 for RHS + int stripside = ((stripIdx % 2) == 0); + int stripnum = (int)(stripIdx / 2); + double xpos = (stripside) ? -LappdStripLength : LappdStripLength; + double ypos = (stripnum * LappdStripSeparation) - ((LappdNumStrips * LappdStripSeparation) / 2.); + ++lvl1_chn; + ++hv_chn; + if (lvl1_chn >= ACDC_CHANNELS_PER_CARD) { lvl1_chn = 0; ++lvl1_crd; ++lvl2_chn; } + if (lvl1_crd >= ACDC_CARDS_PER_CRATE) { lvl1_crd = 0; ++lvl1_crt; } + if (lvl2_chn >= ACC_CHANNELS_PER_CARD) { lvl2_chn = 0; ++lvl2_crd; } + if (lvl2_crd >= ACC_CARDS_PER_CRATE) { lvl2_crd = 0; ++lvl2_crt; } + if (hv_chn >= LAPPD_HV_CHANNELS_PER_CARD) { hv_chn = 0; ++hv_crd; } + if (hv_crd >= LAPPD_HV_CARDS_PER_CRATE) { hv_crd = 0; ++hv_crt; } + + Channel channel(uniquechannelkey, Position(xpos, ypos, 0.), + stripside, stripnum, + lvl1_crt, lvl1_crd, lvl1_chn, + lvl2_crt, lvl2_crd, lvl2_chn, + hv_crt, hv_crd, hv_chn, + channelstatus::ON); + // Add this channel to the geometry - if(verbosity>4) cout<<"Adding channel "<4) cout<<"Adding detector "<AddDetector(adet); - - // Create a paddle - // FIXME remove dependency on MRDSpecs - if(mrdpmti!=(apmt.GetTubeNo()-1)){ - cerr<<"LoadWCSim: mrdpmti!=pmt.GetTubeNo-1!"<AddDetector(adet); + + if (verbosity >= v_debuggier) anniegeom->PrintChannels(); + + // Finished for the Tank and LAPPDs + // Gotta create paddles for the MRD and Veto + // ====================================================== + if (system == "MRD" || system == "Veto") { + + // Things to fill for paddles + int pad_x, pad_y, pad_z; + int orientation = 0; + Position origin; + std::pair ext_x; + std::pair ext_y; + std::pair ext_z; + + if (system == "MRD") { + if (detIdx != (apmt.GetTubeNo()-1) ) { + logmessage = "LoadWCSim::ConstructDetectors: " + system; + logmessage += "detIdx != pmt.GetTubeNo-1! "; + logmessage += "detIdx = " + std::to_string(detIdx); + logmessage += ", pmt.GetTubeNo = " + std::to_string(apmt.GetTubeNo()); + assert(false); } + // calculate MRD_x_y_z ... MRDSpecs doesn't provide a nice way to do this - int layernum=0; - while ((mrdpmti+1) > MRDSpecs::layeroffsets.at(layernum+1)){ layernum++; } - Mrd_Chankey_Layer.emplace(uniquechannelkey,layernum); - int in_layer_pmtnum = mrdpmti - MRDSpecs::layeroffsets.at(layernum); + int layernum = 0; + while ((detIdx + 1) > MRDSpecs::layeroffsets.at(layernum + 1)) ++layernum; + + Mrd_Chankey_Layer.emplace(uniquechannelkey, layernum); + int in_layer_pmtnum = detIdx - MRDSpecs::layeroffsets.at(layernum); + // paddles in each layer alternate on sides; i.e. paddles 0 and 1 are on opposite sides - int side = in_layer_pmtnum%2; - in_layer_pmtnum = std::floor(in_layer_pmtnum/2); - int orientation = MRDSpecs::paddle_orientations.at(mrdpmti); - int MRD_x = (orientation) ? in_layer_pmtnum : side; - int MRD_y = (orientation) ? side : in_layer_pmtnum; - int MRD_z = layernum+2; // first MRD z layer num is 2 (veto are 0,1) - - Paddle apaddle( uniquedetectorkey, - MRD_x, - MRD_y, - MRD_z, - orientation, - Position( MRDSpecs::paddle_originx.at(mrdpmti)/1000., - MRDSpecs::paddle_originy.at(mrdpmti)/1000., - MRDSpecs::paddle_originz.at(mrdpmti)/1000.), - std::pair{MRDSpecs::paddle_extentsx.at(mrdpmti).first/1000., - MRDSpecs::paddle_extentsx.at(mrdpmti).second/1000.}, - std::pair{MRDSpecs::paddle_extentsy.at(mrdpmti).first/1000., - MRDSpecs::paddle_extentsy.at(mrdpmti).second/1000.}, - std::pair{MRDSpecs::paddle_extentsz.at(mrdpmti).first/1000., - MRDSpecs::paddle_extentsz.at(mrdpmti).second/1000.}); - if(verbosity>4) cout<<"Setting paddle for detector "<SetDetectorPaddle(uniquedetectorkey,apaddle); - - if(verbosity>4) cout<<"printing geometry"<4) anniegeom->PrintChannels(); - } + int side = in_layer_pmtnum % 2; + in_layer_pmtnum = std::floor(in_layer_pmtnum / 2); + orientation = MRDSpecs::paddle_orientations.at(detIdx); + pad_x = (orientation) ? in_layer_pmtnum : side; + pad_y = (orientation) ? side : in_layer_pmtnum; + pad_z = layernum + 2; // first MRD z layer num is 2 (veto are 0,1) + + origin = Position(MRDSpecs::paddle_originx.at(detIdx), + MRDSpecs::paddle_originy.at(detIdx), + MRDSpecs::paddle_originz.at(detIdx)); + + origin *= 1./1000.; + + ext_x = std::make_pair(MRDSpecs::paddle_extentsx.at(detIdx).first/1000., + MRDSpecs::paddle_extentsx.at(detIdx).second/1000.); + + ext_y = std::make_pair(MRDSpecs::paddle_extentsy.at(detIdx).first/1000., + MRDSpecs::paddle_extentsy.at(detIdx).second/1000.); + + ext_z = std::make_pair(MRDSpecs::paddle_extentsz.at(detIdx).first/1000., + MRDSpecs::paddle_extentsz.at(detIdx).second/1000.); + }// endif MRD + + else { // system == "Veto" + pad_z = (uniquedetectorkey>12); // 13 paddles per layer + pad_x = pad_z; // i believe PMTs are on LHS for layer 0, RHS for layer 1 + pad_y = uniquedetectorkey - (13*pad_z); + + // numbers from geofile.txt + origin = Position(0, facc_paddle_yorigins.at(uniquedetectorkey)/100., (pad_z) ? 0.0728 : 0.0508); + + // numbers from WCSim source files / measurements + ext_x = std::make_pair(-1.6,1.6); + ext_y = std::make_pair(origin.Y() - 0.1525, origin.Y() + 0.1525); + ext_z = std::make_pair(origin.Z() - 0.01, origin.Z() + 0.01 ); + }// endif Veto + + Paddle apaddle(uniquedetectorkey, pad_x, pad_y, pad_z, + orientation, origin, ext_x, ext_y, ext_z); + + logmessage = "LoadWCSim::ConstructDetectors: " + system; + logmessage += " Setting paddle for detector: " + std::to_string(uniquedetectorkey); + Log(logmessage, v_debuggier, verbosity); + + anniegeom->SetDetectorPaddle(uniquedetectorkey,apaddle); + }// end paddle things for MRD/Veto + }// end loop over detectors +} +//////////////////////////////////////////////////////////////////////////////// +void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) +{ + // MCParticles will only be loaded during trigger 0 of a given entry + // This means that particles from subtriggers are loaded as well + // When we reach those subtriggers then we'll simply update the relevant info + if (MCTriggerNum == 0) { + // Reset the things + MCParticles->clear(); + trackid_to_mcparticleindex->clear(); + mapNeutronIsPrim->clear(); + primaryMuonIndex = -1; + + // Loop over ALL triggers so we can load all MCParticles + for (int trigIdx = 0; trigIdx < WCSimEntry->wcsimrootevent->GetNumberOfEvents(); ++trigIdx) { + + WCSimRootTrigger* aTrigTank = WCSimEntry->wcsimrootevent->GetTrigger(trigIdx); + + // Loop over tracks within this trigger and extract MCParticles + logmessage = "LoadWCSim::LoadMCParticles: Getting " + std::to_string(aTrigTank->GetNtrack()); + logmessage += " tracks from trigger # " + std::to_string(trigIdx); + Log(logmessage, v_message, verbosity); + + for (int trackIdx = 0; trackIdx < aTrigTank->GetNtrack(); trackIdx++) { + logmessage = "LoadWCSim::LoadMCParticles: Getting WCSim track # " + std::to_string(trackIdx); + Log(logmessage, v_message, verbosity); + + auto* nextTrack = (WCSimRootTrack*)aTrigTank->GetTracks()->At(trackIdx); + + tracktype startStopType = tracktype::UNDEFINED; + + // Extract the neutrino information + if (nextTrack->GetFlag() == -1) { + double startTime = AdjustTime((double)nextTrack->GetTime()); + double stopTime = AdjustTime((double)nextTrack->GetStopTime()); + Position startPos(nextTrack->GetStart(0), nextTrack->GetStart(1), nextTrack->GetStart(2)); + Position stopPos(nextTrack->GetStop(0), nextTrack->GetStop(1), nextTrack->GetStop(2)); + startPos.UnitToMeter(); + stopPos.UnitToMeter(); + double length = (stopPos-startPos).Mag(); + + MCParticle neutrino(nextTrack->GetIpnu(), nextTrack->GetE(), nextTrack->GetEndE(), + startPos, stopPos, startTime, stopTime, + Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), + length, startStopType, + nextTrack->GetId(), + nextTrack->GetParenttype(), + nextTrack->GetFlag(), + trigIdx); + + // Save the neutrino own particle in the store + m_data->Stores["ANNIEEvent"]->Set("NeutrinoParticle", neutrino); + + // Don't also create a particle for the neutrino + continue; + + }// done extracting the neutrino information + + // Now load the other particles + double startTime = AdjustTime((double)nextTrack->GetTime()); + double stopTime = AdjustTime((double)nextTrack->GetStopTime()); + Position startPos(nextTrack->GetStart(0), nextTrack->GetStart(1), nextTrack->GetStart(2)); + Position stopPos(nextTrack->GetStop(0), nextTrack->GetStop(1), nextTrack->GetStop(2)); + startPos.UnitToMeter(); + stopPos.UnitToMeter(); + double length = (stopPos-startPos).Mag(); - // veto PMTs - for(int faccpmti=0; faccpmtiGetFACCPMT(faccpmti); - - // Construct the detector associated with this PMT - //unsigned long uniquedetectorkey = anniegeom->ConsumeNextFreeDetectorKey(); - unsigned long uniquedetectorkey = fmvid_to_channelkey[faccpmti]; - std::string CylLocString; - switch (apmt.GetCylLoc()){ - case 0: CylLocString = "TopCap"; break; - case 2: CylLocString = "BottomCap"; break; - case 1: CylLocString = "Barrel"; break; - case 4: CylLocString = "MRD"; break; // TODO set this as H or V paddle? And layer? - case 5: CylLocString = "Veto"; break; // TODO set layer? - default: CylLocString = "NA"; break; // unknown - } - Detector adet( uniquedetectorkey, - "Veto", - CylLocString, - Position( apmt.GetPosition(0)/100., - apmt.GetPosition(1)/100., - apmt.GetPosition(2)/100.), - Direction(apmt.GetOrientation(0), - apmt.GetOrientation(1), - apmt.GetOrientation(2)), - apmt.GetName(), - detectorstatus::ON, - 0.); - if (verbosity > 3) std::cout <<"LoadWCSim tool: FACC tube channelkey: "<ConsumeNextFreeChannelKey(); - unsigned long uniquechannelkey = uniquedetectorkey; - facc_tubeid_to_channelkey.emplace(apmt.GetTubeNo(), uniquechannelkey); - channelkey_to_faccpmtid.emplace(uniquechannelkey, apmt.GetTubeNo()); - - // fill up TDC cards and channels monotonically, they're arbitrary for simulation - TDC_Chan_Num++; - if(TDC_Chan_Num>=TDC_CHANNELS_PER_CARD) { TDC_Chan_Num=0; TDC_Card_Num++; } - if(TDC_Card_Num>=TDC_CARDS_PER_CRATE) { TDC_Card_Num=0; TDC_Crate_Num++; } - // same for HV - LeCroy_HV_Chan_Num++; - if(LeCroy_HV_Chan_Num>=LECROY_HV_CHANNELS_PER_CARD) { LeCroy_HV_Chan_Num=0; LeCroy_HV_Card_Num++; } - if(LeCroy_HV_Card_Num>=LECROY_HV_CARDS_PER_CRATE) { LeCroy_HV_Card_Num=0; LeCroy_HV_Crate_Num++; } - - Channel pmtchannel( uniquechannelkey, - Position(0,0,0.), - 0, // stripside - 0, // stripnum - TDC_Crate_Num, - TDC_Card_Num, - TDC_Chan_Num, - -1, // TDC has no level 2 signal handling - -1, - -1, - LeCroy_HV_Crate_Num, - LeCroy_HV_Card_Num, - LeCroy_HV_Chan_Num, - channelstatus::ON); - - // Add this channel to the geometry - if(verbosity>4) cout<<"Adding channel "<4) cout<<"Adding detector "<AddDetector(adet); - - // Create a paddle - // FIXME even MRDSpecs can't save us here; instead hard-code the values for now, - // this will all be replaced eventually anyway - int MRD_z = (uniquedetectorkey>12); // 13 paddles per layer - int MRD_x = MRD_z; // i believe PMTs are on LHS for layer 0, RHS for layer 1 - int MRD_y = uniquedetectorkey - (13*MRD_z); - double paddle_zorigin = (MRD_z) ? 0.0728 : 0.0508; // numbers from geofile.txt - double paddle_yorigin = facc_paddle_yorigins.at(uniquedetectorkey)/100.; - - Paddle apaddle( uniquedetectorkey, - MRD_x, - MRD_y, - MRD_z, - 0, // orientation 0=horizontal, 1=vertical - Position(0,paddle_yorigin,paddle_zorigin), - std::pair{-1.6,1.6}, // numbers from WCSim source files / measurements - std::pair{paddle_yorigin-0.1525,paddle_yorigin+0.1525}, - std::pair{paddle_zorigin-0.01,paddle_zorigin+0.01}); - - //std::cout <<"FACC detkey: "<4) cout<<"Setting paddle for detector "<SetDetectorPaddle(uniquedetectorkey,apaddle); - - if(verbosity>4) cout<<"printing geometry"<4) anniegeom->PrintChannels(); - } + logmessage = "LoadWCSim::LoadMCParticles: Loaded particle with PDG: " + std::to_string(nextTrack->GetIpnu()); + logmessage += ", stop time: " + std::to_string(stopTime); + logmessage += ", end process: " + nextTrack->GetEndProcess(); + Log(logmessage, v_debug, verbosity); + + // Record neutron primary/secondary + if (nextTrack->GetIpnu() == 2112) + mapNeutronIsPrim->emplace(nextTrack->GetId(), (nextTrack->GetParenttype() == 0)); - // lappds - // lappds moved to end such that all the other channels are assigned first - for(int lappdi=0; lappdiGetLAPPD(lappdi); - - // Construct the detector associated with this tile - unsigned long uniquedetectorkey = anniegeom->ConsumeNextFreeDetectorKey(); - std::cout <<"LAPPD unique detectorkey: "<GetIpnu(), nextTrack->GetE(), nextTrack->GetEndE(), + startPos, stopPos, startTime, stopTime, + Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), + length, startStopType, + nextTrack->GetId(), + nextTrack->GetParenttype(), + nextTrack->GetFlag(), + trigIdx); + + // Exit point is not currently in constructor call so set it separately + // Older WCSim files do not recor this info. This breaks backward compatibility + Position exitPoint(nextTrack->GetTankExitPoint(0), + nextTrack->GetTankExitPoint(1), + nextTrack->GetTankExitPoint(2)); + exitPoint.UnitToMeter(); + thisparticle.SetTankExitPoint(exitPoint); + + // Check if this is a primary muon. Only record the first one + if (nextTrack->GetIpnu() == 13 && nextTrack->GetParenttype() == 0 && + nextTrack->GetFlag() == 0 && primaryMuonIndex < 0 ) + primaryMuonIndex = MCParticles->size(); + + // Some print outs for "interesting" particles + if (abs(nextTrack->GetIpnu()) == 13 || abs(nextTrack->GetIpnu()) == 211 || nextTrack->GetIpnu() == 111){ + logmessage = "LoadWCSim::LoadMCParticles: Found " + std::to_string(nextTrack->GetIpnu()); + logmessage += " with flag: " + std::to_string(nextTrack->GetFlag()); + logmessage += ", parent type " + std::to_string(nextTrack->GetParenttype()); + logmessage += ", Id " + std::to_string(nextTrack->GetId()); + logmessage += ", start vertex (" + std::to_string(nextTrack->GetStart(0)/100.); + logmessage += ", " + std::to_string(nextTrack->GetStart(1)/100.); + logmessage += ", " + std::to_string(nextTrack->GetStart(2)/100.); + logmessage += "), and end vertex (" + std::to_string(nextTrack->GetStop(0)/100.); + logmessage += ", " + std::to_string(nextTrack->GetStop(1)/100.); + logmessage += ", " + std::to_string(nextTrack->GetStop(2)/100.) + ")"; + Log(logmessage, v_debug, verbosity); } - Detector adet( uniquedetectorkey, - "LAPPD", - CylLocString, - Position( anlappd.GetPosition(0)/100., - anlappd.GetPosition(1)/100., - anlappd.GetPosition(2)/100.), - Direction(anlappd.GetOrientation(0), - anlappd.GetOrientation(1), - anlappd.GetOrientation(2)), - anlappd.GetName(), - detectorstatus::ON, - 0.); - - // construct all the channels associated with this LAPPD - for(int stripi=0; stripiConsumeNextFreeChannelKey(); - - int stripside = ((stripi%2)==0); // StripSide=0 for LHS (x<0), StripSide=1 for RHS (x>0) - int stripnum = (int)(stripi/2); // Strip number: add 2 channels per strip as we go - double xpos = (stripside) ? -LappdStripLength : LappdStripLength; - double ypos = (stripnum*LappdStripSeparation) - ((LappdNumStrips*LappdStripSeparation)/2.); - - // fill up ADC cards and channels monotonically, they're arbitrary for simulation - ACDC_Chan_Num++; - if(ACDC_Chan_Num>=ACDC_CHANNELS_PER_CARD) { ACDC_Chan_Num=0; ACDC_Card_Num++; ACC_Chan_Num++; } - if(ACDC_Card_Num>=ACDC_CARDS_PER_CRATE) { ACDC_Card_Num=0; ACDC_Crate_Num++; } - if(ACC_Chan_Num>=ACC_CHANNELS_PER_CARD) { ACC_Chan_Num=0; ACC_Card_Num++; } - if(ACC_Card_Num>=ACC_CARDS_PER_CRATE) { ACC_Card_Num=0; ACC_Crate_Num++; } - // same for HV - LAPPD_HV_Chan_Num++; - if(LAPPD_HV_Chan_Num>=LAPPD_HV_CHANNELS_PER_CARD) { LAPPD_HV_Chan_Num=0; LAPPD_HV_Card_Num++; } - if(LAPPD_HV_Card_Num>=LAPPD_HV_CARDS_PER_CRATE) { LAPPD_HV_Card_Num=0; LAPPD_HV_Crate_Num++; } - - Channel lappdchannel( uniquechannelkey, - Position(xpos,ypos,0.), - stripside, - stripnum, - ACDC_Crate_Num, - ACDC_Card_Num, - ACDC_Chan_Num, - ACC_Crate_Num, - ACC_Card_Num, - ACC_Chan_Num, - LAPPD_HV_Crate_Num, - LAPPD_HV_Card_Num, - LAPPD_HV_Chan_Num, - channelstatus::ON); + + trackid_to_mcparticleindex->emplace(nextTrack->GetId(),MCParticles->size()); + MCParticles->push_back(thisparticle); + }// end loop over tracks + + logmessage = "LoadWCSim::LoadMCParticles: Loaded " + std::to_string(MCParticles->size()) + " MCParticles"; + Log(logmessage, v_debug, verbosity); + } // end loop over events + }// endif MCTriggerNum == 0 + else { + // if MCTrigger > 0 we need to update all the particle times + // since particle times are relative to the trigger time + double timeDiff = EventTimeNs - firstTrig->GetHeader()->GetDate(); + for (MCParticle& aparticle : *MCParticles) { + if (splitSubtriggers) { + aparticle.SetStartTime(aparticle.GetStartTime() - timeDiff); + aparticle.SetStopTime (aparticle.GetStopTime() - timeDiff); + } + } + } // end updating particle times +} + +//////////////////////////////////////////////////////////////////////////////// +void LoadWCSim::LoadNeutronCaptures(WCSimRootTrigger* aTrig) +{ + int numCaptures = aTrig ? aTrig->GetNcaptures() : 0; + + logmessage = "LoadWCSim::LoadNeutronCaptures: looping over "; + logmessage += std::to_string(numCaptures) + "neutron captures"; + Log(logmessage, v_message, verbosity); + + for (int captIdx = 0; captIdx < numCaptures; ++captIdx) { + logmessage = "LoadWCSim::Execute: getting capture # " + std::to_string(captIdx); + Log(logmessage, v_debug, verbosity); + + auto* capt = (WCSimRootCapture*)aTrig->GetCaptures()->At(captIdx); - // Add this channel to the geometry - if(verbosity>4) cout<<"Adding channel "<4) cout<<"Adding detector "<AddDetector(adet); - if(verbosity>4) cout<<"printing geometry"<4) anniegeom->PrintChannels(); - } - - // for other WCSim tools that may need the WCSim Tube IDs - m_data->CStore.Set("lappd_tubeid_to_detectorkey",lappd_tubeid_to_detectorkey); - m_data->CStore.Set("pmt_tubeid_to_channelkey",pmt_tubeid_to_channelkey); - m_data->CStore.Set("mrd_tubeid_to_channelkey",mrd_tubeid_to_channelkey); - m_data->CStore.Set("facc_tubeid_to_channelkey",facc_tubeid_to_channelkey); - // inverse - m_data->CStore.Set("detectorkey_to_lappdid",detectorkey_to_lappdid); - m_data->CStore.Set("channelkey_to_pmtid",channelkey_to_pmtid); - m_data->CStore.Set("channelkey_to_mrdpmtid",channelkey_to_mrdpmtid); - m_data->CStore.Set("channelkey_to_faccpmtid",channelkey_to_faccpmtid); + int captParent = capt->GetCaptureParent(); + double captVtxX = capt->GetCaptureVtx(0); + double captVtxY = capt->GetCaptureVtx(1); + double captVtxZ = capt->GetCaptureVtx(2); + int captNGamma = capt->GetNGamma(); + double captTotalE = capt->GetTotalGammaE(); + double captT = capt->GetCaptureT(); + int captNucleus = capt->GetCaptureNucleus(); + + // Adjust the time if we're splitting subtriggers + if (splitSubtriggers) captT -= EventTimeNs; + + // Check if capture is from a primary + double doubleIsPrimary = -9999; + if (mapNeutronIsPrim->count(captParent) > 0) { + bool isPrimary = mapNeutronIsPrim->at(captParent); - return anniegeom; + logmessage = "LoadWCSim::Execute: NeutCap object at " + std::to_string(captT); + logmessage += " is from a"; + logmessage += isPrimary ? " primary" : " secondary" ; + logmessage += " neutron parent with ID: " + std::to_string(captParent); + logmessage += " on nucleus: " + std::to_string(captNucleus); + Log(logmessage, v_debug, verbosity); + } else { + // The parent neutron to this capture was not recorded by WCSim + logmessage = "LoadWCSim::Execute: NeutCap object at time " + std::to_string(captT); + logmessage += " is from an unrecorded neutron parent"; + logmessage += " with ID: " + std::to_string(captParent); + logmessage += " on nucleus: " + std::to_string(captNucleus); + Log(logmessage, v_message, verbosity); + } + + // Record the photons from the neutron capture + std::vector gammaEnergies; + for (int gammaIdx = 0; gammaIdx < captNGamma; ++gammaIdx) { + auto* captGamma = (WCSimRootCaptureGamma*) capt->GetGammas()->At(gammaIdx); + gammaEnergies.push_back(captGamma->GetE()); + } + + // Push back the things + if (MCNeutCap.size()==0){ + MCNeutCap.emplace("CaptParent",std::vector{double(captParent)}); + MCNeutCap.emplace("CaptVtxX",std::vector{captVtxX}); + MCNeutCap.emplace("CaptVtxY",std::vector{captVtxY}); + MCNeutCap.emplace("CaptVtxZ",std::vector{captVtxZ}); + MCNeutCap.emplace("CaptNGamma",std::vector{double(captNGamma)}); + MCNeutCap.emplace("CaptTotalE",std::vector{captTotalE}); + MCNeutCap.emplace("CaptTime",std::vector{captT}); + MCNeutCap.emplace("CaptNucleus",std::vector{double(captNucleus)}); + MCNeutCapGammas.emplace("CaptGammas",std::vector>{gammaEnergies}); + MCNeutCap.emplace("CaptPrimary",std::vector{doubleIsPrimary}); + } else { + MCNeutCap.at("CaptParent").push_back(captParent); + MCNeutCap.at("CaptVtxX").push_back(captVtxX); + MCNeutCap.at("CaptVtxY").push_back(captVtxY); + MCNeutCap.at("CaptVtxZ").push_back(captVtxZ); + MCNeutCap.at("CaptNGamma").push_back(captNGamma); + MCNeutCap.at("CaptTotalE").push_back(captTotalE); + MCNeutCap.at("CaptTime").push_back(captT); + MCNeutCap.at("CaptNucleus").push_back(captNucleus); + MCNeutCapGammas.at("CaptGammas").push_back(gammaEnergies); + MCNeutCap.at("CaptPrimary").push_back(doubleIsPrimary); + } + }// end loop over neutron captures } -void LoadWCSim::MakeParticleToPmtMap(WCSimRootTrigger* thistrig, WCSimRootTrigger* firstTrig, std::map>* ParticleId_to_TubeIds, std::map* ParticleId_to_Charge, std::map tubeid_to_channelkey){ - if(thistrig==nullptr) return; - ParticleId_to_TubeIds->clear(); - ParticleId_to_Charge->clear(); - // scan through the parents IDs of the photons contributing to each digit - // make note of which parent contributes to which digit, and which digits are associated with each parent - Log("Making Particle to PMT Map",v_message,verbosity); - // technically the charge will be a lower limit as this sums the charge from all digits - // that a given particle contributed to, but not all this digit's charge may have been - // from this particle. +//////////////////////////////////////////////////////////////////////////////// +bool LoadWCSim::LoadHits(WCSimRootTrigger* thisTrig, WCSimRootTrigger* firstTrig, std::string system) +{ + // Check the "system" string. Only Tank, MRD, and Veto are allowed + if (system != "Tank" && system != "MRD" && system != "Veto") { + logmessage = "LoadWCSim::LoadHits: \"system\" must be Tank, MRD, or Veto; but you passed: "; + logmessage += system; + Log(logmessage, v_error, verbosity); + return false; + } + + int numDigiHits = thisTrig ? thisTrig->GetCherenkovDigiHits()->GetEntries() : 0; + logmessage = "LoadWCSim::LoadHits: Looping over " + std::to_string(numDigiHits); + logmessage += " " + system + " digi. hits"; + Log(logmessage, v_message, verbosity); + + for (int hitIdx = 0; hitIdx < numDigiHits; ++hitIdx) { + logmessage = "LoadWCSim::LoadHits: Getting hit: " + std::to_string(hitIdx); + Log(logmessage, v_debug, verbosity); + + auto* digiHit = (WCSimRootCherenkovDigiHit*)thisTrig->GetCherenkovDigiHits()->At(hitIdx); + int tubeID = digiHit->GetTubeId(); + + if ((system == "Tank" && pmt_tubeid_to_channelkey.count(tubeID) == 0) || + (system == "MRD" && mrd_tubeid_to_channelkey.count(tubeID) == 0) || + (system == "Veto" && facc_tubeid_to_channelkey.count(tubeID) == 0) ) { + logmessage = "LoadWCSim::LoadHits: NO PMT ASSOCIATED WITH ID: " + std::to_string(tubeID); + Log(logmessage, v_error, verbosity); + return false; + } + + uint64_t key; + if (system == "Tank") key = pmt_tubeid_to_channelkey.at(tubeID); + if (system == "MRD") key = mrd_tubeid_to_channelkey.at(tubeID); + if (system == "Veto") key = facc_tubeid_to_channelkey.at(tubeID); + + logmessage = "LoadWCSim::LoadHits: tubeid: " + std::to_string(tubeID); + logmessage += " has ChannelKey: " + std::to_string(key); + Log(logmessage, v_debug, verbosity); + + // Omit masked tank PMTs + if (system == "Tank" && PMTMask != "None" && + std::find(masked_ids.begin(), masked_ids.end(), tubeID) != masked_ids.end()) { + logmessage = "LoadWCSim::LoadHits: Skipping masked PMT: " + std::to_string(tubeID); + Log(logmessage, v_debug, verbosity); + continue; + } + + // Grab the hit charge and time + float digiQ = digiHit->GetQ(); + double digiTime; + if (use_smeared_digit_time) { + // relative to trigger + digiTime = static_cast(digiHit->GetT()-HistoricTriggeroffset); + } else { + // Use the true time of the first photon in the hit + double earliestTime = 999999999999; + + // Loop over indices of the digit's photons not "IDs" + std::vector photonIdxs = digiHit->GetPhotonIds(); + for (int& photonIdx : photonIdxs) { + auto* theHitTimeObject = (WCSimRootCherenkovHitTime*)firstTrig->GetCherenkovHitTimes()->At(photonIdx); + + if (theHitTimeObject == nullptr) { + logmessage = "LoadWCSim::LoadHits: Retrieval of photon from digi. hit returned nullptr!"; + Log(logmessage, v_error, verbosity); + continue; + } + + double thisTime = static_cast(theHitTimeObject->GetTruetime()); + if (thisTime < earliestTime) earliestTime = thisTime; + }// end loop over photons + digiTime = earliestTime; + }// endif use_smeared_digit_time + + // Adjust hit time as necessary based on settings and system + if ( (system == "Tank" && !splitSubtriggers && use_smeared_digit_time) || + (system == "MRD" && !splitSubtriggers) || + (system == "Veto") ) { + digiTime += EventTimeNs; + } + + logmessage = "LoadWCSim::LoadHits: digi. hit time is " + std::to_string(digiTime); + logmessage +=" [ns] from the start of the Trigger"; + logmessage += " and charge is: " + std::to_string(digiQ); + Log(logmessage, v_debug, verbosity); + + // Create the hit and put it in the correct map + MCHit nextHit(key, digiTime, digiQ, GetHitParentIdxs(digiHit, firstTrig)); + + if (system == "Tank") { + if (MCHits->count(key) == 0) MCHits->emplace(key, std::vector{nextHit}); + else MCHits->at(key).push_back(nextHit); + } + + if (system == "MRD") { + if (TDCData->count(key) == 0) TDCData->emplace(key, std::vector{nextHit}); + else TDCData->at(key).push_back(nextHit); + + if (Mrd_Chankey_Layer.at(key) == 0) mrd_firstlayer = true; + if (Mrd_Chankey_Layer.at(key) == 10) mrd_lastlayer = true; + } + + if (system == "Veto") { + if (TDCData->count(key) == 0) TDCData->emplace(key, std::vector{nextHit}); + else TDCData->at(key).push_back(nextHit); + } + + logmessage = "LoadWCSim::LoadHits: Tank hit added for " + system; + Log(logmessage, v_debug, verbosity); + + }// end loop over hits + + logmessage = "LoadWCSim::LoadHits: Finished loading " + system + " hits"; + Log(logmessage, v_message, verbosity); + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +void LoadWCSim::MakeParticleToPmtMap(WCSimRootTrigger* thistrig, + WCSimRootTrigger* firstTrig, + std::map>* ParticleId_to_TubeIds, + std::map* ParticleId_to_Charge, + std::map tubeid_to_channelkey) +{ + if (thistrig==nullptr) return; + ParticleId_to_TubeIds->clear(); + ParticleId_to_Charge->clear(); + // scan through the parents IDs of the photons contributing to each digit + // make note of which parent contributes to which digit, and which digits are associated with each parent + Log("Making Particle to PMT Map",v_message, verbosity); + // technically the charge will be a lower limit as this sums the charge from all digits + // that a given particle contributed to, but not all this digit's charge may have been + // from this particle. - // To match digits to their parent particles we need the corresponding CherenkovHitTimes - // both CherenkovHits and CherenkovHitTimes are stored in the first trigger - //-------------------------------------------------------------------------------------- + // To match digits to their parent particles we need the corresponding CherenkovHitTimes + // both CherenkovHits and CherenkovHitTimes are stored in the first trigger + //-------------------------------------------------------------------------------------- - // Loop over all digits - int numdigits = thistrig->GetCherenkovDigiHits()->GetEntries(); - for(int digiti=0; digitiGetCherenkovDigiHits()->At(digiti); - int tubeid = digihit->GetTubeId(); - // loop over the photons in this digit - std::vector truephotonindices = digihit->GetPhotonIds(); - for(int truephoton=0; truephoton<(int)truephotonindices.size(); truephoton++){ - int thephotonsid = truephotonindices.at(truephoton); - // get the index of the photon CherenkovHit object in the TClonesArray - if(WCSimVersion<2){ - if(timeArrayOffsetMap.size()==0) BuildTimeArrayOffsetMap(firstTrig); - thephotonsid+=timeArrayOffsetMap.at(tubeid); } - // Get the CherenkovHitTime object that records the photon's Parent ID - WCSimRootCherenkovHitTime *thehittimeobject = - (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(thephotonsid)); - if(thehittimeobject==nullptr) cerr<<"HITTIME IS NULL"<GetParentID() : -1; - // We'll want a map of particle ID to channel keys, so convert WCSim TubeId to channelkey - int channelkey = tubeid_to_channelkey.at(tubeid); + // Loop over all digits + int numdigits = thistrig->GetCherenkovDigiHits()->GetEntries(); + for (int digitIdx=0; digitIdxGetCherenkovDigiHits()->At(digitIdx); + int tubeID = digihit->GetTubeId(); + double hitQ = digihit->GetQ(); + + + // loop over the photons in this digit + std::vector truephotonindices = digihit->GetPhotonIds(); + for (int thephotonsid : truephotonindices) { + + // get the index of the photon CherenkovHit object in the TClonesArray + if (WCSimVersion < 2) { + if (!timeArrayOffsetMap.size()) BuildTimeArrayOffsetMap(firstTrig); + thephotonsid += timeArrayOffsetMap.at(tubeID); + } + + // Get the CherenkovHitTime object that records the photon's Parent ID + auto* thehittimeobject = (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(thephotonsid)); + + // get the parent ID from the CherenkovHitTime + Int_t parentID = (thehittimeobject) ? thehittimeobject->GetParentID() : -1; + + // We'll want a map of particle ID to channel keys, so convert WCSim TubeID to channelkey + int chankey = tubeid_to_channelkey.at(tubeID); - // Finally we can record this pairing of Tube to Particle - if(ParticleId_to_TubeIds->count(thephotonsparenttrackid)==0){ - // we've not recorded any hits for this particle: make an empty map for it - ParticleId_to_TubeIds->emplace(thephotonsparenttrackid,std::map{}); - } - if(ParticleId_to_TubeIds->at(thephotonsparenttrackid).count(channelkey)==0){ - // in the map for this particle record that this tube was hit - ParticleId_to_TubeIds->at(thephotonsparenttrackid).emplace(channelkey,digihit->GetQ()); -// -// double particletime = -1; -// int particletrigger = -1; -// if(trackid_to_mcparticleindex->count(thephotonsparenttrackid)){ -// int particleindex = trackid_to_mcparticleindex->at(thephotonsparenttrackid); -// MCParticle theparticle = MCParticles->at(particleindex); -// particletime = theparticle.GetStartTime(); -// particletrigger = theparticle.GetMCTriggerNum(); -// } -// double hittime = digihit->GetT(); -// double triggertime = thistrig->GetHeader()->GetDate(); -// std::cout<<"Particle "<at(thephotonsparenttrackid).at(channelkey)+=digihit->GetQ(); - } - if(ParticleId_to_Charge->count(thephotonsparenttrackid)==0){ - // first time seeing this particle - ParticleId_to_Charge->emplace(thephotonsparenttrackid,digihit->GetQ()); - } else { - ParticleId_to_Charge->at(thephotonsparenttrackid)+=digihit->GetQ(); - } - } - } // end loop over digits + // Finally we can record the things + // Check for parentID not in outer map, place it there empty + auto outerIt = ParticleId_to_TubeIds->find(parentID); + if (outerIt == ParticleId_to_TubeIds->end()) + outerIt = ParticleId_to_TubeIds->emplace(parentID, std::map()).first; + + // Now emplace the charge value + // Will succeed (result.second == true) if tubeID is not present + // Will fail if tubeID is present. In which case we'll increment the charge + auto result_PartToTubeIDs = outerIt->second.emplace(tubeID, hitQ); + if (!result_PartToTubeIDs.second) + result_PartToTubeIDs.first->second += hitQ; + + // Similar for the total charge map, but this is not nested + auto result_PartToCharge = ParticleId_to_Charge->emplace(parentID, hitQ); + if (!result_PartToCharge.second) + result_PartToCharge.first->second += hitQ; + + }// end loop over photons + }// end loop over digits } -std::vector LoadWCSim::GetHitParentIds(WCSimRootCherenkovDigiHit* digihit, WCSimRootTrigger* firstTrig){ - /* Get the ID of the MCParticle(s) that produced this digit */ - std::vector parentids; // a hit could technically have more than one contrbuting particle +//////////////////////////////////////////////////////////////////////////////// +// Get the ID of the primary MCParticle(s) that produced this digi. hit +std::vector LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) +{ + std::vector parentIDs; // a hit could technically have more than one contrbuting particle - // loop over the photons in this digit - std::vector truephotonindices = digihit->GetPhotonIds(); - for(int truephoton=0; truephoton<(int)truephotonindices.size(); truephoton++){ - int thephotonsid = truephotonindices.at(truephoton); - // get the indices of the digit's photon CherenkovHitTime objects - if(WCSimVersion<2){ - if(timeArrayOffsetMap.size()==0) BuildTimeArrayOffsetMap(firstTrig); - thephotonsid+=timeArrayOffsetMap.at(digihit->GetTubeId()); - } - // get the CherenkovHitTime objects themselves, which contain the photon parent IDs - WCSimRootCherenkovHitTime *thehittimeobject = - (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(thephotonsid)); - if(thehittimeobject==nullptr) cerr<<"HITTIME IS NULL"<GetParentID(); - // check if this parent track was saved. Not all particles are saved. - if(trackid_to_mcparticleindex->count(theparenttrackid)){ - parentids.push_back(trackid_to_mcparticleindex->at(theparenttrackid)); - } // else this photon may have come from e.g. an electron or gamma that wasn't recorded - } - } - return parentids; + // loop over the photons in this digit + std::vector photonIdxs = digiHit->GetPhotonIds(); + for (int photonIdx : photonIdxs) { + // Special offset for older WCSim + if (WCSimVersion < 2) { + if (timeArrayOffsetMap.size() == 0) BuildTimeArrayOffsetMap(firstTrig); + photonIdx += timeArrayOffsetMap.at(digiHit->GetTubeId()); + } + + // Get the CherenkovHitTime objects themselves, which contain the primary parent IDs + auto* theHitTimeObject = (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(photonIdx)); + + if (theHitTimeObject == nullptr) { + logmessage = "LoadWCSim::GetHitParentIDs: HitTime object is NULL!!"; + Log(logmessage, v_error, verbosity); + } + else + parentIDs.push_back(theHitTimeObject->GetParentID()); + + }// end loop over photons + return parentIDs; } -void LoadWCSim::BuildTimeArrayOffsetMap(WCSimRootTrigger* firstTrig){ - if(WCSimVersion<2){ - // The CherenkovHitTimes is a flattened array (over PMTs) of arrays (over photons) - // For WCSimVersion<2, the PhotonIds available from a digit are the indices - // *within the subarray for that PMT* - // we therefore need we need to offset these indices by the start of the pmt's subarray. - // This offset may be found by scanning the CherenkovHits array (over PMTs), - // finding the correct TubeID, and retrieving the 'GetTotalPe(0)' member for this entry. - int ncherenkovhits=firstTrig->GetCherenkovHits()->GetEntries(); //atrigt->GetNcherenkovhits(); - //int nhittimes = firstTrig->GetCherenkovHitTimes()->GetEntries(); - - for(int ihit = 0; ihit < ncherenkovhits; ihit++){ - // each WCSimRootCherenkovHit represents a hit PMT - WCSimRootCherenkovHit* hitobject = - (WCSimRootCherenkovHit*)firstTrig->GetCherenkovHits()->At(ihit); - if(hitobject==nullptr) cerr<<"HITOBJECT IS NULL!"<GetTubeID(); - int timeArrayOffset = hitobject->GetTotalPe(0); - timeArrayOffsetMap.emplace(tubeNumber,timeArrayOffset); - } - } else { - Log("LoadWCSim Tool: BuildTimeArrayOffsetMap called with WCSimVersion<2: This is not needed!?",v_error,verbosity); - } +//////////////////////////////////////////////////////////////////////////////// +// Get the index within the the MCParticle vector of the primaries that produced this digi. hit +std::vector LoadWCSim::GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) +{ + std::vector parentIDs = GetHitParentIDs(digiHit, firstTrig); + std::vector parentIdxs; + + // Check if the parent was recorded, and if so then translate ID to index + for (int parentID : parentIDs) { + if (trackid_to_mcparticleindex->count(parentID)) + parentIdxs.push_back(trackid_to_mcparticleindex->at(parentID)); + } + + return parentIdxs; +} + +//////////////////////////////////////////////////////////////////////////////// +void LoadWCSim::BuildTimeArrayOffsetMap(WCSimRootTrigger* firstTrig) +{ + if (WCSimVersion < 2) { + // The CherenkovHitTimes is a flattened array (over PMTs) of arrays (over photons) + // For WCSimVersion<2, the PhotonIds available from a digit are the indices + // *within the subarray for that PMT* + // we therefore need we need to offset these indices by the start of the pmt's subarray. + // This offset may be found by scanning the CherenkovHits array (over PMTs), + // finding the correct TubeID, and retrieving the 'GetTotalPe(0)' member for this entry. + int ncherenkovhits = firstTrig->GetCherenkovHits()->GetEntries(); + + // each WCSimRootCherenkovHit represents a hit PMT + for (int hitIdx = 0; hitIdx < ncherenkovhits; ++hitIdx){ + auto* hitobject = (WCSimRootCherenkovHit*)firstTrig->GetCherenkovHits()->At(hitIdx); + if (hitobject == nullptr) { + logmessage = "LoadWCSim::BuildTimeArrayOffsetMap: Hit object is NULL!!"; + Log(logmessage, v_error, verbosity); + } + + int tubeNumber = hitobject->GetTubeID(); + int timeArrayOffset = hitobject->GetTotalPe(0); + timeArrayOffsetMap.emplace(tubeNumber, timeArrayOffset); + } + } else { + logmessage = "LoadWCSim::BuildTimeArrayOffsetMap: Called with WCSimVersion < 2. "; + logmessage += "This is not needed!?"; + Log(logmessage, v_error, verbosity); + } } -std::vector LoadWCSim::LoadPMTMask(std::string path_to_pmtmask){ +//////////////////////////////////////////////////////////////////////////////// +std::vector LoadWCSim::LoadPMTMask(std::string path_to_pmtmask) +{ - std::vector mask_vector; - int temp_id; - ifstream maskfile(path_to_pmtmask.c_str()); - while (!maskfile.eof()){ - maskfile >> temp_id; - mask_vector.push_back(temp_id); - if (maskfile.eof()) break; - } + std::vector mask_vector; + int temp_id; + ifstream maskfile(path_to_pmtmask.c_str()); + while (!maskfile.eof()){ + maskfile >> temp_id; + mask_vector.push_back(temp_id); + if (maskfile.eof()) break; + } - return mask_vector; + return mask_vector; +} + +//////////////////////////////////////////////////////////////////////////////// +double LoadWCSim::AdjustTime(double time) +{ + if (splitSubtriggers) return time - EventTimeNs; + else return time; } diff --git a/UserTools/LoadWCSim/LoadWCSim.h b/UserTools/LoadWCSim/LoadWCSim.h index c85e12454..ad2069ec5 100644 --- a/UserTools/LoadWCSim/LoadWCSim.h +++ b/UserTools/LoadWCSim/LoadWCSim.h @@ -23,25 +23,25 @@ namespace{ //PMTs - constexpr int ADC_CHANNELS_PER_CARD=4; - constexpr int ADC_CARDS_PER_CRATE=20; - constexpr int MT_CHANNELS_PER_CARD=4; - constexpr int MT_CARDS_PER_CRATE=20; + constexpr int ADC_CHANNELS_PER_CARD = 4; + constexpr int ADC_CARDS_PER_CRATE = 20 ; + constexpr int MT_CHANNELS_PER_CARD = 4; + constexpr int MT_CARDS_PER_CRATE = 20; //LAPPDs - constexpr int ACDC_CHANNELS_PER_CARD=30; - constexpr int ACDC_CARDS_PER_CRATE=20; - constexpr int ACC_CHANNELS_PER_CARD=8; - constexpr int ACC_CARDS_PER_CRATE=20; + constexpr int ACDC_CHANNELS_PER_CARD = 30; + constexpr int ACDC_CARDS_PER_CRATE = 20; + constexpr int ACC_CHANNELS_PER_CARD = 8; + constexpr int ACC_CARDS_PER_CRATE = 20; //TDCs - constexpr int TDC_CHANNELS_PER_CARD=32; - constexpr int TDC_CARDS_PER_CRATE=6; + constexpr int TDC_CHANNELS_PER_CARD = 32; + constexpr int TDC_CARDS_PER_CRATE = 6; //HV - constexpr int CAEN_HV_CHANNELS_PER_CARD=16; - constexpr int CAEN_HV_CARDS_PER_CRATE=10; - constexpr int LECROY_HV_CHANNELS_PER_CARD=16; - constexpr int LECROY_HV_CARDS_PER_CRATE=16; - constexpr int LAPPD_HV_CHANNELS_PER_CARD=4; // XXX ??? XXX - constexpr int LAPPD_HV_CARDS_PER_CRATE=10; // XXX ??? XXX + constexpr int CAEN_HV_CHANNELS_PER_CARD = 16; + constexpr int CAEN_HV_CARDS_PER_CRATE = 10; + constexpr int LECROY_HV_CHANNELS_PER_CARD = 16; + constexpr int LECROY_HV_CARDS_PER_CRATE = 16; + constexpr int LAPPD_HV_CHANNELS_PER_CARD = 4; // XXX ??? XXX + constexpr int LAPPD_HV_CARDS_PER_CRATE = 10; // XXX ??? XXX } class LoadWCSim: public Tool { @@ -53,6 +53,11 @@ class LoadWCSim: public Tool { bool Execute(); bool Finalise(); std::vector LoadPMTMask(std::string path_to_pmtmask); + void LoadMCParticles(WCSimRootTrigger* firstTrig); + void LoadNeutronCaptures(WCSimRootTrigger* aTrig); + bool LoadHits(WCSimRootTrigger* thisTrig, WCSimRootTrigger* firstTrig, std::string system); + + double AdjustTime(double time); private: @@ -76,9 +81,9 @@ class LoadWCSim: public Tool { ////////////////// //TFile* file; //TTree* wcsimtree; - wcsimT* WCSimEntry; // from makeclass - WCSimRootTrigger* atrigt, *atrigm, *atrigv; - WCSimRootTrigger* firsttrigt, *firsttrigm, *firsttrigv; + wcsimT* WCSimEntry; + WCSimRootTrigger* aTrigTank, *aTrigMRD, *aTrigVeto; + WCSimRootTrigger* firstTrigTank, *firstTrigMRD, *firstTrigVeto; WCSimRootGeom* wcsimrootgeom; WCSimRootOptions* wcsimrootopts; int WCSimVersion; // WCSim version @@ -98,30 +103,48 @@ class LoadWCSim: public Tool { // For constructing ToolChain Geometry ////////////////////////////////////// Geometry* ConstructToolChainGeometry(); - std::map lappd_tubeid_to_detectorkey; - std::map pmt_tubeid_to_channelkey; - std::map mrd_tubeid_to_channelkey; - std::map facc_tubeid_to_channelkey; + void ConstructDetectors(Geometry* anniegeom, int numDets, std::string system); + + std::map lappd_tubeid_to_detectorkey; + std::map pmt_tubeid_to_channelkey; + std::map mrd_tubeid_to_channelkey; + std::map facc_tubeid_to_channelkey; // inverse - std::map detectorkey_to_lappdid; - std::map channelkey_to_pmtid; - std::map channelkey_to_mrdpmtid; - std::map channelkey_to_faccpmtid; - // alternatively? Better? Save the parentage in each MCHit. Each MCHit will contain - // the index of it's parent MCParticle in the MCParticles vector - std::map* trackid_to_mcparticleindex=nullptr; - std::vector GetHitParentIds(WCSimRootCherenkovDigiHit* digihit, WCSimRootTrigger* firstTrig); - std::map timeArrayOffsetMap; + std::map detectorkey_to_lappdid; + std::map channelkey_to_pmtid; + std::map channelkey_to_mrdpmtid; + std::map channelkey_to_faccpmtid; + + // map to store if a neutron is a primary or secondary, + // key = neutron track ID, value == true if prim. + std::map* mapNeutronIsPrim = nullptr; + + // Each MCHit will contain the idx of it's parent MCParticle's + // position within the MCParticles vector + std::map* trackid_to_mcparticleindex = nullptr; + std::vector GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); + std::vector GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); + + std::map timeArrayOffsetMap; void BuildTimeArrayOffsetMap(WCSimRootTrigger* firstTrig); - int triggers_event; - std::map pmtid_to_channelkey; - std::map mrdid_to_channelkey; - std::map fmvid_to_channelkey; + int trigsInEntry; + + std::map pmtid_to_channelkey; + std::map mrdid_to_channelkey; + std::map fmvid_to_channelkey; // FIXME temporary!! remove me when we have a better way to get FACC paddle origins - std::vector facc_paddle_yorigins{-198.699875000, -167.999875000, -137.299875000, -106.599875000, -75.899875000, -45.199875000, -14.499875000, 16.200125000, 46.900125000, 77.600125000, 108.300125000, 139.000125000, 169.700125000, -198.064875000, -167.364875000, -136.664875000, -105.964875000, -75.264875000, -44.564875000, -13.864875000, 16.835125000, 47.535125000, 78.235125000, 108.935125000, 139.635125000, 170.335125000}; // taken from geofile.txt + // taken from geofile.txt + std::vector facc_paddle_yorigins = + {-198.699875000, -167.999875000, -137.299875000, -106.599875000, + -75.899875000, -45.199875000, -14.499875000, 16.200125000, + 46.900125000, 77.600125000, 108.300125000, 139.000125000, + 169.700125000, -198.064875000, -167.364875000, -136.664875000, + -105.964875000, -75.264875000, -44.564875000, -13.864875000, + 16.835125000, 47.535125000, 78.235125000, 108.935125000, + 139.635125000, 170.335125000}; //////////////// // things that will be filled into the store from this WCSim file. @@ -130,8 +153,8 @@ class LoadWCSim: public Tool { // from WCSim to Raw and the proper RawReader Tools // bool MCFlag=true; std::string MCFile; - uint64_t MCEventNum; - uint16_t MCTriggernum; + uint64_t WCSimEntryNum; + uint16_t MCTriggerNum; uint32_t RunNumber; uint32_t SubrunNumber; uint32_t EventNumber; // will need to be tracked separately, since we flatten triggers @@ -140,37 +163,42 @@ class LoadWCSim: public Tool { TimeClass RunStartTime; // as set from user uint64_t EventTimeNs; std::vector* MCParticles; - std::map>* TDCData; - std::map>* MCHits; + std::map>* TDCData; + std::map>* MCHits; std::vector* TriggerData; BeamStatus beamstat; std::map> MCNeutCap; std::map>> MCNeutCapGammas; - //BeamStatusClass* BeamStatus; - - int primarymuonindex; + + int primaryMuonIndex; // additional info - void MakeParticleToPmtMap(WCSimRootTrigger* thisTrig, WCSimRootTrigger* firstTrig, std::map>* ParticleId_to_DigitIds, std::map* ChargeFromParticleId, std::map tubeid_to_channelkey); - std::map>* ParticleId_to_TankTubeIds = nullptr; - std::map>* ParticleId_to_MrdTubeIds = nullptr; - std::map>* ParticleId_to_VetoTubeIds = nullptr; - std::map* ParticleId_to_TankCharge = nullptr; - std::map* ParticleId_to_MrdCharge = nullptr; - std::map* ParticleId_to_VetoCharge = nullptr; - std::map Mrd_Chankey_Layer; + void MakeParticleToPmtMap(WCSimRootTrigger* thisTrig, + WCSimRootTrigger* firstTrig, + std::map>* ParticleId_to_DigitIds, + std::map* ChargeFromParticleId, + std::map tubeid_to_channelkey); + + std::map>* ParticleId_to_TankTubeIds = nullptr; + std::map>* ParticleId_to_MrdTubeIds = nullptr; + std::map>* ParticleId_to_VetoTubeIds = nullptr; + std::map* ParticleId_to_TankCharge = nullptr; + std::map* ParticleId_to_MrdCharge = nullptr; + std::map* ParticleId_to_VetoCharge = nullptr; + std::map Mrd_Chankey_Layer; bool mrd_firstlayer, mrd_lastlayer; - std::string Triggertype; + std::string TriggerType; int TriggerWord; // verbosity levels: if 'verbosity' < this level, the message type will be logged. - int v_error=0; - int v_warning=1; - int v_message=2; - int v_debug=3; + int v_error = 0; + int v_warning = 1; + int v_message = 2; + int v_debug = 3; + int v_debuggier = 4; std::string logmessage; - int get_ok; + int get_ok; }; diff --git a/UserTools/PlotWaveforms/Makefile b/UserTools/PlotWaveforms/Makefile index 92ed92d59..4243d55da 100755 --- a/UserTools/PlotWaveforms/Makefile +++ b/UserTools/PlotWaveforms/Makefile @@ -1,3 +1,5 @@ +.NOTPARALLEL: dict.cc + ifndef CXXFLAGS CXXFLAGS = -std=c++1y -O3 endif @@ -14,7 +16,7 @@ ifneq ($(MAKECMDGOALS),clean) CXX = g++ CXXFLAGS += -Wall -Wextra -Wpedantic CXXFLAGS += -Werror -Wno-error=unused-parameter -Wcast-align - + # Add extra compiler flags for recognized compilers (currently just gcc # and clang) CXXVERSION = $(shell $(CXX) --version) @@ -22,7 +24,7 @@ ifneq ($(MAKECMDGOALS),clean) ifneq (,$(findstring clang,$(CXXVERSION))) # clang $(info Compiling using version $(COMPILER_VERSION) of clang) - + # The ROOT headers trigger clang's no-keyword-macro warning, so disable it. CXXFLAGS += -Wno-keyword-macro else @@ -43,7 +45,7 @@ ifneq ($(MAKECMDGOALS),clean) endif endif endif - + ROOTCONFIG := $(shell command -v root-config 2> /dev/null) # prefer rootcling as the dictionary generator executable name, but use # rootcint if you can't find it @@ -52,7 +54,7 @@ ifneq ($(MAKECMDGOALS),clean) ROOTCLING := $(shell command -v rootcint 2> /dev/null) endif ROOT := $(shell command -v root 2> /dev/null) - + ifndef ROOTCONFIG $(error Could not find a valid ROOT installation.) else @@ -69,12 +71,17 @@ ifneq ($(MAKECMDGOALS),clean) endif -dictionary: viewer_linkdef.hh +dict.cc: viewer_linkdef.hh RawViewer.h @echo "making $@" $(RM) dict.* dict*.pcm - $(ROOTCLING) -f dict.cc -c -I ../../DataModel $(BoostInclude) RawViewer.h viewer_linkdef.hh + MAKEFLAGS= $(ROOTCLING) -f dict.cc -c -I ../../DataModel $(BoostInclude) RawViewer.h viewer_linkdef.hh + +# dict.cc: viewer_linkdef.hh RawViewer.h +# @echo "making $@" +# $(RM) dict.* dict*.pcm +# $(ROOTCLING) -f dict.cc -c -I ../../DataModel $(BoostInclude) RawViewer.h viewer_linkdef.hh -dict.o: dictionary dict.cc +dict.o: dict.cc @echo "making $@" @cp ../../DataModel/ANNIEconstants.h ../../DataModel/Constants.h ../../DataModel/RawReader.h ../../DataModel/RawReadout.h ../../DataModel/RawChannel.h ../../DataModel/RawCard.h ../../DataModel/RawTrigData.h ../../include/ $(CXX) $(CXXFLAGS) -fPIC -o $@ -I../../include -L../lib $(ROOT_CXXFLAGS) \