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/Framework/Services/Optional/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  fGeantLabel(pset.get< std::string >("G4ModuleLabel")) //from TrueEnergy_module.cc
73  {
74  /*
75  art::ServiceHandle<art::TFileService> tfs;
76  fMyTree = tfs->make<TTree>("ContTree", "ContTree");
77  fMyTree->Branch( "PartEn" , &fPartEn );
78  fMyTree->Branch( "TrackID" , &fTrackID );
79  fMyTree->Branch( "PDGCode" , &fPDGCode );
80 
81  fMyTree->Branch( "StartX" , &fStartX );
82  fMyTree->Branch( "StartY" , &fStartY );
83  fMyTree->Branch( "StartZ" , &fStartZ );
84  fMyTree->Branch( "EndX" , &fEndX );
85  fMyTree->Branch( "EndY" , &fEndY );
86  fMyTree->Branch( "EndZ" , &fEndZ );
87 
88  fMyTree->Branch( "MomX" , &fMomX );
89  fMyTree->Branch( "MomY" , &fMomY );
90  fMyTree->Branch( "MomZ" , &fMomZ );
91  fMyTree->Branch( "NormMomX", &fNormMomX );
92  fMyTree->Branch( "NormMomY", &fNormMomY );
93  fMyTree->Branch( "NormMomZ", &fNormMomZ );
94  fMyTree->Branch( "TotMom" , &fTotMom );
95 
96  fMyTree->Branch( "EnterEn" , &fEnterEn );
97  fMyTree->Branch( "EscapeEn", &fEscapeEn );
98  fMyTree->Branch( "TotEscEn", &fTotEscEn );
99  fMyTree->Branch( "TotDepEn", &fTotDepEn );
100  fMyTree->Branch( "CorEscEn", &fCorEscEn );
101  */
102  }
103 
104  //......................................................................
106  /*
107  fPartEn = fTrackID = fPDGCode = -5;
108  fStartX = fStartY = fStartZ = -5;
109  fEndX = fEndY = fEndZ = -5;
110  fNormMomX= fNormMomY= fNormMomZ = -5;
111  fMomX = fMomY = fMomZ = fTotMom = -5;
112  fEnterEn = fEscapeEn= fTotEscEn = fTotDepEn = fCorEscEn = -5;
113  */
114  }
115 
116  //......................................................................
118  {
119  ResetTree();
121  evt.getByLabel(fGeantLabel, firstParticles);
122 
123  art::Handle< std::vector<sim::TrueEnergy> > firstTrueEnergies;
124  evt.getByLabel(fGeantLabel, firstTrueEnergies);
125 
127 
128  bool fPass = false;
129  for(unsigned int i = 0; i < firstTrueEnergies->size(); ++i){
130  sim::TrueEnergy trueE = ((*firstTrueEnergies)[i]);
131  if (trueE.TrackId() == 1) {
132  const sim::Particle *TruePart = bt->TrackIDToParticle(trueE.TrackId());
133 
134  // -- Make some Vars for the TrueEnergy::EscapingEnergy and the End Position of the Primary.
135  float CorEscEn = trueE.EscapingEnergy();
136  float EndX = TruePart->EndX();
137  float EndY = TruePart->EndY();
138  float EndZ = TruePart->EndZ();
139  // -- Loop through the daughters.
140  for (int par=0; par<TruePart->NumberDaughters(); ++par) {
141  const sim::Particle *Daught = bt->TrackIDToParticle(TruePart->Daughter(par));
142  // -- Loop through the TrueEnergies.
143  for(unsigned int te = 0; te < firstTrueEnergies->size(); ++te){
144  // -- If found the TrueEnergy
145  sim::TrueEnergy DaughtTE = ((*firstTrueEnergies)[te]);
146  if ( DaughtTE.TrackId() == Daught->TrackId() ) {
147  float DaughtStX = Daught->Vx();
148  float DaughtStY = Daught->Vy();
149  float DaughtStZ = Daught->Vz();
150  // -- Check that it is a hard scatter.
151  if ( Daught->PdgCode() == TruePart->PdgCode() &&
152  fabs( DaughtStX - EndX ) < 1.0 &&
153  fabs( DaughtStY - EndY ) < 1.0 &&
154  fabs( DaughtStZ - EndZ ) < 1.0 ) {
155  CorEscEn += DaughtTE.EscapingEnergy();
156  }
157  break;
158  } // if hard-scatter
159  } // loop over te
160  } // loop over par
161 
162  // Do the cut - for rational see DocDB 45934.
163  if ( trueE.EnteringEnergy() > 0.0 && CorEscEn < 0.1 ) {
164  fPass = true;
165  }
166 
167  // Fill MyTree.
168  /*
169  fPartEn = TruePart->E();
170  fTrackID = TruePart->TrackId();
171  fPDGCode = TruePart->PdgCode();
172  // Pos
173  fStartX = TruePart->Vx();
174  fStartY = TruePart->Vy();
175  fStartZ = TruePart->Vz();
176  fEndX = TruePart->EndX();
177  fEndY = TruePart->EndY();
178  fEndZ = TruePart->EndZ();
179  // Mom
180  fMomX = TruePart->Px();
181  fMomY = TruePart->Py();
182  fMomZ = TruePart->Pz();
183  fTotMom = TruePart->P();
184  TVector3 UnitMom = TruePart->Momentum().Vect().Unit();
185  fNormMomX = UnitMom.X();
186  fNormMomY = UnitMom.Y();
187  fNormMomZ = UnitMom.Z();
188  // sim::TrueEnergy
189  fEnterEn = trueE.EnteringEnergy();
190  fEscapeEn = trueE.EscapingEnergy();
191  fTotEscEn = trueE.TotalEscEnergy();
192  fTotDepEn = trueE.TotalDepEnergy();
193  fCorEscEn = CorEscEn;
194  // Fill the tree
195  fMyTree->Fill();
196 
197  std::cout << "\n Primary particle with TrackId " << trueE.TrackId() << " is a " << TruePart->PdgCode() << " with " << TruePart->NumberDaughters() << " daughters."
198  << " It had StartingE " << TruePart->E() << ", and EndE " << TruePart->E(TruePart->NumberTrajectoryPoints()-1)
199  << "\n\t It finished at position " << fEndX << ", " << fEndY << ", " << fEndZ
200  << "\n\t It had EnterE " << trueE.EnteringEnergy() << ", EscapeE " << trueE.EscapingEnergy() << ", TotEscE " << trueE.TotalEscEnergy()
201  << std::endl;
202  */
203 
204  // No point looping all the way through the TrueEnergy vector
205  break;
206  } // if Primary
207  } // loop through TrueEnergy
208  return fPass;
209  } // filter
210 
212 } // 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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
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
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
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
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