Public Member Functions | Private Member Functions | Private Attributes | List of all members
util::ConvertEdep2Art Class Reference
Inheritance diagram for util::ConvertEdep2Art:
art::EDProducer art::detail::Producer art::detail::LegacyModule art::Modifier art::ModuleBase art::ProductRegistryHelper

Public Member Functions

 ConvertEdep2Art (fhicl::ParameterSet const &pset)
 Standard constructor and destructor for an FMWK module. More...
 
virtual ~ConvertEdep2Art ()
 
void produce (::art::Event &evt)
 
void beginJob ()
 
void beginRun (::art::Run &run)
 
void endRun (::art::Run &run)
 
- Public Member Functions inherited from art::EDProducer
 EDProducer (fhicl::ParameterSet const &pset)
 
template<typename Config >
 EDProducer (Table< Config > const &config)
 
std::string workerType () const
 
- Public Member Functions inherited from art::detail::Producer
virtual ~Producer () noexcept
 
 Producer (fhicl::ParameterSet const &)
 
 Producer (Producer const &)=delete
 
 Producer (Producer &&)=delete
 
Producer & operator= (Producer const &)=delete
 
Producer & operator= (Producer &&)=delete
 
void doBeginJob (SharedResources const &resources)
 
void doEndJob ()
 
void doRespondToOpenInputFile (FileBlock const &fb)
 
void doRespondToCloseInputFile (FileBlock const &fb)
 
void doRespondToOpenOutputFiles (FileBlock const &fb)
 
void doRespondToCloseOutputFiles (FileBlock const &fb)
 
bool doBeginRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doEndRun (RunPrincipal &rp, ModuleContext const &mc)
 
bool doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)
 
bool doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)
 
- Public Member Functions inherited from art::Modifier
 ~Modifier () noexcept
 
 Modifier ()
 
 Modifier (Modifier const &)=delete
 
 Modifier (Modifier &&)=delete
 
Modifier & operator= (Modifier const &)=delete
 
Modifier & operator= (Modifier &&)=delete
 
- Public Member Functions inherited from art::ModuleBase
virtual ~ModuleBase () noexcept
 
 ModuleBase ()
 
ModuleDescription const & moduleDescription () const
 
void setModuleDescription (ModuleDescription const &)
 
std::array< std::vector< ProductInfo >, NumBranchTypes > const & getConsumables () const
 
void sortConsumables (std::string const &current_process_name)
 
template<typename T , BranchType BT>
ViewToken< T > consumesView (InputTag const &tag)
 
template<typename T , BranchType BT>
ViewToken< T > mayConsumeView (InputTag const &tag)
 

Private Member Functions

unsigned int GetDetNumber (std::string volname) const
 
unsigned int GetStaveNumber (std::string volname) const
 
unsigned int GetModuleNumber (std::string volname) const
 
unsigned int GetLayerNumber (std::string volname) const
 
unsigned int GetSliceNumber (std::string volname) const
 
double VisibleEnergyDeposition (const TG4HitSegment *hit, bool applyBirks) const
 
bool CheckProcess (std::string process_name) const
 
unsigned int GetParentage (unsigned int trkid) const
 
void AddHits (std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_Deposits, std::vector< gar::sdp::CaloDeposit > &fDeposits)
 
std::map< int, size_t > TrackIDToMCTruthIndexMap () const
 
bool isMCPMatch (const simb::MCParticle &p1, const simb::MCParticle &p2) const
 
size_t FindMCTruthIndex (std::vector< simb::MCTruth > *mctruthcol, const simb::MCParticle &part) const
 

Private Attributes

std::string fEDepSimfile
 
std::string fGhepfile
 
bool fOverlay
 
std::string fECALMaterial
 
std::string fTPCMaterial
 
double fEnergyCut
 
bool fApplyBirks
 
double fTotPOT
 Total POT. More...
 
double fTotGoodPOT
 Total good POT. More...
 
unsigned int fSpills
 Total of spills. More...
 
unsigned int fGoodSpills
 Total of good spills. More...
 
TChain * fTreeChain
 
unsigned int nEntries
 
TChain * fGTreeChain
 
unsigned int nEntriesGhep
 
genie::NtpMCEventRecord * fMCRec
 
TG4Event * fEvent
 
double fBirksCoeff
 
bool fkeepEMShowers
 
std::string fEMShowerDaughterMatRegex
 keep EM shower daughters only in these materials More...
 
genie::PDGLibrary * pdglib
 
const gar::geo::GeometryCore * fGeo
 geometry information More...
 
const gar::detinfo::ECALProperties * fEcalProp
 
sim::ParticleList * fParticleList
 
bool fHasGHEP
 
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_ECALDeposits
 
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_TrackerDeposits
 
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_MuIDDeposits
 
std::vector< gar::sdp::CaloDeposit > fECALDeposits
 
std::vector< gar::sdp::EnergyDeposit > fGArDeposits
 
std::vector< gar::sdp::CaloDeposit > fTrackerDeposits
 
std::vector< gar::sdp::CaloDeposit > fMuIDDeposits
 
int fSpillCount
 
std::vector< int > fStartSpill
 
std::vector< int > fStopSpill
 
std::map< unsigned int, unsigned int > fTrkIDParent
 
std::map< int, size_t > fTrackIDToMCTruthIndex
 
gar::geo::BitFieldCoder * fFieldDecoderTrk
 
const gar::geo::seg::MinervaSegmentationAlg * fMinervaSegAlg
 

Additional Inherited Members

- Public Types inherited from art::EDProducer
using ModuleType = EDProducer
 
using WorkerType = WorkerT< EDProducer >
 
- Public Types inherited from art::detail::Producer
template<typename UserConfig , typename KeysToIgnore = void>
using Table = Modifier::Table< UserConfig, KeysToIgnore >
 
- Public Types inherited from art::Modifier
template<typename UserConfig , typename UserKeysToIgnore = void>
using Table = ProducerTable< UserConfig, detail::ModuleConfig, UserKeysToIgnore >
 
- Static Public Member Functions inherited from art::EDProducer
static void commitEvent (EventPrincipal &ep, Event &e)
 
- Protected Member Functions inherited from art::ModuleBase
ConsumesCollector & consumesCollector ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > consumes (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > consumesView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void consumesMany ()
 
template<typename T , BranchType = InEvent>
ProductToken< T > mayConsume (InputTag const &)
 
template<typename Element , BranchType = InEvent>
ViewToken< Element > mayConsumeView (InputTag const &)
 
template<typename T , BranchType = InEvent>
void mayConsumeMany ()
 

Detailed Description

Definition at line 88 of file ConvertEdep2Art_module.cc.

Constructor & Destructor Documentation

util::ConvertEdep2Art::ConvertEdep2Art ( fhicl::ParameterSet const &  pset)
explicit

Standard constructor and destructor for an FMWK module.

Definition at line 181 of file ConvertEdep2Art_module.cc.

182  : art::EDProducer{pset},
183  fEDepSimfile( pset.get< std::string >("EDepSimFile", "") ),
184  fGhepfile( pset.get< std::string >("GhepFile", "") ),
185  fOverlay( pset.get< bool >("OverlayFile", true) ),
186  fECALMaterial( pset.get< std::string >("ECALMaterial", "Scintillator") ),
187  fTPCMaterial( pset.get< std::string >("TPCMaterial", "HP_ArCH4") ),
188  fEnergyCut( pset.get< double >("EnergyCut", 0.000001) ),
189  fApplyBirks( pset.get< bool >("applyBirks", true) ),
190  fTreeChain(new TChain("EDepSimEvents")),
191  fGTreeChain(new TChain("gtree")),
192  fMCRec(nullptr),
193  fEvent(nullptr),
194  fkeepEMShowers( pset.get< bool >("keepEMShowers", true) ),
195  fEMShowerDaughterMatRegex( pset.get< std::string >("EMShowerDaughterMatRegex", ".*") ),
196  fParticleList(new sim::ParticleList()),
197  fHasGHEP(false)
198  {
199  pdglib = genie::PDGLibrary::Instance();
200 
201  if(fEDepSimfile.empty())
202  {
203  throw cet::exception("ConvertEdep2Art")
204  << "Empty edep-sim file";
205  }
206 
207  //If ghep file is provided
208  if(not fGhepfile.empty()) {
209  fGTreeChain->Add(fGhepfile.c_str());
210  nEntriesGhep = fGTreeChain->GetEntries();
211  fGTreeChain->SetBranchAddress("gmcrec", &fMCRec);
212  fHasGHEP = true;
213  }
214 
215  fTreeChain->Add(fEDepSimfile.c_str());
216  nEntries = fTreeChain->GetEntries();
217  fTreeChain->SetBranchAddress("Event", &fEvent);
218 
219  fGeo = gar::providerFrom<gar::geo::GeometryGAr>();
220  fEcalProp = gar::providerFrom<gar::detinfo::ECALPropertiesService>();
221  fBirksCoeff = fEcalProp->ScintBirksConstant(); //CLHEP::mm/CLHEP::MeV;
222 
223  fSpillCount = 0;
224 
225  if(!fkeepEMShowers){
226  MF_LOG_DEBUG("ConvertEdep2Art")
227  << " Will not keep EM shower daughters!";
228  }
229 
230  std::string fEncoding = fGeo->GetMinervaCellIDEncoding();
231  fFieldDecoderTrk = new gar::geo::BitFieldCoder( fEncoding );
232  fMinervaSegAlg = (gar::geo::seg::MinervaSegmentationAlg*)fGeo->MinervaSegmentationAlg();
233 
234  produces< gar::sumdata::RunData, ::art::InRun >();
235  produces< gar::sumdata::POTSummary, ::art::InRun >();
236 
237  produces< std::vector<simb::MCTruth> >();
238  if(fHasGHEP) {
239  produces< std::vector<gar::sdp::GenieParticle> >();
240  produces< std::vector<simb::GTruth> >();
241  produces< art::Assns<simb::MCTruth, simb::GTruth> >();
242  }
243  produces< art::Assns<simb::MCTruth, simb::MCParticle> >();
244  produces< std::vector<simb::MCParticle> >();
245 
246  produces< std::vector<gar::sdp::EnergyDeposit> >();
247  produces< std::vector<gar::sdp::CaloDeposit> >("ECAL");
248  produces< std::vector<gar::sdp::CaloDeposit> >("TrackerSc");
249  produces< std::vector<gar::sdp::CaloDeposit> >("MuID");
250  // produces< std::vector<gar::sdp::LArDeposit> >();
251 
252  produces< art::Assns<gar::sdp::EnergyDeposit, simb::MCParticle> >();
253  produces< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> >("ECAL");
254  produces< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> >("TrackerSc");
255  produces< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> >("MuID");
256  // produces< art::Assns<gar::sdp::LArDeposit, simb::MCParticle> >();
257  }
string
std::string string
Definition: nybbler.cc:12
util::ConvertEdep2Art::fEcalProp
const gar::detinfo::ECALProperties * fEcalProp
Definition: ConvertEdep2Art_module.cc:149
util::ConvertEdep2Art::fParticleList
sim::ParticleList * fParticleList
Definition: ConvertEdep2Art_module.cc:151
util::ConvertEdep2Art::fMinervaSegAlg
const gar::geo::seg::MinervaSegmentationAlg * fMinervaSegAlg
Definition: ConvertEdep2Art_module.cc:172
util::ConvertEdep2Art::fFieldDecoderTrk
gar::geo::BitFieldCoder * fFieldDecoderTrk
Definition: ConvertEdep2Art_module.cc:171
util::ConvertEdep2Art::fGeo
const gar::geo::GeometryCore * fGeo
geometry information
Definition: ConvertEdep2Art_module.cc:148
util::ConvertEdep2Art::fMCRec
genie::NtpMCEventRecord * fMCRec
Definition: ConvertEdep2Art_module.cc:139
gar::geo::BitFieldCoder
Helper class for decoding and encoding a bit field of 64bits for convenient declaration.
Definition: BitFieldCoder.h:130
util::ConvertEdep2Art::fEMShowerDaughterMatRegex
std::string fEMShowerDaughterMatRegex
keep EM shower daughters only in these materials
Definition: ConvertEdep2Art_module.cc:145
gar::geo::GeometryCore::MinervaSegmentationAlg
gar::geo::seg::SegmentationAlg const * MinervaSegmentationAlg() const
Returns the object handling the Sc Tracker segmentation.
Definition: GeometryCore.h:936
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
gar::geo::GeometryCore::GetMinervaCellIDEncoding
const std::string GetMinervaCellIDEncoding() const
Definition: GeometryCore.cxx:847
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
util::ConvertEdep2Art::~ConvertEdep2Art ( )
virtual

Definition at line 261 of file ConvertEdep2Art_module.cc.

261  {
262  if ( fTreeChain ) delete fTreeChain;
263  if ( fGTreeChain ) delete fGTreeChain;
264  if ( fFieldDecoderTrk ) delete fFieldDecoderTrk;
265  }
util::ConvertEdep2Art::fFieldDecoderTrk
gar::geo::BitFieldCoder * fFieldDecoderTrk
Definition: ConvertEdep2Art_module.cc:171

Member Function Documentation

void util::ConvertEdep2Art::AddHits ( std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > >  m_Deposits,
std::vector< gar::sdp::CaloDeposit > &  fDeposits 
)
private

Definition at line 432 of file ConvertEdep2Art_module.cc.

433  {
434  //Loop over the hits in the map and add them together
435 
436  for(auto &it : m_Deposits) {
437 
438  gar::raw::CellID_t cellID = it.first;
439  std::vector<gar::sdp::CaloDeposit> vechit = it.second;
440  std::sort(vechit.begin(), vechit.end()); //sort per time
441 
442  float esum = 0.;
443  float stepLength = 0.;
444  float time = vechit.at(0).Time();
445  int trackID = vechit.at(0).TrackID();
446  double pos[3] = { vechit.at(0).X(), vechit.at(0).Y(), vechit.at(0).Z() };
447 
448  for(auto const &hit : vechit) {
449  esum += hit.Energy();
450  stepLength += hit.StepLength();
451  }
452 
453  fDeposits.emplace_back( trackID, time, esum, pos, cellID, stepLength );
454  //remove the element from the map now
455  // m_Deposits.erase(it.first);
456  }
457  }
gar::raw::CellID_t
long long int CellID_t
Definition: CaloRawDigit.h:24
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
void util::ConvertEdep2Art::beginJob ( )
virtual

Reimplemented from art::EDProducer.

Definition at line 268 of file ConvertEdep2Art_module.cc.

268  {
269  fTotPOT = 0.;
270  fTotGoodPOT = 0.;
271  fSpills = 0;
272  fGoodSpills = 0;
273 
274  if(fOverlay && fHasGHEP){
275  //need to get where is the end of the spill
276  for(int ientry = 0; ientry < fGTreeChain->GetEntries(); ientry++)
277  {
278  fGTreeChain->GetEntry(ientry);
279 
280  genie::NtpMCRecHeader rec_header = fMCRec->hdr;
281  genie::EventRecord *event = fMCRec->event;
282 
283  MF_LOG_DEBUG("ConvertEdep2Art")
284  << rec_header
285  << *event;
286 
287  //Need to get the Rootino showing the end of spill
288  genie::GHepParticle *neutrino = event->Probe();
289  if(nullptr == neutrino){
290  if(fSpillCount == 0)
291  fStartSpill.insert(fStartSpill.begin()+fSpillCount, 0);
292  else{
293  fStartSpill.insert(fStartSpill.begin()+fSpillCount, fStopSpill.at(fSpillCount-1)+1);
294  }
295 
296  fStopSpill.insert(fStopSpill.begin()+fSpillCount, ientry);
297  fSpillCount++;
298  }
299  }
300 
301  MF_LOG_INFO("ConvertEdep2Art::beginJob()")
302  << "Number of spills in the ghep file "
303  << fSpillCount;
304 
305  fSpills = fSpillCount;
306  fGoodSpills = fSpillCount;
307  }
308 
309  return;
310  }
genie::NtpMCRecordI::hdr
NtpMCRecHeader hdr
record header
Definition: NtpMCRecordI.h:38
util::ConvertEdep2Art::fSpills
unsigned int fSpills
Total of spills.
Definition: ConvertEdep2Art_module.cc:131
util::ConvertEdep2Art::fMCRec
genie::NtpMCEventRecord * fMCRec
Definition: ConvertEdep2Art_module.cc:139
util::ConvertEdep2Art::fGoodSpills
unsigned int fGoodSpills
Total of good spills.
Definition: ConvertEdep2Art_module.cc:132
MF_LOG_INFO
#define MF_LOG_INFO(category)
Definition: MessageLogger.h:225
genie::EventRecord
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
util::ConvertEdep2Art::fStartSpill
std::vector< int > fStartSpill
Definition: ConvertEdep2Art_module.cc:164
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
genie::NtpMCRecHeader
MINOS-style Ntuple Class to hold an MC Event Record Header.
Definition: NtpMCRecHeader.h:32
util::ConvertEdep2Art::fTotGoodPOT
double fTotGoodPOT
Total good POT.
Definition: ConvertEdep2Art_module.cc:130
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
genie::NtpMCEventRecord::event
EventRecord * event
event
Definition: NtpMCEventRecord.h:55
void util::ConvertEdep2Art::beginRun ( ::art::Run &  run)

Definition at line 313 of file ConvertEdep2Art_module.cc.

313  {
314  //Get RunData
315  std::unique_ptr<gar::sumdata::RunData> runcol(new gar::sumdata::RunData(fGeo->DetectorName()));
316  run.put(std::move(runcol));
317 
318  //Get POT
319  if(fHasGHEP)
320  {
321  fTotPOT = fGTreeChain->GetWeight();
322  fTotGoodPOT = fTotPOT;
323  }
324 
325  MF_LOG_INFO("ConvertEdep2Art") << "POT for this file is " << fTotPOT
326  << " The number of spills is " << fSpills;
327 
328  std::unique_ptr<gar::sumdata::POTSummary> p(new gar::sumdata::POTSummary());
329  p->totpot = fTotPOT;
330  p->totgoodpot = fTotGoodPOT;
331  p->totspills = fSpills;
332  p->goodspills = fGoodSpills;
333  run.put(std::move(p));
334 
335  return;
336  }
util::ConvertEdep2Art::fGeo
const gar::geo::GeometryCore * fGeo
geometry information
Definition: ConvertEdep2Art_module.cc:148
util::ConvertEdep2Art::fSpills
unsigned int fSpills
Total of spills.
Definition: ConvertEdep2Art_module.cc:131
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
util::ConvertEdep2Art::fGoodSpills
unsigned int fGoodSpills
Total of good spills.
Definition: ConvertEdep2Art_module.cc:132
test.p
p
Definition: test.py:223
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
MF_LOG_INFO
#define MF_LOG_INFO(category)
Definition: MessageLogger.h:225
gar::geo::GeometryCore::DetectorName
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
Definition: GeometryCore.h:495
util::ConvertEdep2Art::fTotGoodPOT
double fTotGoodPOT
Total good POT.
Definition: ConvertEdep2Art_module.cc:130
bool util::ConvertEdep2Art::CheckProcess ( std::string  process_name) const
private

Definition at line 406 of file ConvertEdep2Art_module.cc.

406  {
407  bool isEMShowerProcess = false;
408 
409  if( process_name.find("EM") != std::string::npos )
410  isEMShowerProcess = true;
411 
412  return isEMShowerProcess;
413  }
void util::ConvertEdep2Art::endRun ( ::art::Run &  run)

Definition at line 339 of file ConvertEdep2Art_module.cc.

339  {
340  return;
341  }
size_t util::ConvertEdep2Art::FindMCTruthIndex ( std::vector< simb::MCTruth > *  mctruthcol,
const simb::MCParticle &  part 
) const
private

Definition at line 470 of file ConvertEdep2Art_module.cc.

470  {
471 
472  size_t mctruthindex = 0;
473  bool found = false;
474  for(size_t index = 0; index < mctruthcol->size(); index++) {
475  simb::MCTruth mctrh = mctruthcol->at(index);
476  for( size_t ipart = 0; ipart < (size_t)mctrh.NParticles(); ipart++ ) {
477  if( isMCPMatch(mctrh.GetParticle(ipart), part) ) {
478 
479  MF_LOG_DEBUG("ConvertEdep2Art") << "FindMCTruthIndex() \n"
480  << mctrh.GetParticle(ipart) << "\n"
481  << part << "\n";
482 
483  mctruthindex = index;
484  found = true;
485  break;
486  }
487  }
488  if( found ) break;
489  }
490 
491  return mctruthindex;
492  }
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
simb::MCTruth::GetParticle
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:76
util::ConvertEdep2Art::isMCPMatch
bool isMCPMatch(const simb::MCParticle &p1, const simb::MCParticle &p2) const
Definition: ConvertEdep2Art_module.cc:460
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
unsigned int util::ConvertEdep2Art::GetDetNumber ( std::string  volname) const
private

Definition at line 354 of file ConvertEdep2Art_module.cc.

354  {
355  unsigned int det_id = 0;
356  if( volname.find("Barrel") != std::string::npos )
357  det_id = 1;
358  if( volname.find("Endcap") != std::string::npos )
359  det_id = 2;
360  return det_id;
361  }
unsigned int util::ConvertEdep2Art::GetLayerNumber ( std::string  volname) const
private

Definition at line 344 of file ConvertEdep2Art_module.cc.

344  {
345  return std::atoi( (volname.substr( volname.find("layer_") + 6, 2)).c_str() );
346  }
unsigned int util::ConvertEdep2Art::GetModuleNumber ( std::string  volname) const
private

Definition at line 369 of file ConvertEdep2Art_module.cc.

369  {
370  return std::atoi( (volname.substr( volname.find("_module") + 7, 2)).c_str() );
371  }
unsigned int util::ConvertEdep2Art::GetParentage ( unsigned int  trkid) const
private

Definition at line 416 of file ConvertEdep2Art_module.cc.

416  {
417  unsigned int parentid = -1;
418  std::map<unsigned int, unsigned int>::const_iterator it = fTrkIDParent.find(trkid);
419  while ( it != fTrkIDParent.end() )
420  {
421  MF_LOG_DEBUG("ConvertEdep2Art::GetParentage")
422  << "parentage for track id " << trkid
423  << " is " << (*it).second;
424 
425  parentid = (*it).second;
426  it = fTrkIDParent.find(parentid);
427  }
428  return parentid;
429  }
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
util::ConvertEdep2Art::fTrkIDParent
std::map< unsigned int, unsigned int > fTrkIDParent
Definition: ConvertEdep2Art_module.cc:167
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
unsigned int util::ConvertEdep2Art::GetSliceNumber ( std::string  volname) const
private

Definition at line 349 of file ConvertEdep2Art_module.cc.

349  {
350  return std::atoi( (volname.substr( volname.find("slice") + 5, 1)).c_str() );
351  }
unsigned int util::ConvertEdep2Art::GetStaveNumber ( std::string  volname) const
private

Definition at line 364 of file ConvertEdep2Art_module.cc.

364  {
365  return std::atoi( (volname.substr( volname.find("_stave") + 6, 2)).c_str() );
366  }
bool util::ConvertEdep2Art::isMCPMatch ( const simb::MCParticle &  p1,
const simb::MCParticle &  p2 
) const
private

Definition at line 460 of file ConvertEdep2Art_module.cc.

461  {
462  int ulp = 2;
463  if( std::fabs(p1.E() - p2.E()) <= std::numeric_limits<float>::epsilon() * std::fabs(p1.E()+p2.E()) * ulp && p1.PdgCode() == p2.PdgCode() )
464  return true;
465 
466  return false;
467  }
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
void util::ConvertEdep2Art::produce ( ::art::Event &  evt)

The main routine of this module: Fetch the primary particles from the event, simulate their evolution in the detctor, and produce the detector response.

Definition at line 495 of file ConvertEdep2Art_module.cc.

495  {
496  MF_LOG_DEBUG("ConvertEdep2Art") << "produce()";
497  art::EventNumber_t eventnumber = evt.id().event();
498 
499  if( eventnumber > nEntries ){
500  //stop...
501  return;
502  }
503 
504  std::unique_ptr< std::vector<simb::MCTruth> > mctruthcol( new std::vector<simb::MCTruth> );
505  std::unique_ptr< std::vector<simb::GTruth> > gtruthcol( new std::vector<simb::GTruth> );
506  std::unique_ptr< art::Assns<simb::MCTruth, simb::GTruth> > tgassn( new art::Assns<simb::MCTruth, simb::GTruth> );
507 
508  std::unique_ptr< std::vector<gar::sdp::GenieParticle> > geniepartcol( new std::vector<gar::sdp::GenieParticle> );
509 
510  art::PtrMaker<simb::MCTruth> makeMCTruthPtr(evt);
511  std::vector< ::art::Ptr<simb::MCTruth> > mctPtrs;
512 
513  //--------------------------------------------------------------------------
514  //Get the event
515  //Starts at 0, evt starts at 1
516  fTreeChain->GetEntry(eventnumber-1);
517 
518  unsigned int_idx = 0;
519 
520  //--------------------------------------------------------------------------
521  if(fHasGHEP)
522  {
523  if(fOverlay) {
524  for(int ientry = fStartSpill[eventnumber-1]; ientry < fStopSpill[eventnumber-1]; ientry++)
525  {
526  fGTreeChain->GetEntry(ientry);
527 
528  genie::NtpMCRecHeader rec_header = fMCRec->hdr;
529  genie::EventRecord *event = fMCRec->event;
530  // genie::Interaction *interaction = event->Summary();
531 
532  MF_LOG_DEBUG("ConvertEdep2Art") << rec_header;
533  MF_LOG_DEBUG("ConvertEdep2Art") << *event;
534  // MF_LOG_INFO("ConvertEdep2Art") << *interaction;
535 
536  genie::GHepParticle *neutrino = event->Probe();
537  //avoid rootino events
538  if(nullptr == neutrino) {
539  int_idx = 0;
540  continue;
541  }
542 
543  //Store the list of GENIE particles
544  genie::GHepParticle* p = nullptr;
545  TObjArrayIter piter(event);
546  unsigned int idx = 0;
547 
548  while( (p = (genie::GHepParticle*) piter.Next()) ) {
549  geniepartcol->emplace_back(int_idx, idx, p->Pdg(), p->Status(), p->Name(), p->FirstMother(), p->LastMother(), p->FirstDaughter(), p->LastDaughter(), p->Px(), p->Py(), p->Pz(), p->E(), p->Mass(), p->RescatterCode());
550  idx++;
551  }
552 
553  simb::MCTruth mctruth;
554  simb::GTruth gtruth;
555 
556  evgb::FillMCTruth(event, 0., mctruth);
557  evgb::FillGTruth(event, gtruth);
558 
559  mctruthcol->push_back(mctruth);
560  gtruthcol->push_back(gtruth);
561 
562  //Make a vector of mctruth art ptr
563  art::Ptr<simb::MCTruth> MCTruthPtr = makeMCTruthPtr(mctruthcol->size() - 1);
564  mctPtrs.push_back(MCTruthPtr);
565 
566  evgb::util::CreateAssn(*this, evt, *mctruthcol, *gtruthcol, *tgassn, gtruthcol->size()-1, gtruthcol->size());
567  int_idx++;
568  }
569  }
570  else {
571  //Starts at 0, evt starts at 1
572  fGTreeChain->GetEntry(eventnumber-1);
573 
574  genie::NtpMCRecHeader rec_header = fMCRec->hdr;
575  genie::EventRecord *event = fMCRec->event;
576  // genie::Interaction *interaction = event->Summary();
577 
578  MF_LOG_DEBUG("ConvertEdep2Art") << rec_header;
579  MF_LOG_DEBUG("ConvertEdep2Art") << *event;
580  // MF_LOG_INFO("ConvertEdep2Art") << *interaction;
581 
582  //Store the list of GENIE particles
583  genie::GHepParticle* p = nullptr;
584  TObjArrayIter piter(event);
585  unsigned int idx = 0;
586 
587  while( (p = (genie::GHepParticle*) piter.Next()) ) {
588  geniepartcol->emplace_back(int_idx, idx, p->Pdg(), p->Status(), p->Name(), p->FirstMother(), p->LastMother(), p->FirstDaughter(), p->LastDaughter(), p->Px(), p->Py(), p->Pz(), p->E(), p->Mass(), p->RescatterCode());
589  idx++;
590  }
591 
592  simb::MCTruth mctruth;
593  simb::GTruth gtruth;
594 
595  evgb::FillMCTruth(event, 0., mctruth);
596  evgb::FillGTruth(event, gtruth);
597 
598  mctruthcol->push_back(mctruth);
599  gtruthcol->push_back(gtruth);
600 
601  //Make a vector of mctruth art ptr
602  art::Ptr<simb::MCTruth> MCTruthPtr = makeMCTruthPtr(mctruthcol->size() - 1);
603  mctPtrs.push_back(MCTruthPtr);
604 
605  evgb::util::CreateAssn(*this, evt, *mctruthcol, *gtruthcol, *tgassn, gtruthcol->size()-1, gtruthcol->size());
606  }
607  }
608  else {
609 
610  //Case where things are made from particle gun! no ghep file provided. Need to create MCTruth object / GTruth object
611  simb::MCTruth truth;
612  truth.SetOrigin(simb::kSingleParticle);
613 
614  for (std::vector<TG4PrimaryVertex>::const_iterator t = fEvent->Primaries.begin(); t != fEvent->Primaries.end(); ++t)
615  {
616  TLorentzVector pos(t->Position.X() / CLHEP::cm, t->Position.Y() / CLHEP::cm, t->Position.Z() / CLHEP::cm, t->Position.T());
617  double evWeight = t->Weight;
618 
619  for (std::vector<TG4PrimaryParticle>::const_iterator p = t->Particles.begin(); p != t->Particles.end(); ++p) {
620  int trackid = p->GetTrackId();
621  std::string primary("primary");
622 
623  TLorentzVector pvec(p->Momentum.Px() * CLHEP::MeV / CLHEP::GeV, p->Momentum.Py() * CLHEP::MeV / CLHEP::GeV, p->Momentum.Pz() * CLHEP::MeV / CLHEP::GeV, p->Momentum.E() * CLHEP::MeV / CLHEP::GeV);
624 
625  simb::MCParticle part(trackid, p->GetPDGCode(), primary);
626  part.AddTrajectoryPoint(pos, pvec);
627  part.SetWeight(evWeight);
628 
629  MF_LOG_DEBUG("ConvertEdep2Art") << "Adding primary particle with "
630  << " momentum " << part.P()
631  << " position " << part.Vx() << " " << part.Vy() << " " << part.Vz();
632 
633  truth.Add(part);
634  }
635  }
636 
637  MF_LOG_DEBUG("ConvertEdep2Art") << "Adding mctruth with "
638  << " nParticles " << truth.NParticles()
639  << " Origin " << truth.Origin();
640 
641  mctruthcol->push_back(truth);
642 
643  //Make a vector of mctruth art ptr
644  art::Ptr<simb::MCTruth> MCTruthPtr = makeMCTruthPtr(mctruthcol->size() - 1);
645  mctPtrs.push_back(MCTruthPtr);
646  }
647 
648  //-----------------------------------
649 
650  std::unique_ptr< std::vector<simb::MCParticle> > partCol( new std::vector<simb::MCParticle> );
651  std::unique_ptr< art::Assns<simb::MCTruth, simb::MCParticle> > tpassn( new art::Assns<simb::MCTruth, simb::MCParticle> );
652 
653  fParticleList->clear();
654  fTrkIDParent.clear();
655  fTrackIDToMCTruthIndex.clear();
656 
657  for (std::vector<TG4Trajectory>::const_iterator t = fEvent->Trajectories.begin(); t != fEvent->Trajectories.end(); ++t)
658  {
659  int trackID = t->GetTrackId();
660  int parentID = t->GetParentId();
661  int pdg = t->GetPDGCode();
662  std::string name = t->GetName();
663 
664  //Avoid breaking.... some pdg don't exist in the library...
665  TParticlePDG *part = pdglib->Find(pdg);
666  double mass = 0.;
667  if(nullptr != part) {
668  mass = part->Mass();//in GeV
669  }
670  else{
671  MF_LOG_INFO("ConvertEdep2Art")
672  << " Could not find TParticlePDG for pdg "
673  << pdg
674  << " Mass is ut to 0 GeV";
675  }
676 
677  int process = 0;
678  int subprocess = 0;
679  std::string process_name = "unknown";
680 
681  if(parentID == -1) {
682  process_name = "primary";
683  // parentID = 0;
684  }
685  else {
686  //Get the first point that created this particle
687  // process_name = t->Points.at(0).GetProcessName();
688  process = t->Points.at(0).GetProcess();
689  subprocess = t->Points.at(0).GetSubprocess();
690  process_name = gar::util::FindProcessName( process, subprocess );
691  }// end if not a primary particle
692 
693  simb::MCParticle *fParticle = new simb::MCParticle(trackID, pdg, process_name, parentID, mass);
694 
695  for (std::vector<TG4TrajectoryPoint>::const_iterator p = t->Points.begin(); p != t->Points.end(); ++p)
696  {
697  TLorentzVector position = p->GetPosition();
698 
699  TLorentzVector fourPos(position.X() / CLHEP::cm, position.Y() / CLHEP::cm, position.Z() / CLHEP::cm, position.T() );
700  TVector3 momentum = p->GetMomentum();
701  double px = momentum.x() * CLHEP::MeV / CLHEP::GeV;
702  double py = momentum.y() * CLHEP::MeV / CLHEP::GeV;
703  double pz = momentum.z() * CLHEP::MeV / CLHEP::GeV;
704  TLorentzVector fourMom(px, py, pz, std::sqrt( px*px + py*py + pz*pz + mass*mass ));
705  // std::string process = p->GetProcessName();
706  int pt_process = p->GetProcess();
707  int pt_subprocess = p->GetSubprocess();
708  std::string pt_process_name = gar::util::FindProcessName( pt_process, pt_subprocess );
709 
710  if(p == t->Points.begin()) pt_process_name = "Start";
711  fParticle->AddTrajectoryPoint(fourPos, fourMom, pt_process_name);
712  }
713 
714  // std::string end_process = t->Points.at(t->Points.size()-1).GetProcessName();
715  process = t->Points.at(t->Points.size()-1).GetProcess();
716  subprocess = t->Points.at(t->Points.size()-1).GetSubprocess();
717  std::string end_process = gar::util::FindProcessName( process, subprocess );
718  fParticle->SetEndProcess(end_process);
719 
720  fParticleList->Add( fParticle );
721  }
722 
723  size_t mcTruthIndex = 0;
724  int fCurrentTrackID = 0;
725  //Work on MCParticle list
726  for (std::map<int, simb::MCParticle*>::iterator itPart = fParticleList->begin(); itPart != fParticleList->end(); ++itPart)
727  {
728  simb::MCParticle &part = *(itPart->second);
729  int parentID = part.Mother();
730  int trackID = part.TrackId();
731  fCurrentTrackID = trackID;
732  std::string process_name = part.Process();
733 
734  if( process_name == "primary" )
735  {
736  if(not fOverlay) {
737  mcTruthIndex = 0;
738  } else {
739  //Need to get the index correctly... Check the list of mctruth and g4 particles, match them according to type and energy and pos?
740  mcTruthIndex = FindMCTruthIndex(mctruthcol.get(), part);
741  }
742  }
743  else {
744 
745  TLorentzVector part_start = part.Trajectory().Position(0);
746  TGeoNode *node = fGeo->FindNode(part_start.X(), part_start.Y(), part_start.Z());
747  std::string material_name = "";
748 
749  if(node)
750  material_name = node->GetMedium()->GetMaterial()->GetName();
751 
752  //Skip the creation of the mcp if it is part of shower (based on process name) and only if it is not matching the material
753  //TODO debug as it does not work for anatree at the moment, problem related to pdg mom exception
754  //Maybe need to keep relation trackid, parentid
755  if(not fkeepEMShowers){
756  bool isEMShowerProcess = CheckProcess( process_name );
757  std::regex const re_material(fEMShowerDaughterMatRegex);
758  if( isEMShowerProcess && not std::regex_match(material_name, re_material) ) {
759 
760  //link trkid and parent id to be able to go back in history only for these and stop at the parent that created the shower
761  fTrkIDParent[trackID] = parentID;
762  fCurrentTrackID = -1 * GetParentage(trackID);
763 
764  MF_LOG_DEBUG("ConvertEdep2Art")
765  << " Skipping EM shower daughter "
766  << " with trackID " << trackID
767  << " with parent id " << parentID
768  << " Ultimate parentage " << GetParentage(trackID)
769  << " created with process [ " << process_name << " ]";
770 
771  fParticleList->Archive(itPart->second);
772  continue;
773  }
774  }//end not keep EM Shower particles
775 
776  if( not fParticleList->KnownParticle(parentID) ) {
777  // do add the particle to the parent id map
778  // just in case it makes a daughter that we have to track as well
779  fTrkIDParent[trackID] = parentID;
780  int pid = GetParentage(parentID);
781 
782  // if we still can't find the parent in the particle navigator,
783  // we have to give up
784  if( not fParticleList->KnownParticle(pid) ) {
785  MF_LOG_DEBUG("ConvertEdep2Art")
786  << "can't find parent id: "
787  << parentID << " in the particle list, or fTrkIDParent."
788  << " Make " << parentID << " the mother ID for track ID "
789  << fCurrentTrackID << " in the hope that it will aid debugging.";
790  }
791  else
792  parentID = pid;
793  }
794 
795  // Attempt to find the MCTruth index corresponding to the
796  // current particle. If the fCurrentTrackID is not in the
797  // map try the parent ID, if that is not there, throw an
798  // exception
799  try {
800  if(fTrackIDToMCTruthIndex.count(fCurrentTrackID) > 0 )
801  mcTruthIndex = fTrackIDToMCTruthIndex.at(fCurrentTrackID);
802  else if(fTrackIDToMCTruthIndex.count(parentID) > 0 )
803  mcTruthIndex = fTrackIDToMCTruthIndex.at(parentID);
804  }
805  catch (std::exception& e) {
806  MF_LOG_DEBUG("ConvertEdep2Art")
807  << "Cannot find MCTruth index for track id "
808  << fCurrentTrackID << " or " << parentID
809  << " exception " << e.what();
810  throw;
811  }
812  } //end not primary particle
813 
814  fTrackIDToMCTruthIndex[fCurrentTrackID] = mcTruthIndex;
815  }
816 
817  // Make link between MCTruth and MCParticles
818  size_t nGeneratedParticles = 0;
819  const std::map<int, size_t> fTrackIDToMCTruthIndex_local = this->TrackIDToMCTruthIndexMap(); //Need to make trackID to MCTruth index map
820 
821  for (std::map<int, simb::MCParticle*>::iterator itPart = fParticleList->begin(); itPart != fParticleList->end(); ++itPart)
822  {
823  simb::MCParticle& p = *(itPart->second);
824 
825  MF_LOG_DEBUG("ConvertEdep2Art")
826  << "adding mc particle with track id: "
827  << p.TrackId();
828 
829  int trackID = p.TrackId();
830 
831  MF_LOG_DEBUG("ConvertEdep2Art")
832  << " Particle with pdg " << p.PdgCode()
833  << " trackID " << p.TrackId()
834  << " parent id " << p.Mother()
835  << " created with process [ " << p.Process() << " ]"
836  << " is EM " << CheckProcess( p.Process() )
837  << " with energy " << p.E();
838 
839  partCol->push_back(std::move(p));
840 
841  try {
842  if( fTrackIDToMCTruthIndex_local.count(trackID) > 0) {
843  size_t mctidx = fTrackIDToMCTruthIndex_local.find(trackID)->second;
844  evgb::util::CreateAssn(*this, evt, *partCol, mctPtrs.at(mctidx), *tpassn, nGeneratedParticles);
845  }
846  }
847  catch ( std::exception& e ) {
848  MF_LOG_DEBUG("ConvertEdep2Art")
849  << "Cannot find MCTruth for Track Id: " << trackID
850  << " to create association between Particle and MCTruth"
851  << " exception " << e.what();
852  throw;
853  }
854 
855  fParticleList->Archive(itPart->second);
856  ++nGeneratedParticles;
857  }
858 
859  MF_LOG_DEBUG("ConvertEdep2Art") << "Finished linking MCTruth and MCParticles";
860 
861  //--------------------------------------------------------------------------
862  m_ECALDeposits.clear();
863  m_TrackerDeposits.clear();
864  m_MuIDDeposits.clear();
865  fGArDeposits.clear();
866  fECALDeposits.clear();
867  fTrackerDeposits.clear();
868  fMuIDDeposits.clear();
869 
870  //Fill simulated hits
871  for (auto d = fEvent->SegmentDetectors.begin(); d != fEvent->SegmentDetectors.end(); ++d)
872  {
873  if( d->first == "TPC_Drift1" || d->first == "TPC_Drift2" )
874  {
875  //GAr deposits
876  for (std::vector<TG4HitSegment>::const_iterator h = d->second.begin(); h != d->second.end(); ++h)
877  {
878  const TG4HitSegment *hit = &(*h);
879 
880  int trackID = hit->GetPrimaryId();
881  double edep = hit->GetEnergyDeposit() * CLHEP::MeV / CLHEP::GeV;
882  double time = (hit->GetStart().T() + hit->GetStop().T())/2 / CLHEP::ns;
883  double x = (hit->GetStart().X() + hit->GetStop().X())/2 /CLHEP::cm;
884  double y = (hit->GetStart().Y() + hit->GetStop().Y())/2 /CLHEP::cm;
885  double z = (hit->GetStart().Z() + hit->GetStop().Z())/2 /CLHEP::cm;
886  double stepLength = hit->GetTrackLength() / CLHEP::cm;
887 
888  if(edep < fEnergyCut)
889  continue;
890 
891  TGeoNode *node = fGeo->FindNode(x, y, z);//Node in cm...
892  std::string VolumeName = node->GetVolume()->GetName();
893  std::string volmaterial = node->GetMedium()->GetMaterial()->GetName();
894  if ( ! std::regex_match(volmaterial, std::regex(fTPCMaterial)) ) continue;
895 
896  fGArDeposits.emplace_back(trackID, time, edep, x, y, z, stepLength, (trackID > 0));
897  }
898  }
899  else if( d->first == "BarrelECal_vol" || d->first == "EndcapECal_vol"){
900  //ECAL deposits
901  for (std::vector<TG4HitSegment>::const_iterator h = d->second.begin(); h != d->second.end(); ++h)
902  {
903  const TG4HitSegment *hit = &(*h);
904 
905  int trackID = hit->GetPrimaryId();
906  double edep = VisibleEnergyDeposition(hit, fApplyBirks) * CLHEP::MeV / CLHEP::GeV;
907  double stepLength = hit->GetTrackLength() / CLHEP::cm;
908  double time = (hit->GetStart().T() + hit->GetStop().T())/2 / CLHEP::s;
909  double x = (hit->GetStart().X() + hit->GetStop().X())/2 /CLHEP::cm;
910  double y = (hit->GetStart().Y() + hit->GetStop().Y())/2 /CLHEP::cm;
911  double z = (hit->GetStart().Z() + hit->GetStop().Z())/2 /CLHEP::cm;
912 
913  if(edep < fEnergyCut)
914  continue;
915 
916  //Check if it is in the active material of the ECAL
917  TGeoNode *node = fGeo->FindNode(x, y, z);//Node in cm...
918  std::string VolumeName = node->GetVolume()->GetName();
919  std::string volmaterial = node->GetMedium()->GetMaterial()->GetName();
920  if ( ! std::regex_match(volmaterial, std::regex(fECALMaterial)) ) continue;
921 
922  unsigned int layer = GetLayerNumber(VolumeName); //get layer number
923  unsigned int slice = GetSliceNumber(VolumeName); // get slice number
924  unsigned int det_id = GetDetNumber(VolumeName); // 1 == Barrel, 2 = Endcap
925  unsigned int stave = GetStaveNumber(VolumeName); //get the stave number
926  unsigned int module = GetModuleNumber(VolumeName); //get the module number
927 
928  std::array<double, 3> GlobalPosCM = {x, y, z};
929  std::array<double, 3> LocalPosCM;
930  gar::geo::LocalTransformation<TGeoHMatrix> trans;
931  fGeo->WorldToLocal(GlobalPosCM, LocalPosCM, trans);
932 
933  MF_LOG_DEBUG("ConvertEdep2Art")
934  << "Sensitive volume " << d->first
935  << " Hit " << hit
936  << " in volume " << VolumeName
937  << " in material " << volmaterial
938  << " det_id " << det_id
939  << " module " << module
940  << " stave " << stave
941  << " layer " << layer
942  << " slice " << slice;
943 
944  gar::raw::CellID_t cellID = fGeo->GetCellID(node, det_id, stave, module, layer, slice, LocalPosCM);//encoding the cellID on 64 bits
945 
946  double G4Pos[3] = {0., 0., 0.}; // in cm
947  G4Pos[0] = GlobalPosCM[0];
948  G4Pos[1] = GlobalPosCM[1];
949  G4Pos[2] = GlobalPosCM[2];
950 
951  gar::sdp::CaloDeposit calohit( trackID, time, edep, G4Pos, cellID, stepLength);
952  if(m_ECALDeposits.find(cellID) != m_ECALDeposits.end())
953  m_ECALDeposits[cellID].push_back(calohit);
954  else {
955  std::vector<gar::sdp::CaloDeposit> vechit;
956  vechit.push_back(calohit);
957  m_ECALDeposits.emplace(cellID, vechit);
958  }
959  }
960  }
961  else if( d->first == "Tracker_vol" ) {
962  //Minerva Style Sc Tracker for temporary det -> triangle of base 4 cm and height 2 cm
963  for (std::vector<TG4HitSegment>::const_iterator h = d->second.begin(); h != d->second.end(); ++h)
964  {
965  const TG4HitSegment *hit = &(*h);
966 
967  int trackID = hit->GetPrimaryId();
968  double edep = VisibleEnergyDeposition(hit, fApplyBirks) * CLHEP::MeV / CLHEP::GeV;
969  double stepLength = hit->GetTrackLength() /CLHEP::cm;
970  double time = (hit->GetStart().T() + hit->GetStop().T())/2 / CLHEP::s;
971  double x = (hit->GetStart().X() + hit->GetStop().X())/2 /CLHEP::cm;
972  double y = (hit->GetStart().Y() + hit->GetStop().Y())/2 /CLHEP::cm;
973  double z = (hit->GetStart().Z() + hit->GetStop().Z())/2 /CLHEP::cm;
974 
975  if(edep < fEnergyCut)
976  continue;
977 
978  //Check if it is in the active material of the ECAL
979  TGeoNode *node = fGeo->FindNode(x, y, z);//Node in cm...
980  std::string VolumeName = node->GetVolume()->GetName();
981  std::string volmaterial = node->GetMedium()->GetMaterial()->GetName();
982  if ( ! std::regex_match(volmaterial, std::regex(fECALMaterial)) ) continue;
983 
984  unsigned int layer = GetLayerNumber(VolumeName); //get layer number
985  unsigned int slice = GetSliceNumber(VolumeName); // get slice number
986  unsigned int det_id = 3;
987 
988  std::array<double, 3> GlobalPosCM = {x, y, z};
989  std::array<double, 3> LocalPosCM;
990  gar::geo::LocalTransformation<TGeoHMatrix> trans;
991  fGeo->WorldToLocal(GlobalPosCM, LocalPosCM, trans);
992 
993  MF_LOG_DEBUG("ConvertEdep2Art")
994  << "Sensitive volume " << d->first
995  << " Hit " << hit
996  << " in volume " << VolumeName
997  << " in material " << volmaterial
998  << " det_id " << det_id
999  << " layer " << layer
1000  << " slice " << slice;
1001 
1002  gar::raw::CellID_t cellID = fGeo->GetCellID(node, det_id, 0, 0, layer, slice, LocalPosCM);//encoding the cellID on 64 bits
1003 
1004  MF_LOG_DEBUG("ConvertEdep2Art")
1005  << "Sensitive volume " << d->first
1006  << " TrackLength " << stepLength
1007  << " Energy " << edep
1008  << " cellID " << cellID
1009  << " local ( " << LocalPosCM[0] << " , " << LocalPosCM[1] << " , " << LocalPosCM[2] << " )"
1010  << " global ( " << GlobalPosCM[0] << " , " << GlobalPosCM[1] << " , " << GlobalPosCM[2] << " )";
1011 
1012  double G4Pos[3] = {0., 0., 0.}; // in cm
1013  G4Pos[0] = GlobalPosCM[0];
1014  G4Pos[1] = GlobalPosCM[1];
1015  G4Pos[2] = GlobalPosCM[2];
1016 
1017  gar::sdp::CaloDeposit calohit( trackID, time, edep, G4Pos, cellID, stepLength );
1018  if(m_TrackerDeposits.find(cellID) != m_TrackerDeposits.end())
1019  m_TrackerDeposits[cellID].push_back(calohit);
1020  else {
1021  std::vector<gar::sdp::CaloDeposit> vechit;
1022  vechit.push_back(calohit);
1023  m_TrackerDeposits.emplace(cellID, vechit);
1024  }
1025  }
1026  }
1027  else if( d->first == "MuID_vol" ) {
1028  //MuonID detector in the SPY
1029  for (std::vector<TG4HitSegment>::const_iterator h = d->second.begin(); h != d->second.end(); ++h)
1030  {
1031  const TG4HitSegment *hit = &(*h);
1032 
1033  int trackID = hit->GetPrimaryId();
1034  double stepLength = hit->GetTrackLength() /CLHEP::cm;
1035  double edep = VisibleEnergyDeposition(hit, fApplyBirks) * CLHEP::MeV / CLHEP::GeV;
1036  double time = (hit->GetStart().T() + hit->GetStop().T())/2 / CLHEP::s;
1037  double x = (hit->GetStart().X() + hit->GetStop().X())/2 /CLHEP::cm;
1038  double y = (hit->GetStart().Y() + hit->GetStop().Y())/2 /CLHEP::cm;
1039  double z = (hit->GetStart().Z() + hit->GetStop().Z())/2 /CLHEP::cm;
1040 
1041  if(edep < fEnergyCut)
1042  continue;
1043 
1044  //Check if it is in the active material of the ECAL
1045  TGeoNode *node = fGeo->FindNode(x, y, z);//Node in cm...
1046  std::string VolumeName = node->GetVolume()->GetName();
1047  std::string volmaterial = node->GetMedium()->GetMaterial()->GetName();
1048  if ( ! std::regex_match(volmaterial, std::regex(fECALMaterial)) ) continue;
1049 
1050  unsigned int layer = GetLayerNumber(VolumeName); //get layer number
1051  unsigned int slice = GetSliceNumber(VolumeName); // get slice number
1052  unsigned int det_id = 4;
1053  unsigned int stave = GetStaveNumber(VolumeName);
1054  unsigned int module = GetModuleNumber(VolumeName);
1055 
1056  std::array<double, 3> GlobalPosCM = {x, y, z};
1057  std::array<double, 3> LocalPosCM;
1058  gar::geo::LocalTransformation<TGeoHMatrix> trans;
1059  fGeo->WorldToLocal(GlobalPosCM, LocalPosCM, trans);
1060 
1061  MF_LOG_DEBUG("ConvertEdep2Art")
1062  << "Sensitive volume " << d->first
1063  << " Hit " << hit
1064  << " in volume " << VolumeName
1065  << " in material " << volmaterial
1066  << " det_id " << det_id
1067  << " layer " << layer
1068  << " slice " << slice
1069  << " stave " << stave
1070  << " module " << module;
1071 
1072  gar::raw::CellID_t cellID = fGeo->GetCellID(node, det_id, stave, module, layer, slice, LocalPosCM);//encoding the cellID on 64 bits
1073 
1074  double G4Pos[3] = {0., 0., 0.}; // in cm
1075  G4Pos[0] = GlobalPosCM[0];
1076  G4Pos[1] = GlobalPosCM[1];
1077  G4Pos[2] = GlobalPosCM[2];
1078 
1079  gar::sdp::CaloDeposit calohit( trackID, time, edep, G4Pos, cellID, stepLength );
1080  if(m_MuIDDeposits.find(cellID) != m_MuIDDeposits.end())
1081  m_MuIDDeposits[cellID].push_back(calohit);
1082  else {
1083  std::vector<gar::sdp::CaloDeposit> vechit;
1084  vechit.push_back(calohit);
1085  m_MuIDDeposits.emplace(cellID, vechit);
1086  }
1087  }
1088  }
1089  else{
1090  MF_LOG_DEBUG("ConvertEdep2Art")
1091  << "Ignoring hits for sensitive material: "
1092  << d->first;
1093  continue;
1094  }
1095  }
1096 
1097  MF_LOG_DEBUG("ConvertEdep2Art") << "Finished collection sensitive hits";
1098 
1099  //--------------------------------------------------------------------------
1100 
1101  std::unique_ptr< std::vector< gar::sdp::EnergyDeposit> > TPCCol(new std::vector<gar::sdp::EnergyDeposit> );
1102  std::unique_ptr< std::vector< gar::sdp::CaloDeposit > > ECALCol(new std::vector<gar::sdp::CaloDeposit> );
1103  std::unique_ptr< std::vector< gar::sdp::CaloDeposit > > TrackerCol(new std::vector<gar::sdp::CaloDeposit> );
1104  std::unique_ptr< std::vector< gar::sdp::CaloDeposit > > MuIDCol(new std::vector<gar::sdp::CaloDeposit> );
1105  std::unique_ptr< std::vector< gar::sdp::LArDeposit > > LArCol(new std::vector<gar::sdp::LArDeposit> );
1106 
1107  bool hasGAr = false;
1108  bool hasECAL = false;
1109  bool hasTrackerSc = false;
1110  bool hasMuID = false;
1111  bool hasLAr = false;
1112  if(fGArDeposits.size() > 0) hasGAr = true;
1113  if(m_ECALDeposits.size() > 0) hasECAL = true;
1114  if(m_TrackerDeposits.size() > 0) hasTrackerSc = true;
1115  if(m_MuIDDeposits.size() > 0) hasMuID = true;
1116 
1117  if(hasGAr) {
1118  std::sort(fGArDeposits.begin(), fGArDeposits.end());
1119 
1120  for(auto const& garhit : fGArDeposits)
1121  {
1122  MF_LOG_DEBUG("ConvertEdep2Art")
1123  << "adding GAr deposits for track id: "
1124  << garhit.TrackID();
1125  TPCCol->emplace_back(garhit);
1126  }
1127  }
1128 
1129  if(hasECAL) {
1130  this->AddHits(m_ECALDeposits, fECALDeposits);
1131  std::sort(fECALDeposits.begin(), fECALDeposits.end());
1132 
1133  for(auto const& ecalhit : fECALDeposits)
1134  {
1135  MF_LOG_DEBUG("ConvertEdep2Art")
1136  << "adding calo deposits for track id: "
1137  << ecalhit.TrackID();
1138 
1139  ECALCol->emplace_back(ecalhit);
1140  }
1141  }
1142 
1143  if(hasTrackerSc) {
1144  fMinervaSegAlg->AddHitsMinerva(m_TrackerDeposits, fTrackerDeposits);
1145  std::sort(fTrackerDeposits.begin(), fTrackerDeposits.end());
1146 
1147  for(auto const& trkhit : fTrackerDeposits)
1148  {
1149  MF_LOG_DEBUG("ConvertEdep2Art")
1150  << "adding tracker Sc deposits for track id: "
1151  << trkhit.TrackID();
1152 
1153  TrackerCol->emplace_back(trkhit);
1154  }
1155  }
1156 
1157  if(hasMuID) {
1158  this->AddHits(m_MuIDDeposits, fMuIDDeposits);
1159  std::sort(fMuIDDeposits.begin(), fMuIDDeposits.end());
1160 
1161  for(auto const& muidhit : fMuIDDeposits)
1162  {
1163  MF_LOG_DEBUG("ConvertEdep2Art")
1164  << "adding muID deposits for track id: "
1165  << muidhit.TrackID();
1166 
1167  MuIDCol->emplace_back(muidhit);
1168  }
1169  }
1170 
1171  std::unique_ptr< art::Assns<gar::sdp::EnergyDeposit, simb::MCParticle> > ghmcassn(new art::Assns<gar::sdp::EnergyDeposit, simb::MCParticle>);
1172  std::unique_ptr< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> > ehmcassn(new art::Assns<gar::sdp::CaloDeposit, simb::MCParticle>); //ECAL
1173  std::unique_ptr< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> > thmcassn(new art::Assns<gar::sdp::CaloDeposit, simb::MCParticle>); //TrackerSc
1174  std::unique_ptr< art::Assns<gar::sdp::CaloDeposit, simb::MCParticle> > mhmcassn(new art::Assns<gar::sdp::CaloDeposit, simb::MCParticle>); //MuID
1175  // std::unique_ptr< art::Assns<gar::sdp::LArDeposit, simb::MCParticle> > lhmcassn(new art::Assns<gar::sdp::LArDeposit, simb::MCParticle>); //LAr
1176 
1177  //Create assn between hits and mcp
1178  art::PtrMaker<simb::MCParticle> makeMCPPtr(evt);
1179  art::PtrMaker<gar::sdp::EnergyDeposit> makeEnergyDepositPtr(evt);
1180  art::PtrMaker<gar::sdp::CaloDeposit> makeCaloDepositPtr(evt, "ECAL");
1181  art::PtrMaker<gar::sdp::CaloDeposit> makeTrackerDepositPtr(evt, "TrackerSc");
1182  art::PtrMaker<gar::sdp::CaloDeposit> makeMuIDDepositPtr(evt, "MuID");
1183 
1184  unsigned int imcp = 0;
1185  for(auto const &part : *partCol)
1186  {
1187  int mpc_trkid = part.TrackId();
1188  art::Ptr<simb::MCParticle> partPtr = makeMCPPtr(imcp);
1189 
1190  unsigned int igashit = 0;
1191  unsigned int iecalhit = 0;
1192  unsigned int itrkhit = 0;
1193  unsigned int imuidhit = 0;
1194 
1195  if(hasGAr) {
1196  for(auto const& gashit : *TPCCol)
1197  {
1198  if(mpc_trkid == gashit.TrackID()){
1199  art::Ptr<gar::sdp::EnergyDeposit> gashitPtr = makeEnergyDepositPtr(igashit);
1200  ghmcassn->addSingle(gashitPtr, partPtr);
1201  }
1202  igashit++;
1203  }
1204  }
1205 
1206  if(hasECAL) {
1207  for(auto const& ecalhit : *ECALCol)
1208  {
1209  if(mpc_trkid == ecalhit.TrackID()){
1210  art::Ptr<gar::sdp::CaloDeposit> ecalhitPtr = makeCaloDepositPtr(iecalhit);
1211  ehmcassn->addSingle(ecalhitPtr, partPtr);
1212  }
1213  iecalhit++;
1214  }
1215  }
1216 
1217  if(hasTrackerSc) {
1218  for(auto const& trkhit : *TrackerCol)
1219  {
1220  if(mpc_trkid == trkhit.TrackID()){
1221  art::Ptr<gar::sdp::CaloDeposit> trkhitPtr = makeTrackerDepositPtr(itrkhit);
1222  thmcassn->addSingle(trkhitPtr, partPtr);
1223  }
1224  itrkhit++;
1225  }
1226  }
1227 
1228  if(hasMuID) {
1229  for(auto const& muidhit : *MuIDCol)
1230  {
1231  if(mpc_trkid == muidhit.TrackID()){
1232  art::Ptr<gar::sdp::CaloDeposit> muIDhitPtr = makeMuIDDepositPtr(imuidhit);
1233  mhmcassn->addSingle(muIDhitPtr, partPtr);
1234  }
1235  imuidhit++;
1236  }
1237  }
1238 
1239  imcp++;
1240  }
1241 
1242  evt.put(std::move(mctruthcol));
1243  if(fHasGHEP) {
1244  evt.put(std::move(gtruthcol));
1245  evt.put(std::move(tgassn));
1246  evt.put(std::move(geniepartcol));
1247  }
1248  evt.put(std::move(tpassn));
1249  evt.put(std::move(partCol));
1250  evt.put(std::move(TPCCol));
1251  evt.put(std::move(ghmcassn));
1252  evt.put(std::move(ECALCol), "ECAL");
1253  evt.put(std::move(ehmcassn), "ECAL");
1254  evt.put(std::move(TrackerCol), "TrackerSc");
1255  evt.put(std::move(thmcassn), "TrackerSc");
1256  evt.put(std::move(MuIDCol), "MuID");
1257  evt.put(std::move(mhmcassn), "MuID");
1258  if(hasLAr) {
1259  // evt.put(std::move(LArCol));
1260  // evt.put(std::move(lhmcassn));
1261  }
1262 
1263  return;
1264  }
genie::units::cm
static constexpr double cm
Definition: Units.h:68
name
static QCString name
Definition: declinfo.cpp:673
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:233
TG4Event::Primaries
TG4PrimaryVertexContainer Primaries
Definition: TG4Event.h:28
gar::geo::GeometryCore::FindNode
TGeoNode * FindNode(T const &x, T const &y, T const &z) const
util::ConvertEdep2Art::fTrackIDToMCTruthIndex
std::map< int, size_t > fTrackIDToMCTruthIndex
Definition: ConvertEdep2Art_module.cc:168
iterator
intermediate_table::iterator iterator
Definition: intermediate_table.cc:27
util::ConvertEdep2Art::TrackIDToMCTruthIndexMap
std::map< int, size_t > TrackIDToMCTruthIndexMap() const
Definition: ConvertEdep2Art_module.cc:116
util::ConvertEdep2Art::GetStaveNumber
unsigned int GetStaveNumber(std::string volname) const
Definition: ConvertEdep2Art_module.cc:364
genie::GHepParticle::RescatterCode
int RescatterCode(void) const
Definition: GHepParticle.h:65
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:212
reco_momentum_tuples.py
py
Definition: reco_momentum_tuples.py:219
gar::geo::GeometryCore::WorldToLocal
bool WorldToLocal(std::array< double, 3 > const &world, std::array< double, 3 > &local, gar::geo::LocalTransformation< TGeoHMatrix > &trans) const
Definition: GeometryCore.cxx:394
util::ConvertEdep2Art::m_TrackerDeposits
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_TrackerDeposits
Definition: ConvertEdep2Art_module.cc:155
genie::NtpMCRecordI::hdr
NtpMCRecHeader hdr
record header
Definition: NtpMCRecordI.h:38
simb::MCParticle::AddTrajectoryPoint
void AddTrajectoryPoint(TLorentzVector const &position, TLorentzVector const &momentum)
Definition: MCParticle.h:257
genie::GHepParticle::E
double E(void) const
Get energy.
Definition: GHepParticle.h:91
simb::MCTruth::SetOrigin
void SetOrigin(simb::Origin_t origin)
Definition: MCTruth.h:82
string
std::string string
Definition: nybbler.cc:12
simb::MCParticle::Trajectory
const simb::MCTrajectory & Trajectory() const
Definition: MCParticle.h:253
genie::GHepParticle::FirstDaughter
int FirstDaughter(void) const
Definition: GHepParticle.h:68
simb::MCParticle::Mother
int Mother() const
Definition: MCParticle.h:213
TG4Event::Trajectories
TG4TrajectoryContainer Trajectories
Definition: TG4Event.h:35
TG4HitSegment::GetEnergyDeposit
double GetEnergyDeposit() const
The total energy deposit in this hit.
Definition: TG4HitSegment.h:43
simb::MCTruth::Origin
simb::Origin_t Origin() const
Definition: MCTruth.h:74
bump_copyright.d
string d
Definition: bump_copyright.py:70
const_iterator
intermediate_table::const_iterator const_iterator
Definition: intermediate_table.cc:28
gar::geo::GeometryCore::GetCellID
gar::raw::CellID_t GetCellID(const TGeoNode *node, const unsigned int &det_id, const unsigned int &stave, const unsigned int &module, const unsigned int &layer, const unsigned int &slice, const std::array< double, 3 > &localPosition) const
Definition: GeometryCore.cxx:805
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
simb::MCTruth::NParticles
int NParticles() const
Definition: MCTruth.h:75
reco_momentum_tuples.px
px
Definition: reco_momentum_tuples.py:218
genie::units::MeV
static constexpr double MeV
Definition: Units.h:129
gar::geo::seg::MinervaSegmentationAlg::AddHitsMinerva
void AddHitsMinerva(std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > &m_Deposits, std::vector< gar::sdp::CaloDeposit > &fDeposits) const
Definition: MinervaSegmentationAlg.cxx:303
simb::MCParticle::Process
std::string Process() const
Definition: MCParticle.h:215
util::ConvertEdep2Art::fParticleList
sim::ParticleList * fParticleList
Definition: ConvertEdep2Art_module.cc:151
util::ConvertEdep2Art::fMinervaSegAlg
const gar::geo::seg::MinervaSegmentationAlg * fMinervaSegAlg
Definition: ConvertEdep2Art_module.cc:172
genie::GHepParticle::Pz
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
genie::GHepParticle::Status
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
gar::geo::LocalTransformation
Class to transform between world and local coordinates.
Definition: LocalTransformation.h:41
simb::MCParticle::TrackId
int TrackId() const
Definition: MCParticle.h:210
util::ConvertEdep2Art::fGeo
const gar::geo::GeometryCore * fGeo
geometry information
Definition: ConvertEdep2Art_module.cc:148
util::ConvertEdep2Art::FindMCTruthIndex
size_t FindMCTruthIndex(std::vector< simb::MCTruth > *mctruthcol, const simb::MCParticle &part) const
Definition: ConvertEdep2Art_module.cc:470
genie::GHepParticle::Px
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
genie::GHepParticle::LastMother
int LastMother(void) const
Definition: GHepParticle.h:67
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::GHepParticle::FirstMother
int FirstMother(void) const
Definition: GHepParticle.h:66
genie::GHepParticle::Name
string Name(void) const
Name that corresponds to the PDG code.
Definition: GHepParticle.cxx:135
search.process
def process(f, kind)
Definition: search.py:254
util::ConvertEdep2Art::fTrkIDParent
std::map< unsigned int, unsigned int > fTrkIDParent
Definition: ConvertEdep2Art_module.cc:167
reco_momentum_tuples.h
h
Definition: reco_momentum_tuples.py:66
genie::GHepParticle::LastDaughter
int LastDaughter(void) const
Definition: GHepParticle.h:69
util::ConvertEdep2Art::fMCRec
genie::NtpMCEventRecord * fMCRec
Definition: ConvertEdep2Art_module.cc:139
push_back
fInnerVessel push_back(Point(-578.400000, 0.000000, 0.000000))
e
const double e
Definition: gUpMuFluxGen.cxx:165
TG4Event::SegmentDetectors
TG4HitSegmentDetectors SegmentDetectors
Definition: TG4Event.h:39
util::ConvertEdep2Art::GetParentage
unsigned int GetParentage(unsigned int trkid) const
Definition: ConvertEdep2Art_module.cc:416
TG4HitSegment::GetPrimaryId
int GetPrimaryId() const
Definition: TG4HitSegment.h:40
genie::units::GeV
static constexpr double GeV
Definition: Units.h:28
gar::FillGTruth
void FillGTruth(const simb::GTruth &gt, garana::GTruth &outtruth)
Definition: AnaUtils.cxx:151
simb::kSingleParticle
single particles thrown at the detector
Definition: MCTruth.h:26
reco_momentum_tuples.t
t
Definition: reco_momentum_tuples.py:25
wirecell.gen.depos.move
def move(depos, offset)
Definition: depos.py:107
gar::util::FindProcessName
std::string FindProcessName(int process, int subprocess)
Definition: ProcessTable.hh:107
util::ConvertEdep2Art::fEMShowerDaughterMatRegex
std::string fEMShowerDaughterMatRegex
keep EM shower daughters only in these materials
Definition: ConvertEdep2Art_module.cc:145
test.p
p
Definition: test.py:223
util::ConvertEdep2Art::GetLayerNumber
unsigned int GetLayerNumber(std::string volname) const
Definition: ConvertEdep2Art_module.cc:344
art::DataViewImpl::put
ProductID put(std::unique_ptr< PROD > &&edp, std::string const &instance={})
Definition: DataViewImpl.h:686
util::CreateAssn
bool CreateAssn(PRODUCER const &prod, art::Event &evt, std::vector< T > const &a, art::Ptr< U > const &b, art::Assns< U, T > &assn, std::string a_instance, size_t indx=UINT_MAX)
Creates a single one-to-one association.
Definition: AssociationUtil.h:601
gar::raw::CellID_t
long long int CellID_t
Definition: CaloRawDigit.h:24
simb::MCTrajectory::Position
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
Definition: MCTrajectory.cxx:36
MF_LOG_INFO
#define MF_LOG_INFO(category)
Definition: MessageLogger.h:225
util::ConvertEdep2Art::GetDetNumber
unsigned int GetDetNumber(std::string volname) const
Definition: ConvertEdep2Art_module.cc:354
util::ConvertEdep2Art::fGArDeposits
std::vector< gar::sdp::EnergyDeposit > fGArDeposits
Definition: ConvertEdep2Art_module.cc:159
hit
Detector simulation of raw signals on wires.
Definition: NumberOfHitsFilter_module.cc:25
util::ConvertEdep2Art::GetSliceNumber
unsigned int GetSliceNumber(std::string volname) const
Definition: ConvertEdep2Art_module.cc:349
util::ConvertEdep2Art::GetModuleNumber
unsigned int GetModuleNumber(std::string volname) const
Definition: ConvertEdep2Art_module.cc:369
simb::MCParticle::SetEndProcess
void SetEndProcess(std::string s)
Definition: MCParticle.cxx:105
TG4HitSegment::GetTrackLength
double GetTrackLength() const
Definition: TG4HitSegment.h:59
util::ConvertEdep2Art::VisibleEnergyDeposition
double VisibleEnergyDeposition(const TG4HitSegment *hit, bool applyBirks) const
Definition: ConvertEdep2Art_module.cc:374
genie::EventRecord
Generated Event Record. It is a GHepRecord object that can accept / be visited by EventRecordVisitorI...
Definition: EventRecord.h:37
util::ConvertEdep2Art::m_MuIDDeposits
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_MuIDDeposits
Definition: ConvertEdep2Art_module.cc:156
util::ConvertEdep2Art::fMuIDDeposits
std::vector< gar::sdp::CaloDeposit > fMuIDDeposits
Definition: ConvertEdep2Art_module.cc:161
create_imaplist.material_name
material_name
Definition: create_imaplist.py:40
util::ConvertEdep2Art::fStartSpill
std::vector< int > fStartSpill
Definition: ConvertEdep2Art_module.cc:164
util::ConvertEdep2Art::fECALDeposits
std::vector< gar::sdp::CaloDeposit > fECALDeposits
Definition: ConvertEdep2Art_module.cc:158
reco_momentum_tuples.pz
pz
Definition: reco_momentum_tuples.py:220
simb::MCTruth::Add
void Add(simb::MCParticle const &part)
Definition: MCTruth.h:80
util::ConvertEdep2Art::CheckProcess
bool CheckProcess(std::string process_name) const
Definition: ConvertEdep2Art_module.cc:406
MF_LOG_DEBUG
#define MF_LOG_DEBUG(id)
Definition: MessageLogger.h:268
art::EventNumber_t
IDNumber_t< Level::Event > EventNumber_t
Definition: IDNumber.h:118
run_cnn_1class.pdg
pdg
Definition: run_cnn_1class.py:54
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
train.x
list x
Definition: train.py:276
util::ConvertEdep2Art::m_ECALDeposits
std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_ECALDeposits
Definition: ConvertEdep2Art_module.cc:154
util::ConvertEdep2Art::AddHits
void AddHits(std::map< gar::raw::CellID_t, std::vector< gar::sdp::CaloDeposit > > m_Deposits, std::vector< gar::sdp::CaloDeposit > &fDeposits)
Definition: ConvertEdep2Art_module.cc:432
genie::NtpMCRecHeader
MINOS-style Ntuple Class to hold an MC Event Record Header.
Definition: NtpMCRecHeader.h:32
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
simb::MCTruth
Event generator information.
Definition: MCTruth.h:32
reco_momentum_tuples.momentum
def momentum(x1, x2, x3, scale=1.)
Definition: reco_momentum_tuples.py:198
TG4HitSegment::GetStop
const TLorentzVector & GetStop() const
The stopping position of the segment.
Definition: TG4HitSegment.h:65
ns
QAsciiDict< Entry > ns
Definition: tclscanner.cpp:1147
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
s
static QCString * s
Definition: config.cpp:1042
art::Event::id
EventID id() const
Definition: Event.cc:34
genie::NtpMCEventRecord::event
EventRecord * event
event
Definition: NtpMCEventRecord.h:55
fhicl::exception
cet::coded_exception< error, detail::translate > exception
Definition: exception.h:33
util::ConvertEdep2Art::fTrackerDeposits
std::vector< gar::sdp::CaloDeposit > fTrackerDeposits
Definition: ConvertEdep2Art_module.cc:160
event
Event finding and building.
Definition: EventCheater_module.cc:32
TG4HitSegment::GetStart
const TLorentzVector & GetStart() const
The starting position of the segment.
Definition: TG4HitSegment.h:62
genie::GHepParticle::Py
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
std::map<int, size_t> util::ConvertEdep2Art::TrackIDToMCTruthIndexMap ( ) const
inlineprivate

Definition at line 116 of file ConvertEdep2Art_module.cc.

116 { return fTrackIDToMCTruthIndex; }
util::ConvertEdep2Art::fTrackIDToMCTruthIndex
std::map< int, size_t > fTrackIDToMCTruthIndex
Definition: ConvertEdep2Art_module.cc:168
double util::ConvertEdep2Art::VisibleEnergyDeposition ( const TG4HitSegment *  hit,
bool  applyBirks 
) const
private

Definition at line 374 of file ConvertEdep2Art_module.cc.

374  {
375  double edep = hit->GetEnergyDeposit();//MeV
376  double niel = hit->GetSecondaryDeposit();//MeV
377  double length = hit->GetTrackLength();//mm
378  double evis = edep;
379 
380  if(applyBirks)
381  {
382  if(fBirksCoeff > 0)
383  {
384  double nloss = niel;
385  if(nloss < 0.0) nloss = 0.0;
386  double eloss = edep - nloss;
387 
388  if(eloss < 0.0 || length <= 0.0) {
389  nloss = edep;
390  eloss = 0.0;
391  }
392  if(eloss > 0.0) { eloss /= (1.0 + fBirksCoeff*eloss/length); }
393 
394  if(nloss > 0){
395  nloss /= (1.0 + fBirksCoeff*nloss/length);
396  }
397 
398  evis = eloss + nloss;
399  }
400  }
401 
402  return evis;//MeV
403  }
TG4HitSegment::GetEnergyDeposit
double GetEnergyDeposit() const
The total energy deposit in this hit.
Definition: TG4HitSegment.h:43
TG4HitSegment::GetTrackLength
double GetTrackLength() const
Definition: TG4HitSegment.h:59
TG4HitSegment::GetSecondaryDeposit
double GetSecondaryDeposit() const
Definition: TG4HitSegment.h:54

Member Data Documentation

bool util::ConvertEdep2Art::fApplyBirks
private

Definition at line 127 of file ConvertEdep2Art_module.cc.

double util::ConvertEdep2Art::fBirksCoeff
private

Definition at line 142 of file ConvertEdep2Art_module.cc.

std::vector<gar::sdp::CaloDeposit> util::ConvertEdep2Art::fECALDeposits
private

Definition at line 158 of file ConvertEdep2Art_module.cc.

std::string util::ConvertEdep2Art::fECALMaterial
private

Definition at line 124 of file ConvertEdep2Art_module.cc.

const gar::detinfo::ECALProperties* util::ConvertEdep2Art::fEcalProp
private

Definition at line 149 of file ConvertEdep2Art_module.cc.

std::string util::ConvertEdep2Art::fEDepSimfile
private

Definition at line 121 of file ConvertEdep2Art_module.cc.

std::string util::ConvertEdep2Art::fEMShowerDaughterMatRegex
private

keep EM shower daughters only in these materials

Definition at line 145 of file ConvertEdep2Art_module.cc.

double util::ConvertEdep2Art::fEnergyCut
private

Definition at line 126 of file ConvertEdep2Art_module.cc.

TG4Event* util::ConvertEdep2Art::fEvent
private

Definition at line 141 of file ConvertEdep2Art_module.cc.

gar::geo::BitFieldCoder* util::ConvertEdep2Art::fFieldDecoderTrk
private

Definition at line 171 of file ConvertEdep2Art_module.cc.

std::vector<gar::sdp::EnergyDeposit> util::ConvertEdep2Art::fGArDeposits
private

Definition at line 159 of file ConvertEdep2Art_module.cc.

const gar::geo::GeometryCore* util::ConvertEdep2Art::fGeo
private

geometry information

Definition at line 148 of file ConvertEdep2Art_module.cc.

std::string util::ConvertEdep2Art::fGhepfile
private

Definition at line 122 of file ConvertEdep2Art_module.cc.

unsigned int util::ConvertEdep2Art::fGoodSpills
private

Total of good spills.

Definition at line 132 of file ConvertEdep2Art_module.cc.

TChain* util::ConvertEdep2Art::fGTreeChain
private

Definition at line 137 of file ConvertEdep2Art_module.cc.

bool util::ConvertEdep2Art::fHasGHEP
private

Definition at line 152 of file ConvertEdep2Art_module.cc.

bool util::ConvertEdep2Art::fkeepEMShowers
private

Definition at line 144 of file ConvertEdep2Art_module.cc.

genie::NtpMCEventRecord* util::ConvertEdep2Art::fMCRec
private

Definition at line 139 of file ConvertEdep2Art_module.cc.

const gar::geo::seg::MinervaSegmentationAlg* util::ConvertEdep2Art::fMinervaSegAlg
private

Definition at line 172 of file ConvertEdep2Art_module.cc.

std::vector<gar::sdp::CaloDeposit> util::ConvertEdep2Art::fMuIDDeposits
private

Definition at line 161 of file ConvertEdep2Art_module.cc.

bool util::ConvertEdep2Art::fOverlay
private

Definition at line 123 of file ConvertEdep2Art_module.cc.

sim::ParticleList* util::ConvertEdep2Art::fParticleList
private

Definition at line 151 of file ConvertEdep2Art_module.cc.

int util::ConvertEdep2Art::fSpillCount
private

Definition at line 163 of file ConvertEdep2Art_module.cc.

unsigned int util::ConvertEdep2Art::fSpills
private

Total of spills.

Definition at line 131 of file ConvertEdep2Art_module.cc.

std::vector<int> util::ConvertEdep2Art::fStartSpill
private

Definition at line 164 of file ConvertEdep2Art_module.cc.

std::vector<int> util::ConvertEdep2Art::fStopSpill
private

Definition at line 165 of file ConvertEdep2Art_module.cc.

double util::ConvertEdep2Art::fTotGoodPOT
private

Total good POT.

Definition at line 130 of file ConvertEdep2Art_module.cc.

double util::ConvertEdep2Art::fTotPOT
private

Total POT.

Definition at line 129 of file ConvertEdep2Art_module.cc.

std::string util::ConvertEdep2Art::fTPCMaterial
private

Definition at line 125 of file ConvertEdep2Art_module.cc.

std::vector<gar::sdp::CaloDeposit> util::ConvertEdep2Art::fTrackerDeposits
private

Definition at line 160 of file ConvertEdep2Art_module.cc.

std::map<int, size_t> util::ConvertEdep2Art::fTrackIDToMCTruthIndex
private

Definition at line 168 of file ConvertEdep2Art_module.cc.

TChain* util::ConvertEdep2Art::fTreeChain
private

Definition at line 134 of file ConvertEdep2Art_module.cc.

std::map<unsigned int, unsigned int> util::ConvertEdep2Art::fTrkIDParent
private

Definition at line 167 of file ConvertEdep2Art_module.cc.

std::map< gar::raw::CellID_t, std::vector<gar::sdp::CaloDeposit> > util::ConvertEdep2Art::m_ECALDeposits
private

Definition at line 154 of file ConvertEdep2Art_module.cc.

std::map< gar::raw::CellID_t, std::vector<gar::sdp::CaloDeposit> > util::ConvertEdep2Art::m_MuIDDeposits
private

Definition at line 156 of file ConvertEdep2Art_module.cc.

std::map< gar::raw::CellID_t, std::vector<gar::sdp::CaloDeposit> > util::ConvertEdep2Art::m_TrackerDeposits
private

Definition at line 155 of file ConvertEdep2Art_module.cc.

unsigned int util::ConvertEdep2Art::nEntries
private

Definition at line 135 of file ConvertEdep2Art_module.cc.

unsigned int util::ConvertEdep2Art::nEntriesGhep
private

Definition at line 138 of file ConvertEdep2Art_module.cc.

genie::PDGLibrary* util::ConvertEdep2Art::pdglib
private

Definition at line 147 of file ConvertEdep2Art_module.cc.

The documentation for this class was generated from the following file:
Generated by   doxygen 1.8.11