GenieTruth_module.cc
Go to the documentation of this file.
1 ///
2 /// \brief
3 /// \author Xuebing Bu - xbbu@fnal.gov
4 ///
6 
7 #include "Geometry/Geometry.h"
9 
10 // ROOT includes
11 #include "TVector3.h"
12 #include "TH1F.h"
13 #include "TMath.h"
14 #include "TTree.h"
15 #include "TProfile2D.h"
16 #include "time.h"
17 #include "TProfile.h"
18 
22 #include "Simulation/Particle.h"
24 #include "SummaryData/POTSum.h"
25 #include "SummaryData/SpillData.h"
26 #include "RecoBase/CellHit.h"
27 #include "RecoBase/Cluster.h"
28 #include "RecoBase/RecoHit.h"
29 #include "RecoBase/Track.h"
30 #include "RecoBase/Vertex.h"
31 #include "RecoBase/Prong.h"
32 #include "RecoBase/HoughResult.h"
33 
34 #include "Utilities/AssociationUtil.h"
35 #include "MCCheater/BackTracker.h"
36 #include "Utilities/func/MathUtil.h"
37 
38 // Framework includes
47 #include "fhiclcpp/ParameterSet.h"
53 
54 #include <cmath>
55 #include <vector>
56 #include <string>
57 #include <vector>
58 
59 //include TMVA
60 #include "TMVA/Tools.h"
61 #include "TMVA/Reader.h"
62 #include "TMVA/MethodCuts.h"
63 
64 class TVector3;
65 class TH1F;
66 class TTree;
67 
68 namespace rb{class Cluster; class Track;}
69 namespace sim{class Particle;}
70 namespace simb{class MCFlux; class MCTruth; class MCNeutrino;}
71 
72 namespace ncs {
73  class GenieTruth : public art::EDAnalyzer {
74  public:
75  explicit GenieTruth(fhicl::ParameterSet const& pset);
76  virtual ~GenieTruth();
77 
78  void beginJob();
79  void analyze(const art::Event& evt);
80  void reconfigure(const fhicl::ParameterSet& p);
81  void endSubRun(const art::SubRun& sr);
82  void endJob();
83  TTree *_btree;
84 
85  private:
86 
87  const rb::RecoHit MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
88  double const& w);
89 
90  //folders to get objects from
91  std::string fMCFluxLabel; ///label for generator information
92  std::string fweightLabel; ///genie weights
93 
94  //switch to print the decay chain
96 
97  //spill-level info
98  double spillPOT;
99 
100  TTree *_otree;
101 
102  /*** true info ***/
103  float nu_Vx;//true vertex X
104  float nu_Vy;//true vertex Y
105  float nu_Vz;//true vertex Z
106  float nu_truePx;//true Px
107  float nu_truePy;//true Py
108  float nu_truePz;//true Pz
109  float nu_trueE;//true E
110  int nu_ccnc;
111  int nu_PDG;
114  int nu_mode;
117  float nu_HadX;
118  float nu_HadY;
119  float nu_HadW;
120  float nu_Q2;
121  //lepton info
122  int lep_PDG;
123  float lep_Px;
124  float lep_Py;
125  float lep_Pz;
126  float lep_E;
127 
128  //4-momentum for particles originating from the nu interaction vertex
130  int partPDG[100];
131  float partPx[100];
132  float partPy[100];
133  float partPz[100];
134  float partE[100];
135 
136  //GENIE reweights or cross section uncertainties
137  int Nweights;
138  float plus1sigma[68];
139  float plus2sigma[68];
140  float minus1sigma[68];
141  float minus2sigma[68];
142 
143  };
144 }
145 
146 namespace ncs{
147 
148  //......................................................................
149  GenieTruth::GenieTruth(fhicl::ParameterSet const& pset)
150  : EDAnalyzer(pset)
151  {
152  mf::LogInfo("GenieTruth")<<__PRETTY_FUNCTION__<<"\n";
153  reconfigure(pset);
154  }
155 
156  //......................................................................
158  {
159  }
160 
161  //......................................................................
163  {
164  fMCFluxLabel = pset.get< std::string >("MCFluxLabel");
165  fweightLabel = pset.get< std::string >("weightLabel");
166  printDecay = pset.get< int >("printdecay");
167  }
168 
169  //......................................................................
171  {
173 
174  //personal tuples
175  _otree = tfs->make<TTree>("GenieTruth","GenieTruth particle");
176 
177  _otree->Branch("nu_Vx",&nu_Vx,"nu_Vx/F");
178  _otree->Branch("nu_Vy",&nu_Vy,"nu_Vy/F");
179  _otree->Branch("nu_Vz",&nu_Vz,"nu_Vz/F");
180  _otree->Branch("nu_trueE",&nu_trueE,"nu_trueE/F");
181  _otree->Branch("nu_truePx",&nu_truePx,"nu_truePx/F");
182  _otree->Branch("nu_truePy",&nu_truePy,"nu_truePy/F");
183  _otree->Branch("nu_truePz",&nu_truePz,"nu_truePz/F");
184  _otree->Branch("nu_ccnc",&nu_ccnc,"nu_ccnc/I");
185  _otree->Branch("nu_PDG",&nu_PDG,"nu_PDG/I");
186  _otree->Branch("nu_orig_parentPDG",&nu_orig_parentPDG,"nu_orig_parentPDG/I");
187  _otree->Branch("nu_origPDG",&nu_origPDG,"nu_origPDG/I");
188  _otree->Branch("nu_mode",&nu_mode,"nu_mode/I");
189  _otree->Branch("nu_intType",&nu_intType,"nu_intType/I");
190  _otree->Branch("nu_hitNucl",&nu_hitNucl,"nu_hitNucl/I");
191  _otree->Branch("nu_HadX",&nu_HadX,"nu_HadX/F");
192  _otree->Branch("nu_HadY",&nu_HadY,"nu_HadY/F");
193  _otree->Branch("nu_HadW",&nu_HadW,"nu_HadW/F");
194  _otree->Branch("nu_Q2",&nu_Q2,"nu_Q2/F");
195  _otree->Branch("lep_PDG",&lep_PDG,"lep_PDG/I");
196  _otree->Branch("lep_Px",&lep_Px,"lep_Px/F");
197  _otree->Branch("lep_Py",&lep_Py,"lep_Py/F");
198  _otree->Branch("lep_Pz",&lep_Pz,"lep_Pz/F");
199  _otree->Branch("lep_E",&lep_E,"lep_E/F");
200 
201  _otree->Branch("Nparticles",&Nparticles,"Nparticles/I");
202  _otree->Branch("partPDG",&partPDG,"partPDG[Nparticles]/I");
203  _otree->Branch("partPx",&partPx,"partPx[Nparticles]/F");
204  _otree->Branch("partPy",&partPy,"partPy[Nparticles]/F");
205  _otree->Branch("partPz",&partPz,"partPz[Nparticles]/F");
206  _otree->Branch("partE",&partE,"partE[Nparticles]/F");
207 
208  _otree->Branch("Nweights",&Nweights,"Nweights/I");
209  _otree->Branch("plus1sigma",&plus1sigma,"plus1sigma[Nweights]/F");
210  _otree->Branch("plus2sigma",&plus2sigma,"plus2sigma[Nweights]/F");
211  _otree->Branch("minus1sigma",&minus1sigma,"minus1sigma[Nweights]/F");
212  _otree->Branch("minus2sigma",&minus2sigma,"minus2sigma[Nweights]/F");
213 
214  _btree=new TTree("spillInfo","spill level info");
215  _btree->Branch("spillPOT",&spillPOT,"spillPOT/D");
216  }
217 
218  //......................................................................
220  {
221  }
222 
223  //......................................................................
225  double const& w)
226  {
228  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
229 
230  return rhit;
231  }
232 
233  //......................................................................
235  {
236 
237  nu_Vx = -999.;
238  nu_Vy = -999.;
239  nu_Vz = -999.;
240  nu_trueE = -1.;
241  nu_truePx = -999.;
242  nu_truePy = -999.;
243  nu_truePz = -999.;
244 
245  nu_ccnc=-1;
246  nu_PDG=-1;
247  nu_origPDG=-1;
249  nu_mode=-1;
250  nu_intType=-1;
251  nu_hitNucl=-1;
252 
253  nu_HadX=-1.;
254  nu_HadY=-1.;
255  nu_HadW=-1.;
256  nu_Q2=-1.;
257 
258  lep_PDG=-1;
259  lep_Px=-1.;
260  lep_Py=-1.;
261  lep_Pz=-1.;
262  lep_E=-1.;
263 
264  Nparticles=0;
265  for( int i=0; i<100; ++i ){
266  partPDG[i]=-1;
267  partPx[i]=-999.;
268  partPy[i]=-999.;
269  partPz[i]=-999.;
270  partE[i]=-1.0;
271  }
272 
273  Nweights=0;
274  for( int i=0; i<68; ++i ){
275  plus1sigma[i]=1.;
276  plus2sigma[i]=1.;
277  minus1sigma[i]=1.;
278  minus2sigma[i]=1.;
279  }
280 
281  //art::ServiceHandle<calib::Calibrator> cal;
282  //art::ServiceHandle<geo::Geometry> geom;
284 
285  //beam intensity
287  evt.getByLabel(fMCFluxLabel, spillPot);
288  spillPOT=(spillPot->spillpot)/1.0e+12;
289 
290  _btree->Fill();//fill the tree for the spill-level info
291 
292  //////////////
293  //Truth info//
294  //////////////
296  evt.getByLabel(fMCFluxLabel,mclist);
298  evt.getByLabel(fMCFluxLabel, fllist);
299 
300  int Nneutrino=mclist->size();
301  //std::cout<<"# of nu interactions are "<<Nneutrino<<std::endl;
302 
303  if( Nneutrino>0 ){
304  for( int ii=0; ii<Nneutrino; ++ii ){
305 
306  art::Ptr<simb::MCTruth> mctruth(mclist,ii);
307  const simb::MCNeutrino& nu = mctruth->GetNeutrino();
308 
309  //remove the interactions outside the detector
310  if( fabs(nu.Nu().Vx())>200.0 ) continue;
311  if( fabs(nu.Nu().Vy())>200.0 ) continue;
312  if( nu.Nu().Vz()<0.0 ) continue;
313 
314  nu_trueE = nu.Nu().E();
315  nu_truePx = nu.Nu().Px();
316  nu_truePy = nu.Nu().Py();
317  nu_truePz = nu.Nu().Pz();
318 
319  nu_Vx=nu.Nu().Vx();
320  nu_Vy=nu.Nu().Vy();
321  nu_Vz=nu.Nu().Vz();
322 
323  nu_PDG = nu.Nu().PdgCode();
324  nu_ccnc = nu.CCNC();
325  nu_mode = nu.Mode();
327 
328  nu_hitNucl = nu.HitNuc();
329  nu_HadX = nu.X();
330  nu_HadY = nu.Y();
331  nu_HadW = nu.W();
332  nu_Q2 = nu.QSqr();
333 
334  lep_PDG = nu.Lepton().PdgCode();
335  lep_Px = nu.Lepton().Px();
336  lep_Py = nu.Lepton().Py();
337  lep_Pz = nu.Lepton().Pz();
338  lep_E = nu.Lepton().E();
339 
340  art::Ptr<simb::MCFlux> flux(fllist,ii);
341 
342  nu_origPDG = flux->fntype;
343  nu_orig_parentPDG = flux->fptype;
344 
345  //get genie cross section reweight table
347  pv.push_back(mctruth);
349 
350  const std::vector<art::Ptr<rwgt::GENIEReweightTable>> XSprods = gweights.at(0);
351  if( !XSprods.empty()){
352  const unsigned int N = XSprods[0]->NShifts();
353  if( N==68 ){
354  Nweights=N;
355  for(unsigned int ir = 0; ir < N; ++ir){
356  plus1sigma[ir]=XSprods[0]->Plus1Sigma(ir);
357  plus2sigma[ir]=XSprods[0]->Plus2Sigma(ir);
358  minus1sigma[ir]=XSprods[0]->Minus1Sigma(ir);
359  minus2sigma[ir]=XSprods[0]->Minus2Sigma(ir);
360  }
361  }//number of weights>0
362  }//having genie weights
363 
364  sim::ParticleNavigator const& pn = bt->ParticleNavigator();
365  int numberOfParticles = pn.size();
366 
367  std::vector<int> particlePDG;
368  std::vector<float> particlePx;
369  std::vector<float> particlePy;
370  std::vector<float> particlePz;
371  std::vector<float> particleE;
372  particlePDG.clear();
373  particlePx.clear();
374  particlePy.clear();
375  particlePz.clear();
376  //build a map for indexing and track ID
377  std::map<std::string,int> IndexMap;
378  IndexMap.clear();
379  for( int i=0; i<numberOfParticles; ++i ){
380  sim::Particle* particle = pn.Particle(i);
381  std::stringstream ss;//create a stringstream
382  ss << particle->TrackId();//add number to the stream
383  IndexMap[ss.str()] = i;
384 
385  double Dx=fabs(nu.Lepton().Vx()-particle->Vx());
386  double Dy=fabs(nu.Lepton().Vy()-particle->Vy());
387  double Dz=fabs(nu.Lepton().Vz()-particle->Vz());
388  if( particle->Mother()==0 && Dx<1.0E-6 && Dy<1.0E-6 && Dz<1.0E-6 ){
389  //std::cout<<particle->PdgCode()<<" "<<particle->E()<<" "<<particle->Mother()<<std::endl;
390  particlePDG.push_back(particle->PdgCode());
391  particleE.push_back(particle->E());
392  particlePx.push_back(particle->Px());
393  particlePy.push_back(particle->Py());
394  particlePz.push_back(particle->Pz());
395  }
396  }//end loop of particles
397 
398  int Npart=particlePDG.size();
399  if( Npart>100 ) Npart=100;
400  if( Npart>0 ){
401  Nparticles=Npart;
402  for( int i=0; i<Npart; ++i ){
403  partPDG[i]=particlePDG[i];
404  partE[i]=particleE[i];
405  partPx[i]=particlePx[i];
406  partPy[i]=particlePy[i];
407  partPz[i]=particlePz[i];
408  }//end loop of particlePDG
409  }//find the particles originating from the nu interaction vertex
410 
411  if( printDecay ){
412  std::cout<<"Incoming nu pdg is "<<flux->fntype<<std::endl;
413  std::cout<<"Outgoing nu pdg is "<<nu.Nu().PdgCode()<<std::endl;
414  std::cout<<"Interaction mode is "<<nu.Mode()<<std::endl;
415  std::cout<<"Is CCNC or not: "<<nu.CCNC()<<std::endl;
416  std::cout<<"# of particles = "<<numberOfParticles<<std::endl;
417  for( int i=0; i<numberOfParticles; ++i ){
418  sim::Particle* particle = pn.Particle(i);
419  int pdgid=particle->PdgCode();
420  double energy=particle->E();
421  int ndau=particle->NumberDaughters();
422  std::cout<<i<<" "<<pdgid<<" E="<<energy<<" ";
423 
424  //find mother
425  int MotherIndex=particle->Mother();
426  if( MotherIndex>0 ){
427  std::cout<<"mom idx: ";
428  std::stringstream ss;
429  ss << pn[MotherIndex]->TrackId();
430  std::map<std::string, int>::const_iterator mitr=IndexMap.find(ss.str());
431  if( mitr != IndexMap.end() )
432  std::cout<<mitr->second<<" ";
433  }//haveing moms
434 
435  //find daughter
436  if( ndau>0 ){
437  std::cout<<"dau idx: ";
438  for( int j=0; j<ndau; ++j ){
439  int dauID=particle->Daughter(j);
440  std::stringstream ss;
441  ss << pn[dauID]->TrackId();
442  std::map<std::string, int>::const_iterator mitr=IndexMap.find(ss.str());
443  if( mitr != IndexMap.end() )
444  std::cout<<mitr->second<<" ";
445  }//end loop of daughter particles
446  }//having daughters
447  std::cout<<"\n";
448  }//end loop of all particles
449  std::cout<<"=================================="<<std::endl;
450  }//end of printing decay chain
451 
452  _otree->Fill();
453  }//end loop of nu interactions
454  }//have neutrino interactions
455 
456  return;
457  }//end analyze
458 
460  {
461  }
462 
463 }//end namespace
464 
465 namespace ncs{
466 
468 
469 }
double E(const int i=0) const
Definition: MCParticle.h:232
back track the reconstruction to the simulation
int PdgCode() const
Definition: MCParticle.h:211
int CCNC() const
Definition: MCNeutrino.h:148
double QSqr() const
Definition: MCNeutrino.h:157
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
void reconfigure(const fhicl::ParameterSet &p)
double Py(const int i=0) const
Definition: MCParticle.h:230
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::string fMCFluxLabel
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
size_type size() const
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
int Mother() const
Definition: MCParticle.h:212
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
double Px(const int i=0) const
Definition: MCParticle.h:229
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Float_t ss
Definition: plot.C:24
int HitNuc() const
Definition: MCNeutrino.h:152
std::string fweightLabel
label for generator information
Loaders::FluxType flux
int NumberDaughters() const
Definition: MCParticle.h:216
object containing MC flux information
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
const key_type & TrackId(const size_type) const
int InteractionType() const
Definition: MCNeutrino.h:150
DEFINE_ART_MODULE(ROCKMRE)
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
double W() const
Definition: MCNeutrino.h:154
int fptype
Definition: MCFlux.h:63
void beginJob()
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
double Y() const
Definition: MCNeutrino.h:156
T get(std::string const &key) const
Definition: ParameterSet.h:231
int evt
void analyze(const art::Event &evt)
double energy
Definition: plottest35.C:25
caf::StandardRecord * sr
double X() const
Definition: MCNeutrino.h:155
const double j
Definition: BetheBloch.cxx:29
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:16
int printDecay
genie weights
This class describes a particle created in the detector Monte Carlo simulation.
OStream cout
Definition: OStream.cxx:6
double Vx(const int i=0) const
Definition: MCParticle.h:220
T * make(ARGS...args) const
Data resulting from a Hough transform on the cell hit positions.
int fntype
Definition: MCFlux.h:51
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Pz(const int i=0) const
Definition: MCParticle.h:231
double Vz(const int i=0) const
Definition: MCParticle.h:222
void endSubRun(const art::SubRun &sr)
Float_t e
Definition: plot.C:35
Event generator information.
Definition: MCNeutrino.h:18
Float_t w
Definition: plot.C:20
int Mode() const
Definition: MCNeutrino.h:149
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
mapped_type Particle(const size_type) const
enum BeamMode string