TruthContain_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // \file TruthContain_module.cc
3 //\brief Filter to remove events that are not suficiently contained within the detector
4 
5 // \author Beatriz Tapia - btapiaoregui@utexas.edu & Karl Warburton karlwarb@iastate.edu
6 ////////////////////////////////////////////////////////////////////////
7 
8 // NOvA includes
10 
11 // Framework includes
14 //#include "art_root_io/TFileService.h"
17 #include "fhiclcpp/ParameterSet.h"
19 
21 #include "MCCheater/BackTracker.h"
22 #include "Simulation/Particle.h"
23 #include "Simulation/TrueEnergy.h"
24 #include "Utilities/AssociationUtil.h"
27 #include <memory>
28 
29 // ROOT includes
30 #include "TDatabasePDG.h"
31 #include "TTree.h"
32 #include "TH1F.h"
33 
34 
35 namespace filter
36 {
37  class TruthContain;
38 
39  class TruthContain : public art::EDFilter {
40  public:
41  explicit TruthContain(const fhicl::ParameterSet & pset);
42 
43  // Plugins should not be copied or assigned.
44  TruthContain(TruthContain const &) = delete;
45  TruthContain(TruthContain &&) = delete;
46  TruthContain & operator = (TruthContain const &) = delete;
47  TruthContain & operator = (TruthContain &&) = delete;
48 
49  // Required functions.
50  bool filter(art::Event & e) override;
51 
52  // Added Functions
53  void ResetTree();
54 
55  private:
56  std::string fGeantLabel; ///< label for module running G4 and making particles, etc
57 
58  // My TTree
59  /*
60  TTree* fMyTree;
61  float fPartEn , fTrackID , fPDGCode;
62  float fStartX , fStartY , fStartZ;
63  float fEndX , fEndY , fEndZ;
64  float fNormMomX, fNormMomY, fNormMomZ;
65  float fMomX , fMomY , fMomZ , fTotMom;
66  float fEnterEn , fEscapeEn, fTotEscEn, fTotDepEn, fCorEscEn;
67  */
68  };
69 
70  //......................................................................
72  EDFilter(pset),
73  fGeantLabel(pset.get< std::string >("G4ModuleLabel")) //from TrueEnergy_module.cc
74  {
75  /*
76  art::ServiceHandle<art::TFileService> tfs;
77  fMyTree = tfs->make<TTree>("ContTree", "ContTree");
78  fMyTree->Branch( "PartEn" , &fPartEn );
79  fMyTree->Branch( "TrackID" , &fTrackID );
80  fMyTree->Branch( "PDGCode" , &fPDGCode );
81 
82  fMyTree->Branch( "StartX" , &fStartX );
83  fMyTree->Branch( "StartY" , &fStartY );
84  fMyTree->Branch( "StartZ" , &fStartZ );
85  fMyTree->Branch( "EndX" , &fEndX );
86  fMyTree->Branch( "EndY" , &fEndY );
87  fMyTree->Branch( "EndZ" , &fEndZ );
88 
89  fMyTree->Branch( "MomX" , &fMomX );
90  fMyTree->Branch( "MomY" , &fMomY );
91  fMyTree->Branch( "MomZ" , &fMomZ );
92  fMyTree->Branch( "NormMomX", &fNormMomX );
93  fMyTree->Branch( "NormMomY", &fNormMomY );
94  fMyTree->Branch( "NormMomZ", &fNormMomZ );
95  fMyTree->Branch( "TotMom" , &fTotMom );
96 
97  fMyTree->Branch( "EnterEn" , &fEnterEn );
98  fMyTree->Branch( "EscapeEn", &fEscapeEn );
99  fMyTree->Branch( "TotEscEn", &fTotEscEn );
100  fMyTree->Branch( "TotDepEn", &fTotDepEn );
101  fMyTree->Branch( "CorEscEn", &fCorEscEn );
102  */
103  }
104 
105  //......................................................................
107  /*
108  fPartEn = fTrackID = fPDGCode = -5;
109  fStartX = fStartY = fStartZ = -5;
110  fEndX = fEndY = fEndZ = -5;
111  fNormMomX= fNormMomY= fNormMomZ = -5;
112  fMomX = fMomY = fMomZ = fTotMom = -5;
113  fEnterEn = fEscapeEn= fTotEscEn = fTotDepEn = fCorEscEn = -5;
114  */
115  }
116 
117  //......................................................................
119  {
120  ResetTree();
122  evt.getByLabel(fGeantLabel, firstParticles);
123 
124  art::Handle< std::vector<sim::TrueEnergy> > firstTrueEnergies;
125  evt.getByLabel(fGeantLabel, firstTrueEnergies);
126 
128 
129  bool fPass = false;
130  for(unsigned int i = 0; i < firstTrueEnergies->size(); ++i){
131  sim::TrueEnergy trueE = ((*firstTrueEnergies)[i]);
132  if (trueE.TrackId() == 1) {
133  const sim::Particle *TruePart = bt->TrackIDToParticle(trueE.TrackId());
134 
135  // -- Make some Vars for the TrueEnergy::EscapingEnergy and the End Position of the Primary.
136  float CorEscEn = trueE.EscapingEnergy();
137  float EndX = TruePart->EndX();
138  float EndY = TruePart->EndY();
139  float EndZ = TruePart->EndZ();
140  // -- Loop through the daughters.
141  for (int par=0; par<TruePart->NumberDaughters(); ++par) {
142  const sim::Particle *Daught = bt->TrackIDToParticle(TruePart->Daughter(par));
143  // -- Loop through the TrueEnergies.
144  for(unsigned int te = 0; te < firstTrueEnergies->size(); ++te){
145  // -- If found the TrueEnergy
146  sim::TrueEnergy DaughtTE = ((*firstTrueEnergies)[te]);
147  if ( DaughtTE.TrackId() == Daught->TrackId() ) {
148  float DaughtStX = Daught->Vx();
149  float DaughtStY = Daught->Vy();
150  float DaughtStZ = Daught->Vz();
151  // -- Check that it is a hard scatter.
152  if ( Daught->PdgCode() == TruePart->PdgCode() &&
153  fabs( DaughtStX - EndX ) < 1.0 &&
154  fabs( DaughtStY - EndY ) < 1.0 &&
155  fabs( DaughtStZ - EndZ ) < 1.0 ) {
156  CorEscEn += DaughtTE.EscapingEnergy();
157  }
158  break;
159  } // if hard-scatter
160  } // loop over te
161  } // loop over par
162 
163  // Do the cut - for rational see DocDB 45934.
164  if ( trueE.EnteringEnergy() > 0.0 && CorEscEn < 0.1 ) {
165  fPass = true;
166  }
167 
168  // Fill MyTree.
169  /*
170  fPartEn = TruePart->E();
171  fTrackID = TruePart->TrackId();
172  fPDGCode = TruePart->PdgCode();
173  // Pos
174  fStartX = TruePart->Vx();
175  fStartY = TruePart->Vy();
176  fStartZ = TruePart->Vz();
177  fEndX = TruePart->EndX();
178  fEndY = TruePart->EndY();
179  fEndZ = TruePart->EndZ();
180  // Mom
181  fMomX = TruePart->Px();
182  fMomY = TruePart->Py();
183  fMomZ = TruePart->Pz();
184  fTotMom = TruePart->P();
185  TVector3 UnitMom = TruePart->Momentum().Vect().Unit();
186  fNormMomX = UnitMom.X();
187  fNormMomY = UnitMom.Y();
188  fNormMomZ = UnitMom.Z();
189  // sim::TrueEnergy
190  fEnterEn = trueE.EnteringEnergy();
191  fEscapeEn = trueE.EscapingEnergy();
192  fTotEscEn = trueE.TotalEscEnergy();
193  fTotDepEn = trueE.TotalDepEnergy();
194  fCorEscEn = CorEscEn;
195  // Fill the tree
196  fMyTree->Fill();
197 
198  std::cout << "\n Primary particle with TrackId " << trueE.TrackId() << " is a " << TruePart->PdgCode() << " with " << TruePart->NumberDaughters() << " daughters."
199  << " It had StartingE " << TruePart->E() << ", and EndE " << TruePart->E(TruePart->NumberTrajectoryPoints()-1)
200  << "\n\t It finished at position " << fEndX << ", " << fEndY << ", " << fEndZ
201  << "\n\t It had EnterE " << trueE.EnteringEnergy() << ", EscapeE " << trueE.EscapingEnergy() << ", TotEscE " << trueE.TotalEscEnergy()
202  << std::endl;
203  */
204 
205  // No point looping all the way through the TrueEnergy vector
206  break;
207  } // if Primary
208  } // loop through TrueEnergy
209  return fPass;
210  } // filter
211 
213 } // namespace
std::string fGeantLabel
label for module running G4 and making particles, etc
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int TrackId() const
Definition: TrueEnergy.h:36
double EndZ() const
Definition: MCParticle.h:227
Int_t par
Definition: SimpleIterate.C:24
TruthContain & operator=(TruthContain const &)=delete
TruthContain(const fhicl::ParameterSet &pset)
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
double EndY() const
Definition: MCParticle.h:226
float EscapingEnergy() const
Definition: TrueEnergy.h:38
int NumberDaughters() const
Definition: MCParticle.h:216
Module that kips a configurable number of events between each that it allows through. Note that this module really skips (N-1) events, it uses a simple modular division as its critera. This module will cut down the data sample to 1/N of its original size.
float EnteringEnergy() const
Definition: TrueEnergy.h:37
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
bool getByLabel(std::string const &label, std::string const &instance, Handle< PROD > &result) const
Definition: DataViewImpl.h:446
int evt
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
double Vx(const int i=0) const
Definition: MCParticle.h:220
EDFilter(fhicl::ParameterSet const &pset)
Definition: EDFilter.h:22
double Vz(const int i=0) const
Definition: MCParticle.h:222
bool filter(art::Event &e) override
double EndX() const
Definition: MCParticle.h:225
Float_t e
Definition: plot.C:35
double Vy(const int i=0) const
Definition: MCParticle.h:221
enum BeamMode string