Xbeam_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 
21 #include "Simulation/Particle.h"
23 #include "SummaryData/POTSum.h"
24 #include "SummaryData/SpillData.h"
25 #include "RecoBase/CellHit.h"
26 #include "RecoBase/Cluster.h"
27 #include "RecoBase/RecoHit.h"
28 #include "RecoBase/Track.h"
29 #include "RecoBase/Vertex.h"
30 #include "RecoBase/Prong.h"
31 #include "RecoBase/HoughResult.h"
32 
33 #include "Utilities/AssociationUtil.h"
34 #include "MCCheater/BackTracker.h"
35 #include "Utilities/func/MathUtil.h"
36 
37 // Framework includes
46 #include "fhiclcpp/ParameterSet.h"
52 
53 #include <cmath>
54 #include <vector>
55 #include <string>
56 #include <vector>
57 
58 //include TMVA
59 #include "TMVA/Tools.h"
60 #include "TMVA/Reader.h"
61 #include "TMVA/MethodCuts.h"
62 
63 class TVector3;
64 class TH1F;
65 class TTree;
66 
67 namespace rb{class Cluster; class Track;}
68 namespace sim{class Particle;}
69 namespace simb{class MCFlux; class MCTruth; class MCNeutrino;}
70 
71 namespace ncs {
72  class Xbeam : public art::EDAnalyzer {
73  public:
74  explicit Xbeam(fhicl::ParameterSet const& pset);
75  virtual ~Xbeam();
76 
77  void beginJob();
78  void analyze(const art::Event& evt);
79  void reconfigure(const fhicl::ParameterSet& p);
80  void endSubRun(const art::SubRun& sr);
81  void endJob();
82  double fTotalPOT;
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 fSliceLabel; ///label for slices
92  std::string fMCFluxLabel; ///label for generator information
93  std::string fCellHitLabel; ///label for cell hits
94 
95  TProfile *slEff_Etrue;
96 
97  TTree *_otree;
98 
99  /*** true info ***/
100  float nu_Vx;//true vertex X
101  float nu_Vy;//true vertex Y
102  float nu_Vz;//true vertex Z
103  float nu_truePt;//true Pt
104  float nu_trueE;//true E
105  float nu_trueTheta;//true theta
106  int nu_ccnc;
107  int nu_PDG;
110  int nu_mode;
111  bool nu_isQE;
114  float nu_HadX;
115  float nu_HadY;
116  float nu_HadW;
117  float nu_qsqr;
118  float nu_OscWt;
120  float meson_Px;
121  float meson_Py;
122  float meson_Pz;
123  int lep_PDG;
124  float lep_Px;
125  float lep_Py;
126  float lep_Pz;
127  float lep_E;
128 
129  /*** slicer information ***/
130  int Nslices;
131  int hasSlice;
132 
133  };
134 }
135 
136 namespace ncs{
137 
138  //......................................................................
139  Xbeam::Xbeam(fhicl::ParameterSet const& pset)
140  : EDAnalyzer(pset)
141  {
142  mf::LogInfo("Xbeam")<<__PRETTY_FUNCTION__<<"\n";
143  reconfigure(pset);
144  }
145 
146  //......................................................................
148  {
149  }
150 
151  //......................................................................
153  {
154  fSliceLabel = pset.get< std::string >("SliceLabel");
155  fMCFluxLabel = pset.get< std::string >("MCFluxLabel");
156  fCellHitLabel = pset.get< std::string >("CellHitLabel");
157  }
158 
159  //......................................................................
161  {
163 
164  fTotalPOT=0.;
165  fTotalSpill=0;
166 
167  //histo
168  slEff_Etrue=tfs->make<TProfile>("slEff_Etrue","",20,0.,10.,0,1.2);
169 
170  //personal tuples
171  _otree = tfs->make<TTree>("Xbeam","Xbeam particle");
172 
173  _otree->Branch("nu_Vx",&nu_Vx,"nu_Vx/F");
174  _otree->Branch("nu_Vy",&nu_Vy,"nu_Vy/F");
175  _otree->Branch("nu_Vz",&nu_Vz,"nu_Vz/F");
176  _otree->Branch("nu_truePt",&nu_truePt,"nu_truePt/F");
177  _otree->Branch("nu_trueE",&nu_trueE,"nu_trueE/F");
178  _otree->Branch("nu_trueTheta",&nu_trueTheta,"nu_trueTheta/F");
179  _otree->Branch("nu_ccnc",&nu_ccnc,"nu_ccnc/I");
180  _otree->Branch("nu_PDG",&nu_PDG,"nu_PDG/I");
181  _otree->Branch("nu_orig_parentPDG",&nu_orig_parentPDG,"nu_orig_parentPDG/I");
182  _otree->Branch("nu_origPDG",&nu_origPDG,"nu_origPDG/I");
183  _otree->Branch("nu_mode",&nu_mode,"nu_mode/I");
184  _otree->Branch("nu_isQE",&nu_isQE,"nu_isQE/B");
185  _otree->Branch("nu_intType",&nu_intType,"nu_intType/I");
186  _otree->Branch("nu_hitNucl",&nu_hitNucl,"nu_hitNucl/I");
187  _otree->Branch("nu_HadX",&nu_HadX,"nu_HadX/F");
188  _otree->Branch("nu_HadY",&nu_HadY,"nu_HadY/F");
189  _otree->Branch("nu_HadW",&nu_HadW,"nu_HadW/F");
190  _otree->Branch("nu_qsqr",&nu_qsqr,"nu_qsqr/F");
191  _otree->Branch("nu_OscWt",&nu_OscWt,"nu_OscWt/F");
192  _otree->Branch("meson_PDG",&meson_PDG,"meson_PDG/I");
193  _otree->Branch("meson_Px",&meson_Px,"meson_Px/F");
194  _otree->Branch("meson_Py",&meson_Py,"meson_Py/F");
195  _otree->Branch("meson_Pz",&meson_Pz,"meson_Pz/F");
196  _otree->Branch("lep_PDG",&lep_PDG,"lep_PDG/I");
197  _otree->Branch("lep_Px",&lep_Px,"lep_Px/F");
198  _otree->Branch("lep_Py",&lep_Py,"lep_Py/F");
199  _otree->Branch("lep_Pz",&lep_Pz,"lep_Pz/F");
200  _otree->Branch("lep_E",&lep_E,"lep_E/F");
201 
202  _otree->Branch("hasSlice",&hasSlice,"hasSlice/I");
203  _otree->Branch("Nslices",&Nslices,"Nslices/I");
204 
205  }
206 
207  //......................................................................
209  {
211  sr.getByLabel(fMCFluxLabel, pot);
212 
213 
214  if( !pot.failedToGet() ){
215  fTotalPOT += pot->totgoodpot;
216  fTotalSpill += pot->goodspills;
217 
218  std::cout << "POT in subrun " << sr.subRun()
219  << " = " << pot->totgoodpot <<" with "<<fTotalSpill<<" spills"<< std::endl;
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  nu_Vx = -999.;
237  nu_Vy = -999.;
238  nu_Vz = -999.;
239  nu_truePt = -1.;
240  nu_trueE = -1.;
241  nu_trueTheta = -1.;
242 
243  nu_ccnc=-1;
244  nu_PDG=-1;
245  nu_origPDG=-1;
247  nu_mode=-1;
248  nu_isQE=false;
249  nu_intType=-1;
250  nu_hitNucl=-1;
251 
252  nu_HadX=-1.;
253  nu_HadY=-1.;
254  nu_HadW=-1.;
255  nu_qsqr=-1.;
256  nu_OscWt=-1.;
257 
258  meson_PDG=-1;
259  meson_Px=-1.;
260  meson_Py=-1.;
261  meson_Pz=-1.;
262 
263  lep_PDG=-1;
264  lep_Px=-1.;
265  lep_Py=-1.;
266  lep_Pz=-1.;
267  lep_E=-1.;
268 
269  Nslices=-1;
270  hasSlice=-1;
271 
274  //////////////
275  //Truth info//
276  //////////////
278  evt.getByLabel(fMCFluxLabel,mclist);
280  evt.getByLabel(fMCFluxLabel, fllist);
281  //# of interactions
282  //std::cout<<"# of interactions: "<<mclist->size()<<std::endl;
283  if( mclist->size()!=1 ){
284  std::cout<<"ABNORMAL! "<<mclist->size()<<" interactions."<<std::endl;
285  return;
286  }
287 
288  art::Ptr<simb::MCTruth> mctruth(mclist,0);
289  const simb::MCNeutrino& nu = mctruth->GetNeutrino();
290 
291  nu_trueE = nu.Nu().E();
292  nu_PDG = nu.Nu().PdgCode();
293  nu_Vx=nu.Nu().Vx();
294  nu_Vy=nu.Nu().Vy();
295  nu_Vz=nu.Nu().Vz();
296 
297  nu_truePt = nu.Pt();
298  nu_trueTheta = nu.Theta();
299  nu_ccnc = nu.CCNC();
300  nu_mode = nu.Mode();
301  nu_isQE = nu.InteractionType() == simb::kQE ||
304  nu_hitNucl = nu.HitNuc();
305  nu_HadX = nu.X();
306  nu_HadY = nu.Y();
307  nu_HadW = nu.W();
308  nu_qsqr = nu.QSqr();
309 
310  lep_PDG = nu.Lepton().PdgCode();
311  lep_Px = nu.Lepton().Px();
312  lep_Py = nu.Lepton().Py();
313  lep_Pz = nu.Lepton().Pz();
314  lep_E = nu.Lepton().E();
315 
316 
317  art::Ptr<simb::MCFlux> flux(fllist,0);
318 
319  nu_origPDG = flux->fntype;
320  nu_orig_parentPDG = flux->fptype;
321  meson_PDG = flux->ftptype;
322  meson_Px = flux->ftpx;
323  meson_Py = flux->ftpy;
324  meson_Pz = flux->ftpz;
325 
326  //std::cout<<i_intx<<" "<<nu.CCNC()<<" "<<flux->fntype<<" "<<nu.Nu().PdgCode()<<std::endl;
327  //std::cout<<flux->fntype<<" "<<flux->fptype<<" "<<flux->ftptype<<std::endl;
328  //std::cout<<"meson is "<<flux->ftptype<<std::endl;
329  //CC nue only
330  if( nu.CCNC() != 0 ) return;
331  if( nu.Nu().PdgCode() != 12 ) return;
332  if( flux->fntype != 12 ) return;
333 
334  ////////////////////////
335  //get reco slicer info//
336  ////////////////////////
338  evt.getByLabel(fSliceLabel, slices);
339  //get all hits info
341  evt.getByLabel(fCellHitLabel, chits);
342 
343  bool matchSlice=false;
344 
345  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
346 
347  const rb::Cluster& slice = (*slices)[sliceIdx];
348  if(slice.IsNoise()) continue;
351  for(size_t h = 0; h < chits->size(); ++h)
352  allhits.push_back(art::Ptr<rb::CellHit>(chits, h));
353  const art::Ptr<rb::Cluster> allslice(slices,sliceIdx);
354  std::vector<cheat::NeutrinoEffPur> iacts = bt->SliceToNeutrinoInteractions(allslice, allhits);
355 
356  if( iacts.empty() ){
357  std::cout<<"no matching between slice and nu interactions!"<<std::endl;
358  return;
359  }
360  art::Ptr<simb::MCTruth>& truth = iacts[0].neutrinoInt;
362  pv.push_back(truth);
364  std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
366  if(fluxs.size() == 1){
367  flux = *fluxs[0];
368  }
369  else{
370  std::cout<<"There is no assoicated neutrino intereaction !"<<std::endl;
371  }
372  const simb::MCNeutrino& nu = truth->GetNeutrino();
373  if( flux.fntype==12 && nu.Nu().PdgCode()==12 && nu.CCNC()==0 ) matchSlice=true;
374  }
375  if( matchSlice ) hasSlice=1;
376  else hasSlice=0;
377 
378  slEff_Etrue->Fill(nu.Nu().E(),matchSlice?1:0);
379 
380  _otree->Fill();
381 
382  return;
383  }//end analyze
384 
386  {
387  }
388 
389 }//end namespace
390 
391 namespace ncs{
392 
394 
395 }
396 
double E(const int i=0) const
Definition: MCParticle.h:232
std::string fCellHitLabel
label for generator information
Definition: Xbeam_module.cc:93
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
double Theta() const
angle between incoming and outgoing leptons, in radians
Definition: MCNeutrino.cxx:63
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
double Py(const int i=0) const
Definition: MCParticle.h:230
float nu_trueE
double ftpx
Definition: MCFlux.h:79
void analyze(const art::Event &evt)
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
SubRunNumber_t subRun() const
Definition: SubRun.h:44
std::vector< NeutrinoEffPur > SliceToNeutrinoInteractions(const std::vector< const rb::CellHit * > &sliceHits, const std::vector< const rb::CellHit * > &allHits, bool sortPur=false) const
Given a collection of hits (often a slice), returns vector of structures of neutrino interactions...
void endJob()
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
int ftptype
Definition: MCFlux.h:82
const char * p
Definition: xmltok.h:285
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
double Pt() const
transverse momentum of interaction, in GeV/c
Definition: MCNeutrino.cxx:74
double Px(const int i=0) const
Definition: MCParticle.h:229
std::string fMCFluxLabel
label for slices
Definition: Xbeam_module.cc:92
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
int HitNuc() const
Definition: MCNeutrino.h:152
A collection of associated CellHits.
Definition: Cluster.h:47
std::string fSliceLabel
Definition: Xbeam_module.cc:91
charged current quasi-elastic
Definition: MCNeutrino.h:96
int nu_orig_parentPDG
float nu_trueTheta
Loaders::FluxType flux
float meson_Py
object containing MC flux information
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
void beginJob()
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
#define pot
caf::StandardRecord * sr
double X() const
Definition: MCNeutrino.h:155
double ftpz
Definition: MCFlux.h:81
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
float nu_OscWt
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
This class describes a particle created in the detector Monte Carlo simulation.
virtual ~Xbeam()
float meson_Pz
OStream cout
Definition: OStream.cxx:6
double Vx(const int i=0) const
Definition: MCParticle.h:220
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
TTree * _otree
Definition: Xbeam_module.cc:97
TProfile * slEff_Etrue
label for cell hits
Definition: Xbeam_module.cc:95
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
int goodspills
Definition: POTSum.h:31
float meson_Px
void geom(int which=0)
Definition: geom.C:163
double Pz(const int i=0) const
Definition: MCParticle.h:231
float nu_truePt
double fTotalPOT
Definition: Xbeam_module.cc:82
double Vz(const int i=0) const
Definition: MCParticle.h:222
double ftpy
Definition: MCFlux.h:80
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
Event generator information.
Definition: MCNeutrino.h:18
Float_t w
Definition: plot.C:20
void endSubRun(const art::SubRun &sr)
int Mode() const
Definition: MCNeutrino.h:149
double Vy(const int i=0) const
Definition: MCParticle.h:221
Encapsulate the geometry of one entire detector (near, far, ndos)
bool failedToGet() const
Definition: Handle.h:196
void reconfigure(const fhicl::ParameterSet &p)