Xeff_module.cc
Go to the documentation of this file.
1 ///
2 /// \brief
3 /// \author Xuebing Bu - xbbu@fnal.gov
4 ///
5 #include "ShowerLID/EventLID.h"
6 #include "ShowerLID/ShowerLID.h"
7 #include "Calibrator/Calibrator.h"
8 
9 #include "Geometry/Geometry.h"
10 #include "GeometryObjects/CellGeo.h"
11 
12 // ROOT includes
13 #include "TVector3.h"
14 #include "TH1F.h"
15 #include "TMath.h"
16 #include "TTree.h"
17 #include "TProfile2D.h"
18 #include "time.h"
19 
20 #include "nusimdata/SimulationBase/MCFlux.h"
21 #include "nusimdata/SimulationBase/MCTruth.h"
22 #include "Simulation/Particle.h"
23 #include "Simulation/ParticleNavigator.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
39 #include "art/Framework/Core/EDAnalyzer.h"
40 #include "art/Framework/Core/ModuleMacros.h"
41 #include "canvas/Persistency/Common/PtrVector.h"
42 #include "canvas/Persistency/Common/Assns.h"
43 #include "canvas/Persistency/Common/FindManyP.h"
44 #include "canvas/Persistency/Common/FindMany.h"
45 #include "art/Framework/Principal/Event.h"
46 #include "art/Framework/Principal/SubRun.h"
47 #include "fhiclcpp/ParameterSet.h"
48 #include "art/Framework/Principal/Handle.h"
49 #include "art/Framework/Services/Registry/ServiceHandle.h"
50 #include "art/Framework/Services/Optional/TFileService.h"
51 #include "art/Framework/Services/Optional/TFileDirectory.h"
52 #include "messagefacility/MessageLogger/MessageLogger.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 Xeff : public art::EDAnalyzer {
74  public:
75  explicit Xeff(fhicl::ParameterSet const& pset);
76  virtual ~Xeff();
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  double fTotalPOT;
84  int fTotalSpill;
85  TTree *_btree;
86 
87  private:
88 
89  const rb::RecoHit MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
90  double const& w);
91 
92  //folders to get objects from
93  std::string fSliceLabel; ///label for slices
94  std::string fKTrackLabel; ///label for kalman tracks
95  std::string fProngLabel; ///label for prong
96  std::string fInstLabel3D; ///instance label for prongs 3D
97  std::string fInstLabel2D; ///instance label for prongs 2D
98  std::string fMCFluxLabel; ///label for generator information
99  std::string fVertexLabel; ///label for vertex
100  std::string fCellHitLabel; ///label for cell hits
101  std::string flidLabel;
102 
103  //save event info to ntuples
104  int runnum;
105  int subrunnum;
106  int evtnum;
107  int Nslices;
108  double spillPOT;
109  int Ndeaddcm;
110 
111  TTree *_otree;
112  //TTree *_evttree;
113  int sl_runnum;
114  int sl_subrunnum;
115  int sl_evtnum;
116  int sl_Nslices;
117  double sl_spillPOT;
118  int sl_Ndeaddcm;
119 
120  /*** true info ***/
121  float truePt;//true Pt
122  float trueE;//true E
123  float trueTheta;//true theta
124  int ccnc;
125  int PDG;
126  int origPDG;
127  int mode;
128  bool isQE;
129  int intType;
130  int hitNucl;
131  float HadX;
132  float HadY;
133  float HadW;
134  float qsqr;
135  float OscWt;
136  int lepPDG;
137  float lepPx;
138  float lepPy;
139  float lepPz;
140  float lepE;
141  int mesonPDG;
142  float mesonPx;
143  float mesonPy;
144  float mesonPz;
145  float Vx;//true vertex X
146  float Vy;//true vertex Y
147  float Vz;//true vertex Z
148 
149  /*** Reco Vertex ***/
150  int hasVertex;
151  int hasOneVertex;
152  float recoVx;//reco vertex X
153  float recoVy;//reco vertex Y
154  float recoVz;//reco vertex Z
155 
156  /*** slicer information ***/
157  int Ncells;//# of associate cells
158  int Nplanes;//# of extended planes
159  float Eslice;//reco energy
160  int sliceIndex;
161 
162  /*** longest kalman track ***/
163  float KtrackLength;//length of longest track
164  int KtrackNcells;//# of cells associated to longest track
165  float KtrackEnergy;//energy associated to longest track
166  float KtrackNcells_ratio;//# of cells associated to longest track over Ncells
167  float KtrackEnergy_ratio;//energy associated to longest track over recoE
168  float KtrackNplanes;
169  float KtrackTheta;
170  float KtrackPhi;
171  float KtrackStartX;
172  float KtrackStartY;
173  float KtrackStartZ;
174  float KtrackStopX;
175  float KtrackStopY;
176  float KtrackStopZ;
177  int NKtracks;//# of tracks
178  int Ktrack3D;
179 
180  /*** mip information ***/
181  float Nmip;//# of mip cells
182  float Fmip;//# of mip cells over Ncells
183 
184  /*** prong information ***/
185  int Nprongs;
186  float Ebalance;
187  float Dphi;
188 
189  float p1Energy;
190  float p1Phi;
191  float p1Theta;
192  float p1StartX;
193  float p1StartY;
194  float p1StartZ;
195  float p1Length;
196  float p1StopX;
197  float p1StopY;
198  float p1StopZ;
199  float p1EnergyXView;//energy for XView cells
200  float p1EnergyYView;//energy for YView cells
201  int p1NCells;//# of cells
202  int p1NCellsXView;//# of cells for XView
203  int p1NCellsYView;//# of cells for YView
204  int p1Nplanes;
205  float p1Fmip;
206  //bdt inputs
207  float p2Fmip;
208  float emaxfrac_6p;
209  float efrac_10p;
210  float efrac_p2;
211  float efrac_p3;
212  float efrac_p4;
213  float efrac_2sig;
214 
215  //leading prong
216  int prong1cellNcell;
217  float prong1cellX[500];
218  float prong1cellY[500];
219  float prong1cellZ[500];
220  float prong1cellE[500];
221  int prong1cellADC[500];
222  int prong1cellPlane[500];
223  int prong1cellCell[500];
224  float prong1cellPECorr[500];
225  float prong1cellW[500];
226 
227  float lid;
228  };
229 }
230 
231 namespace ncs{
232 
233 /*
234  //BDT reader
235  TMVA::Reader *_bdt=new TMVA::Reader();
236  //input variable values
237  float inputvars[15]={0.};
238  //input variable names
239  TString inputs[15]={"_ncells","_energy","_ltrack","_ltrack_ncells_ratio","_fmip",
240  "_nmip","_efrac_20planes","_efrac_slide2","_efrac_slide6",
241  "_efrac_2sigma","_eiso_3sigma","_nprongs","_ebalance",
242  "_nprongs2D","_ebalance2D"};
243 */
244  //......................................................................
245  Xeff::Xeff(fhicl::ParameterSet const& pset)
246  : EDAnalyzer(pset)
247  {
248  mf::LogInfo("Xeff")<<__PRETTY_FUNCTION__<<"\n";
249  reconfigure(pset);
250  }
251 
252  //......................................................................
253  Xeff::~Xeff()
254  {
255  }
256 
257  //......................................................................
258  void Xeff::reconfigure(const fhicl::ParameterSet& pset)
259  {
260  fSliceLabel = pset.get< std::string >("SliceLabel");
261  fKTrackLabel = pset.get< std::string >("KTrackLabel");
262  fProngLabel = pset.get< std::string >("ProngLabel");
263  fInstLabel3D = pset.get< std::string >("InstLabel3D");
264  fInstLabel2D = pset.get< std::string >("InstLabel2D");
265  fMCFluxLabel = pset.get< std::string >("MCFluxLabel");
266  fCellHitLabel = pset.get< std::string >("CellHitLabel");
267  fVertexLabel = pset.get< std::string >("VertexLabel");
268 
269  flidLabel = pset.get< std::string >("lidLabel");
270  }
271 
272  //......................................................................
273  void Xeff::beginJob()
274  {
275  art::ServiceHandle<art::TFileService> tfs;
276 
277  fTotalPOT=0.;
278  fTotalSpill=0;
279 
280 /*
281  for( int i=0; i<15; ++i ){
282  _bdt->AddVariable(inputs[i],&inputvars[i]);
283  }
284 
285  _bdt->BookMVA("BDT method","nonoscBDTweights_20140304.xml");
286 */
287 
288  //histo
289  //Ereco=tfs->make<TH1F>("Ereco","Ereco;Ereco;# of events",100,0.,5.);
290  //personal tuples
291  _otree = tfs->make<TTree>("Xeff","Xeff particle");
292 
293  _otree->Branch("sl_runnum",&sl_runnum,"sl_runnum/I");
294  _otree->Branch("sl_subrunnum",&sl_subrunnum,"sl_subrunnum/I");
295  _otree->Branch("sl_evtnum",&sl_evtnum,"sl_evtnum/I");
296  _otree->Branch("sl_Nslices",&sl_Nslices,"sl_Nslices/I");
297  _otree->Branch("sl_spillPOT",&sl_spillPOT,"sl_spillPOT/D");
298  _otree->Branch("sl_Ndeaddcm",&sl_Ndeaddcm,"sl_Ndeaddcm/I");
299 
300  _otree->Branch("truePt",&truePt,"truePt/F");
301  _otree->Branch("trueE",&trueE,"trueE/F");
302  _otree->Branch("trueTheta",&trueTheta,"trueTheta/F");
303  _otree->Branch("ccnc",&ccnc,"ccnc/I");
304  _otree->Branch("PDG",&PDG,"PDG/I");
305  _otree->Branch("origPDG",&origPDG,"origPDG/I");
306  _otree->Branch("mode",&mode,"mode/I");
307  _otree->Branch("isQE",&isQE,"isQE/B");
308  _otree->Branch("intType",&intType,"intType/I");
309  _otree->Branch("hitNucl",&hitNucl,"hitNucl/I");
310  _otree->Branch("HadX",&HadX,"HadX/F");
311  _otree->Branch("HadY",&HadY,"HadY/F");
312  _otree->Branch("HadW",&HadW,"HadW/F");
313  _otree->Branch("qsqr",&qsqr,"qsqr/F");
314  _otree->Branch("OscWt",&OscWt,"OscWt/F");
315  _otree->Branch("lepE",&lepE,"lepE/F");
316  _otree->Branch("lepPx",&lepPx,"lepPx/F");
317  _otree->Branch("lepPy",&lepPy,"lepPy/F");
318  _otree->Branch("lepPz",&lepPz,"lepPz/F");
319  _otree->Branch("lepPDG",&lepPDG,"lepPDG/I");
320  _otree->Branch("mesonPDG",&mesonPDG,"mesonPDG/I");
321  _otree->Branch("mesonPx",&mesonPx,"mesonPx/F");
322  _otree->Branch("mesonPy",&mesonPy,"mesonPy/F");
323  _otree->Branch("mesonPz",&mesonPz,"mesonPz/F");
324  _otree->Branch("Vx",&Vx,"Vx/F");
325  _otree->Branch("Vy",&Vy,"Vy/F");
326  _otree->Branch("Vz",&Vz,"Vz/F");
327  _otree->Branch("hasVertex",&hasVertex,"hasVertex/I");
328  _otree->Branch("hasOneVertex",&hasOneVertex,"hasOneVertex/I");
329  _otree->Branch("recoVx",&recoVx,"recoVx/F");
330  _otree->Branch("recoVy",&recoVy,"recoVy/F");
331  _otree->Branch("recoVz",&recoVz,"recoVz/F");
332 
333  _otree->Branch("NKtracks",&NKtracks,"NKtracks/I");
334  _otree->Branch("KtrackLength",&KtrackLength,"KtrackLength/F");
335  _otree->Branch("KtrackNcells",&KtrackNcells,"KtrackNcells/I");
336  _otree->Branch("KtrackEnergy",&KtrackEnergy,"KtrackEnergy/F");
337  _otree->Branch("KtrackNcells_ratio",&KtrackNcells_ratio,"KtrackNcells_ratio/F");
338  _otree->Branch("KtrackEnergy_ratio",&KtrackEnergy_ratio,"KtrackEnergy_ratio/F");
339  _otree->Branch("KtrackNplanes",&KtrackNplanes,"KtrackNplanes/F");
340  _otree->Branch("KtrackPhi",&KtrackPhi,"KtrackPhi/F");
341  _otree->Branch("KtrackTheta",&KtrackTheta,"KtrackTheta/F");
342  _otree->Branch("KtrackStartX",&KtrackStartX,"KtrackStartX/F");
343  _otree->Branch("KtrackStartY",&KtrackStartY,"KtrackStartY/F");
344  _otree->Branch("KtrackStartZ",&KtrackStartZ,"KtrackStartZ/F");
345  _otree->Branch("KtrackStopX",&KtrackStopX,"KtrackStopX/F");
346  _otree->Branch("KtrackStopY",&KtrackStopY,"KtrackStopY/F");
347  _otree->Branch("KtrackStopZ",&KtrackStopZ,"KtrackStopZ/F");
348  _otree->Branch("Ktrack3D",&Ktrack3D,"Ktrack3D/I");
349 
350  _otree->Branch("Ncells",&Ncells,"Ncells/I");
351  _otree->Branch("Nplanes",&Nplanes,"Nplanes/I");
352  _otree->Branch("Eslice",&Eslice,"Eslice/F");
353  _otree->Branch("sliceIndex",&sliceIndex,"sliceIndex/I");
354 
355  _otree->Branch("Nmip",&Nmip,"Nmip/F");
356  _otree->Branch("Fmip",&Fmip,"Fmip/F");
357 
358  _otree->Branch("Nprongs",&Nprongs,"Nprongs/I");
359  _otree->Branch("Ebalance",&Ebalance,"Ebalance/F");
360  _otree->Branch("Dphi",&Dphi,"Dphi/F");
361 
362  _otree->Branch("p1Length",&p1Length,"p1Length/F");
363  _otree->Branch("p1StopX",&p1StopX,"p1StopX/F");
364  _otree->Branch("p1StopY",&p1StopY,"p1StopY/F");
365  _otree->Branch("p1StopZ",&p1StopZ,"p1StopZ/F");
366  _otree->Branch("p1Energy",&p1Energy,"p1Energy/F");
367  _otree->Branch("p1Phi",&p1Phi,"p1Phi/F");
368  _otree->Branch("p1Theta",&p1Theta,"p1Theta/F");
369  _otree->Branch("p1StartX",&p1StartX,"p1StartX/F");
370  _otree->Branch("p1StartY",&p1StartY,"p1StartY/F");
371  _otree->Branch("p1StartZ",&p1StartZ,"p1StartZ/F");
372  _otree->Branch("p1EnergyXView",&p1EnergyXView,"p1EnergyXView/F");
373  _otree->Branch("p1EnergyYView",&p1EnergyYView,"p1EnergyYView/F");
374  _otree->Branch("p1NCells",&p1NCells,"p1NCells/I");
375  _otree->Branch("p1NCellsXView",&p1NCellsXView,"p1NCellsXView/I");
376  _otree->Branch("p1NCellsYView",&p1NCellsYView,"p1NCellsYView/I");
377  _otree->Branch("p1Nplanes",&p1Nplanes,"p1Nplanes/I");
378  _otree->Branch("p1Fmip",&p1Fmip,"p1Fmip/F");
379 
380  _otree->Branch("p2Fmip",&p2Fmip,"p2Fmip/F");
381  _otree->Branch("emaxfrac_6p",&emaxfrac_6p,"emaxfrac_6p/F");
382  _otree->Branch("efrac_10p",&efrac_10p,"efrac_10p/F");
383  _otree->Branch("efrac_p2",&efrac_p2,"efrac_p2/F");
384  _otree->Branch("efrac_p3",&efrac_p3,"efrac_p3/F");
385  _otree->Branch("efrac_p4",&efrac_p4,"efrac_p4/F");
386  _otree->Branch("efrac_2sig",&efrac_2sig,"efrac_2sig/F");
387 
388  _otree->Branch("prong1cellNcell",&prong1cellNcell,"prong1cellNcell/I");
389  _otree->Branch("prong1cellX",&prong1cellX,"prong1cellX[prong1cellNcell]/F");
390  _otree->Branch("prong1cellY",&prong1cellY,"prong1cellY[prong1cellNcell]/F");
391  _otree->Branch("prong1cellZ",&prong1cellZ,"prong1cellZ[prong1cellNcell]/F");
392  _otree->Branch("prong1cellE",&prong1cellE,"prong1cellE[prong1cellNcell]/F");
393  _otree->Branch("prong1cellADC",&prong1cellADC,"prong1cellADC[prong1cellNcell]/I");
394  _otree->Branch("prong1cellPlane",&prong1cellPlane,"prong1cellPlane[prong1cellNcell]/I");
395  _otree->Branch("prong1cellCell",&prong1cellCell,"prong1cellCell[prong1cellNcell]/I");
396  _otree->Branch("prong1cellW",&prong1cellW,"prong1cellW[prong1cellNcell]/F");
397  _otree->Branch("prong1cellPECorr",&prong1cellPECorr,"prong1cellPECorr[prong1cellNcell]/F");
398 
399  _otree->Branch("lid",&lid,"lid/F");
400 
401  _btree=new TTree("beamInfo","beam info");
402  _btree->Branch("runnum",&runnum,"runnum/I");
403  _btree->Branch("subrunnum",&subrunnum,"subrunnum/I");
404  _btree->Branch("evtnum",&evtnum,"evtnum/I");
405  _btree->Branch("Nslices",&Nslices,"Nslices/I");
406  _btree->Branch("spillPOT",&spillPOT,"spillPOT/D");
407  _btree->Branch("Ndeaddcm",&Ndeaddcm,"Ndeaddcm/I");
408 
409  }
410 
411  //......................................................................
412  void Xeff::endSubRun(const art::SubRun& sr)
413  {
414  art::Handle<sumdata::POTSum> pot;
415  sr.getByLabel(fMCFluxLabel, pot);
416 
417 
418  if( !pot.failedToGet() ){
419  fTotalPOT += pot->totgoodpot;
420  fTotalSpill += pot->goodspills;
421 
422  std::cout << "POT in subrun " << sr.subRun()
423  << " = " << pot->totgoodpot <<" with "<<fTotalSpill<<" spills"<< std::endl;
424  }
425  }
426 
427  //......................................................................
428  const rb::RecoHit Xeff::MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
429  double const& w)
430  {
431  art::ServiceHandle<calib::Calibrator> cal;
432  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
433 
434  return rhit;
435  }
436 
437  //......................................................................
438  void Xeff::analyze(const art::Event& evt)
439  {
440  sl_runnum=-1;
441  sl_subrunnum=-1;
442  sl_evtnum=-1;
443  sl_Nslices=-1;
444  sl_spillPOT=-1.;
445  sl_Ndeaddcm=-1;
446  Ndeaddcm=-1;
447 
448  runnum=-1;
449  subrunnum=-1;
450  evtnum=-1;
451  Nslices=-1;
452  spillPOT=0.;
453 
454  lid=-99.;
455 
456  Vx=-999.;
457  Vy=-999.;
458  Vz=-999.;
459  hasVertex=-1;
460  hasOneVertex=-1;
461  recoVx=-999.;
462  recoVy=-999.;
463  recoVz=-999.;
464 
465  truePt=-1.;
466  trueE=-1.;
467  ccnc=-1;
468  PDG=-1;
469  origPDG=-1;
470  mode=-1;
471  isQE=false;
472  intType=-1;
473  hitNucl=-1;
474  HadX=-999.;
475  HadY=-999.;
476  HadW=-999.;
477  qsqr=-999.;
478  OscWt=-999.;
479 
480  lepPDG=-1;
481  lepPx=-1.;
482  lepPy=-1.;
483  lepPz=-1.;
484  mesonPDG=-1;
485  mesonPx=-1.;
486  mesonPy=-1.;
487  mesonPz=-1.;
488 
489  Ncells=0;
490  Nplanes=0;
491  Eslice=-1.;
492  sliceIndex=-1;
493 
494  KtrackLength=0.;
495  KtrackNcells=0;
496  KtrackEnergy=0.;
497  KtrackEnergy_ratio=0.;
498  KtrackNcells_ratio=0.;
499  KtrackNplanes=0;
500  KtrackTheta=-1.;
501  KtrackPhi=-1.;
502  KtrackStartX=-999.;
503  KtrackStartY=-999.;
504  KtrackStartZ=-999.;
505  KtrackStopX=-999.;
506  KtrackStopY=-999.;
507  KtrackStopZ=-999.;
508  NKtracks=0.;
509 
510  Nmip=-1.;
511  Fmip=-1.;
512 
513  Nprongs=0;
514  Ebalance=0.;
515  Dphi=-1.;
516 
517  p1Energy=0.;
518  p1Phi=-1.;
519  p1Theta=-1.;
520  p1StartX=-999.;
521  p1StartY=-999.;
522  p1StartZ=-999.;
523  p1StopX=-999.;
524  p1StopY=-999.;
525  p1StopZ=-999.;
526  p1Length=-1.;
527  p1EnergyXView=0.;
528  p1EnergyYView=0.;
529  p1NCells=0;
530  p1NCellsXView=0;
531  p1NCellsYView=0;
532  p1Nplanes=0;
533  p1Fmip=-1.;
534  p2Fmip=0.;
535  emaxfrac_6p=0.;
536  efrac_10p=0.;
537  efrac_p2=0.;
538  efrac_p3=0.;
539  efrac_p4=0.;
540  efrac_2sig=0.;
541 
542  prong1cellNcell=0;
543  for( int i=0; i<500; ++i ){
544  prong1cellX[i]=-999.;
545  prong1cellY[i]=-999.;
546  prong1cellZ[i]=-999.;
547  prong1cellE[i]=-999.;
548  prong1cellADC[i]=0;
549  prong1cellPlane[i]=-1;
550  prong1cellCell[i]=-1;
551  prong1cellW[i]=-1.;
552  prong1cellPECorr[i]=-1.;
553 
554  }
555 
556  art::ServiceHandle<calib::Calibrator> cal;
557  art::ServiceHandle<geo::Geometry> geom;
558 
559  //event info
560  runnum = evt.run();
561  subrunnum = evt.subRun();
562  evtnum = evt.id().event();
563 
564  //get all hits info
565  art::Handle<std::vector<rb::CellHit> > chits;
566  evt.getByLabel(fCellHitLabel, chits);
567 
568  std::vector<int> sq_plane;
569  sq_plane.clear();
570  std::vector<int> sq_xview;
571  sq_xview.clear();
572  std::vector<float> sq_tns;
573  sq_tns.clear();
574  std::vector<int> sq_cell;
575  sq_cell.clear();
576 
577  int sq_ncells=chits->size();
578  for(int i = 0; i < sq_ncells; ++i){
579  art::Ptr< rb::CellHit> hit(chits, i);
580  double hittime=hit->TNS()/1000.;
581  sq_plane.push_back(hit->Plane());
582  sq_tns.push_back(hittime);
583  sq_cell.push_back(hit->Cell());
584  if( (hit->View())==geo::kX ) sq_xview.push_back(1);
585  else sq_xview.push_back(0);
586  }//loop all hits
587 
588  int Nmissingdcm[14]={0};
589  for( int ic=0; ic<sq_ncells; ++ic ){
590  if( sq_plane[ic]<64 ){
591  if( sq_xview[ic]==1 ){
592  if( sq_cell[ic]>63 ) Nmissingdcm[0]++;
593  else Nmissingdcm[1]++;
594  }//XZ view
595  else {
596  if( sq_cell[ic]>31 ) Nmissingdcm[2]++;
597  else Nmissingdcm[3]++;
598  }//YZ view
599  }//DB1
600  else if( sq_plane[ic]<128 ){
601  if( sq_xview[ic]==1 ){
602  if( sq_cell[ic]>63 ) Nmissingdcm[4]++;
603  else Nmissingdcm[5]++;
604  }//XZ view
605  else {
606  if( sq_cell[ic]>31 ) Nmissingdcm[6]++;
607  else Nmissingdcm[7]++;
608  }//YZ view
609  }//DB2
610  else if( sq_plane[ic]<192 ){
611  if( sq_xview[ic]==1 ){
612  if( sq_cell[ic]>63 ) Nmissingdcm[8]++;
613  else Nmissingdcm[9]++;
614  }//XZ view
615  else{
616  if( sq_cell[ic]>31 ) Nmissingdcm[10]++;
617  else Nmissingdcm[11]++;
618  }//YZ view
619  }//DB3
620  else{
621  if( sq_xview[ic]==1 ) Nmissingdcm[12]++;
622  else Nmissingdcm[13]++;
623  }//muon-catcher
624  }//loop all cells
625 
626  int ndeaddcms=0;
627  for( int id=0; id<14; ++id ){
628  if( Nmissingdcm[id]==0 ) ndeaddcms++;
629  }
630 
631  Ndeaddcm=ndeaddcms;
632 
633  ////////////////////////
634  //get reco slicer info//
635  ////////////////////////
636  art::Handle< std::vector<rb::Cluster> > slices;
637  evt.getByLabel(fSliceLabel, slices);
638 
639  //get tracks, verticies and showers associated with the slices
640  //art::FindManyP<rb::Track> fmDTrack(slices, evt, fDTrackLabel);
641  art::FindManyP<rb::Track> fmKTrack(slices, evt, fKTrackLabel);
642  art::FindManyP<rb::Vertex> fmVertex(slices, evt, fVertexLabel);
643  //art::FindManyP<jmshower::JMShower> fmShower(slices, evt, fShowerLabel);
644 
645  //std::cout<<evt.run()<<" "<<evt.subRun()<<" "<<evt.id().event()<<std::endl;
646  //std::cout<<"# of slices: "<<slices->size()<<std::endl;
647 
648  int nslice=0;
649  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
650 
651  const rb::Cluster& slice = (*slices)[sliceIdx];
652  if(slice.IsNoise()) continue;
653  ++nslice;
654 
655  }//loop slices
656 
657  Nslices=nslice;
658  //std::cout<<"# of slice: "<<nslice<<std::endl;
659 
660  _btree->Fill();
661 
662  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
663 
664  const rb::Cluster& slice = (*slices)[sliceIdx];
665  if(slice.IsNoise()) continue;
666 
667  //get MC truth info
668  art::ServiceHandle<cheat::BackTracker> bt;
669  art::PtrVector<rb::CellHit> allhits;
670  for(size_t h = 0; h < chits->size(); ++h)
671  allhits.push_back(art::Ptr<rb::CellHit>(chits, h));
672  const art::Ptr<rb::Cluster> allslice(slices,sliceIdx);
673  std::vector<cheat::NeutrinoEffPur> iacts = bt->SliceToNeutrinoInteractions(allslice, allhits);
674 
675  if( iacts.empty() ){
676  std::cout<<"no matching between slice and nu interactions!"<<std::endl;
677  return;
678  }
679  // Pick the best-matching one
680  art::Ptr<simb::MCTruth>& truth = iacts[0].neutrinoInt;
681  // Assn interface is needlessly confusing. Have to pack stuff up like this
682  art::PtrVector<simb::MCTruth> pv;
683  pv.push_back(truth);
684  // Get the flux associated with this truth
685  //art::PtrVector<simb::MCFlux> fluxs = util::FindManyP<simb::MCFlux>(pv, evt, "geantgen", 0);
686  art::FindManyP<simb::MCFlux> fmf(pv, evt, fMCFluxLabel);
687  std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
688  simb::MCFlux flux;
689  if(fluxs.size() == 1){
690  flux = *fluxs[0];
691  }
692  else{
693  std::cout<<"There is no assoicated neutrino intereaction !"<<std::endl;
694  }
695  const simb::MCNeutrino& nu = truth->GetNeutrino();
696 
697  //select nue CC only events
698  if( !(abs(flux.fntype)==12 && abs(nu.Nu().PdgCode())==12 ) ) continue;
699  if( !(nu.CCNC()==0) ) continue;
700 
701  trueE = nu.Nu().E();
702  truePt = nu.Pt();
703  trueTheta = nu.Theta();
704  ccnc = nu.CCNC();
705  PDG = nu.Nu().PdgCode();
706  origPDG = flux.fntype;
707  mode = nu.Mode();
708  isQE = nu.InteractionType() == simb::kQE ||
709  nu.InteractionType() == simb::kCCQE;
710  intType = nu.InteractionType();
711  hitNucl = nu.HitNuc();
712  HadX = nu.X();
713  HadY = nu.Y();
714  HadW = nu.W();
715  qsqr = nu.QSqr();
716 
717  Vx=nu.Nu().Vx();
718  Vy=nu.Nu().Vy();
719  Vz=nu.Nu().Vz();
720 
721  if( fabs(nu.Nu().Vx())>200. ) continue;
722  if( fabs(nu.Nu().Vy())>200. ) continue;
723  if( nu.Nu().Vz()<0. ) continue;
724 
725  lepPDG = nu.Lepton().PdgCode();
726  lepE = nu.Lepton().E();
727  lepPx = nu.Lepton().Px();
728  lepPy = nu.Lepton().Py();
729  lepPz = nu.Lepton().Pz();
730 
731  art::Handle<std::vector<simb::MCTruth> > mclist;
732  evt.getByLabel(fMCFluxLabel,mclist);
733  art::Handle< std::vector<simb::MCFlux> > fllist;
734  evt.getByLabel(fMCFluxLabel, fllist);
735  if( mclist->size()>0 && fllist->size()>0 ){
736  for( unsigned int i_intx=0; i_intx<mclist->size(); ++i_intx){
737  art::Ptr<simb::MCTruth> mctruth(mclist,i_intx);
738  const simb::MCNeutrino& mcnu = mctruth->GetNeutrino();
739  TLorentzVector Tslnu(nu.Nu().Px(),nu.Nu().Py(),
740  nu.Nu().Pz(),nu.Nu().E());
741  TLorentzVector Tmcnu(mcnu.Nu().Px(),mcnu.Nu().Py(),
742  mcnu.Nu().Pz(),mcnu.Nu().E());
743  if(mcnu.Nu().PdgCode()==nu.Nu().PdgCode() &&
744  mcnu.CCNC()==nu.CCNC() && mcnu.Mode()==nu.Mode() &&
745  mcnu.Lepton().PdgCode()==nu.Lepton().PdgCode() &&
746  Tslnu==Tmcnu ){
747  art::Ptr<simb::MCFlux> flux(fllist,i_intx);
748  mesonPDG = flux->ftptype;
749  mesonPx = flux->ftpx;
750  mesonPy = flux->ftpy;
751  mesonPz = flux->ftpz;
752  }
753  }
754  }
755 
756  //missing DCMs
757  sl_Ndeaddcm=ndeaddcms;
758 
759  /////////////////
760  //select vertex//
761  /////////////////
762  if (!(fmVertex.isValid())){
763  hasVertex=0;
764  }//no reco vertex
765  else{
766  hasVertex=1;
767  std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
768  if( vert.size() != 1 ){
769  hasOneVertex=0;
770  }// ! 1-vertex
771  else{
772  hasOneVertex=1;
773  ///////////////////////////////////
774  //select the longest kalman track//
775  ///////////////////////////////////
776  std::vector<art::Ptr<rb::Track>> Ktracks = fmKTrack.at(sliceIdx);
777  double longestKTrack = 0;
778  int longestKTrack_ncells=0;
779  double longestKTrack_energy=0.;
780  int longestKTrack_nplanes=0;
781  float longestKTrack_phi=0.;
782  float longestKTrack_theta=0.;
783  float longestKTrack_startx=0.;
784  float longestKTrack_starty=0.;
785  float longestKTrack_startz=0.;
786  float longestKTrack_stopx=0.;
787  float longestKTrack_stopy=0.;
788  float longestKTrack_stopz=0.;
789  int longestKTrack_is3D=0;
790 
791  if( fmKTrack.isValid() ){
792  for(unsigned int n = 0; n < Ktracks.size(); ++n){
793  //if( !(Ktracks[n]->Is3D()) ) continue;
794  //tracks3d.push_back(Ktracks[n]);
795 
796  const double L = Ktracks[n]->TotalLength();
797  if(L > longestKTrack){
798  longestKTrack = L;
799  longestKTrack_ncells=Ktracks[n]->NCell();
800  longestKTrack_energy=Ktracks[n]->TotalGeV();
801  longestKTrack_nplanes=Ktracks[n]->ExtentPlane();
802  TVector3 dir = Ktracks[n]->Dir();
803  longestKTrack_phi=dir.Phi();
804  longestKTrack_theta=dir.Theta();
805  TVector3 start = Ktracks[n]->Start();
806  TVector3 stop = Ktracks[n]->Stop();
807  longestKTrack_startx = start.X();
808  longestKTrack_starty = start.Y();
809  longestKTrack_startz = start.Z();
810  longestKTrack_stopx = stop.X();
811  longestKTrack_stopy = stop.Y();
812  longestKTrack_stopz = stop.Z();
813  if( Ktracks[n]->Is3D() ) longestKTrack_is3D=1;
814  }//for longest track
815  }//loop all tracks
816  }//has tracks
817 
818  //mip cells and reco energy for slicer
819  double nCells = slice.NCell();
820  int lengthPlanes = slice.ExtentPlane();
821  if( slice.ExtentPlane()>192 ) lengthPlanes=192;
822 
823  double recoEnergy=0.;
824  double nmip=0.;
825  double Ntotal=0.;
826  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
827  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
828  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, slice.W(chit.get())));
829 
830  if( !rhit.IsCalibrated() ) continue;
831  //std::cout<<"cell "<<icell<<" energy is "<<rhit.GeV()<<std::endl;
832  recoEnergy += rhit.GeV();
833 
834  if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
835  Ntotal++;
836  }
837 
838  sl_runnum = evt.run();
839  sl_subrunnum = evt.subRun();
840  sl_evtnum = evt.id().event();
841  sl_Nslices=nslice;
842 
843  NKtracks = Ktracks.size();
844  KtrackLength = longestKTrack;
845  KtrackEnergy = longestKTrack_energy;
846  KtrackNcells = longestKTrack_ncells;
847  KtrackEnergy_ratio = longestKTrack_energy/recoEnergy;
848  KtrackNcells_ratio = longestKTrack_ncells/nCells;
849  KtrackNplanes = longestKTrack_nplanes;
850  KtrackPhi = longestKTrack_phi;
851  KtrackTheta = longestKTrack_theta;
852  KtrackStartX = longestKTrack_startx;
853  KtrackStartY = longestKTrack_starty;
854  KtrackStartZ = longestKTrack_startz;
855  KtrackStopX = longestKTrack_stopx;
856  KtrackStopY = longestKTrack_stopy;
857  KtrackStopZ = longestKTrack_stopz;
858  Ktrack3D = longestKTrack_is3D;
859 
860  sliceIndex = sliceIdx;
861  Eslice = recoEnergy;
862  Ncells = nCells;
863  Nplanes = lengthPlanes;
864  Nmip=nmip;
865  Fmip=nmip/Ntotal;
866 
867  recoVx=vert[0]->GetX();
868  recoVy=vert[0]->GetY();
869  recoVz=vert[0]->GetZ();
870 
871  ////////////////////////////////
872  //select most-energetic shower//
873  ////////////////////////////////
874  art::FindManyP<rb::Prong> fmProng3D(vert, evt, art::InputTag(fProngLabel,fInstLabel3D));
875  std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
876 
877  int nprongs=SelectedProng.size();
878  if( nprongs>20 ) nprongs=20;
879  Nprongs=nprongs;
880  std::vector<double> ProngEnergy;
881  std::vector<double> ProngEnergyXView;
882  std::vector<double> ProngEnergyYView;
883  std::vector<double> ProngNCells;
884  std::vector<double> ProngNCellsXView;
885  std::vector<double> ProngNCellsYView;
886  std::vector<double> ProngPhi;
887  std::vector<double> ProngTheta;
888  std::vector<double> ProngLength;
889  std::vector<double> ProngStartX;
890  std::vector<double> ProngStartY;
891  std::vector<double> ProngStartZ;
892  std::vector<double> ProngStopX;
893  std::vector<double> ProngStopY;
894  std::vector<double> ProngStopZ;
895  std::vector<int> ProngNplanes;
896 
897  //energy balance between first two prongs
898  double ratio=1.;
899  //phi angle between first two prongs
900  double dphi=0.;
901  //index for most-energetic prong
902  int Index1=-1;
903  if( SelectedProng.size()>0 ){
904 
905  std::vector<int> prongIndex;
906  std::vector<float> pcellPECorr;
907 
908  for(int ip=0; ip<nprongs; ++ip ){
909  //reco direction
910  TVector3 dir1 = SelectedProng[ip]->Dir();
911  TVector3 start1 = SelectedProng[ip]->Start();
912  TVector3 Pstop = start1 + (SelectedProng[ip]->TotalLength())*dir1;
913 
914  //reco momentum
915  double Eprong=0.;
916  double EprongX=0.;
917  double EprongY=0.;
918 
919  int npcells=0;;
920  int nxpcells=0;
921  int nypcells=0;
922 
923  for( unsigned int icell=0; icell<SelectedProng[ip]->NCell(); ++ icell){
924  const art::Ptr<rb::CellHit>& chit = SelectedProng[ip]->Cell(icell);
925  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[ip]->W(chit.get())));
926 
927  if( !rhit.IsCalibrated() ) continue;
928  Eprong += rhit.GeV();
929  npcells++;
930 
931  if( chit->View() == geo::kX ){
932  EprongX += rhit.GeV();
933  ++nxpcells;
934  }
935  if( chit->View() == geo::kY ){
936  EprongY += rhit.GeV();
937  ++nypcells;
938  }
939  prongIndex.push_back(ip);
940  pcellPECorr.push_back(rhit.PECorr());
941  }//loop prong associated cells
942  ProngTheta.push_back(dir1.Theta());
943  ProngLength.push_back(SelectedProng[ip]->TotalLength());
944 
945  ProngStartX.push_back(start1.X());
946  ProngStartY.push_back(start1.Y());
947  ProngStartZ.push_back(start1.Z());
948  ProngStopX.push_back(Pstop.X());
949  ProngStopY.push_back(Pstop.Y());
950  ProngStopZ.push_back(Pstop.Z());
951 
952  ProngEnergy.push_back(Eprong);
953  ProngEnergyXView.push_back(EprongX);
954  ProngEnergyYView.push_back(EprongY);
955  ProngPhi.push_back(dir1.Phi());
956  ProngNCells.push_back(npcells);
957  ProngNCellsXView.push_back(nxpcells);
958  ProngNCellsYView.push_back(nypcells);
959 
960  ProngNplanes.push_back(SelectedProng[ip]->ExtentPlane());
961  }//loop prongs
962 
963  int Index2=-1;
964  double MaxEnergy1=-999.;
965  double MaxEnergy2=-999.;
966  for( int i=0;i<nprongs; ++i ){
967  if( MaxEnergy1<ProngEnergy[i] ){
968  MaxEnergy1=ProngEnergy[i];
969  Index1=i;
970  }
971  }//leading
972 
973  p1Phi=ProngPhi[Index1];
974  p1Theta=ProngTheta[Index1];
975  p1StartX=ProngStartX[Index1];
976  p1StartY=ProngStartY[Index1];
977  p1StartZ=ProngStartZ[Index1];
978  p1StopX=ProngStopX[Index1];
979  p1StopY=ProngStopY[Index1];
980  p1StopZ=ProngStopZ[Index1];
981  p1Length=ProngLength[Index1];
982 
983  p1Energy=ProngEnergy[Index1];
984  p1EnergyXView=ProngEnergyXView[Index1];
985  p1EnergyYView=ProngNCellsYView[Index1];
986  p1NCells=ProngNCells[Index1];
987  p1NCellsXView=ProngNCellsXView[Index1];
988  p1NCellsYView=ProngNCellsYView[Index1];
989  p1Nplanes=ProngNplanes[Index1];
990 
991  int p1_ncells=0;
992  int p1_mipcells=0;
993  std::vector<float> p1cellX;
994  std::vector<float> p1cellY;
995  std::vector<float> p1cellZ;
996  std::vector<float> p1cellE;
997  std::vector<int> p1cellPlane;
998  std::vector<int> p1cellCell;
999  std::vector<int> p1cellADC;
1000  std::vector<float> p1cellPE;
1001  std::vector<float> p1cellPECorr;
1002  std::vector<int> p1cellXView;
1003  std::vector<float> p1cellTNS;
1004  std::vector<float> p1cellW;
1005 
1006  for( unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
1007  const art::Ptr<rb::CellHit>& chit = SelectedProng[Index1]->Cell(icell);
1008  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[Index1]->W(chit.get())));
1009  double wx=SelectedProng[Index1]->W(chit.get());
1010  double wy=wx;
1011  if( !rhit.IsCalibrated() ) continue;
1012  ++p1_ncells;
1013  if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
1014 
1015  if( chit->View()==geo::kX ){
1016  p1cellW.push_back(wx);
1017  }//XZ view
1018  else{
1019  p1cellW.push_back(wy);
1020  }//YZ view
1021 
1022  p1cellX.push_back(rhit.X());
1023  p1cellY.push_back(rhit.Y());
1024  p1cellZ.push_back(rhit.Z());
1025  p1cellE.push_back(rhit.GeV());
1026 
1027  p1cellADC.push_back(chit->ADC());
1028  p1cellPlane.push_back(chit->Plane());
1029  p1cellCell.push_back(chit->Cell());
1030  p1cellPECorr.push_back(rhit.PECorr());
1031 
1032  }//loop energetic prong associated cells
1033  //if( p1_ncells>500 ){
1034  // std::cout<<"&&& crazy event with "<<p1_ncells<<" cells for leading prong &&&"<<std::endl;
1035  // std::cout<<evt.id().event()<<" "<<sliceIdx<<" "<<nu.Nu().E()<<std::endl;
1036  // }
1037  if( p1_ncells>500 ) p1_ncells=500;
1038  prong1cellNcell=p1_ncells;
1039  int minplane=999;
1040  for( int ic=0; ic<p1_ncells; ++ic ){
1041  prong1cellX[ic]=p1cellX[ic];
1042  prong1cellY[ic]=p1cellY[ic];
1043  prong1cellZ[ic]=p1cellZ[ic];
1044  prong1cellE[ic]=p1cellE[ic];
1045 
1046  prong1cellADC[ic]=p1cellADC[ic];
1047  prong1cellPlane[ic]=p1cellPlane[ic];
1048  prong1cellCell[ic]=p1cellCell[ic];
1049  prong1cellPECorr[ic]=p1cellPECorr[ic];
1050  prong1cellW[ic]=p1cellW[ic];
1051 
1052  if( minplane>p1cellPlane[ic] ){
1053  minplane=p1cellPlane[ic];
1054  }
1055  }//loop selected cells
1056  if( p1_ncells>0 ) p1Fmip=(float)p1_mipcells/p1_ncells;
1057 
1058  int Nplanes_p1=ProngNplanes[Index1];
1059 
1060  if( ProngEnergy[Index1]>0. ){
1061  std::vector<double> p1PlaneEnergy;
1062  p1PlaneEnergy.clear();
1063  double totalE_road20=0.;
1064  for( int ip=0; ip<Nplanes_p1; ++ip ){
1065  double eplane=0.;
1066 
1067  double cellX=0.;//energetic cell X in each plane
1068  double cellY=0.;//energetic cell Y in each plane
1069  double Ecell=0.;
1070  double EXcell=0.;
1071  double EYcell=0.;
1072 
1073  std::vector<double> Xxview;
1074  std::vector<double> Exview;
1075  std::vector<double> Yyview;
1076  std::vector<double> Eyview;
1077  Xxview.clear();
1078  Exview.clear();
1079  Yyview.clear();
1080  Eyview.clear();
1081 
1082  for( int ic=0; ic<p1_ncells; ++ic ){
1083  int Dplane=p1cellPlane[ic]-minplane;
1084  if( Dplane==ip ){
1085  eplane += p1cellE[ic];
1086 
1087  Ecell += p1cellE[ic];
1088  EXcell += p1cellE[ic]*p1cellX[ic];
1089  EYcell += p1cellE[ic]*p1cellY[ic];
1090 
1091  if( p1cellPlane[ic]%2 != 0 ){
1092  Xxview.push_back(p1cellX[ic]);
1093  Exview.push_back(p1cellE[ic]);
1094  }//in X view
1095  else {
1096  Yyview.push_back(p1cellY[ic]);
1097  Eyview.push_back(p1cellE[ic]);
1098  }//in Y view
1099  }
1100  }//loop leading-prong associated cells
1101 
1102  p1PlaneEnergy.push_back(eplane);
1103 
1104  if( Ecell>0.0 ){
1105  cellX = EXcell/Ecell;
1106  cellY = EYcell/Ecell;
1107  }
1108  int Nxviews=Xxview.size();
1109  int Nyviews=Yyview.size();
1110  if( Nxviews>0 ){
1111  for( int ixv=0; ixv<Nxviews; ++ixv ){
1112  if( fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road20+=Exview[ixv];
1113  }
1114  }//in X view
1115  if( Nyviews>0 ){
1116  for( int iyv=0; iyv<Nyviews; ++iyv ){
1117  if( fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road20+=Eyview[iyv];
1118  }
1119  }//in Y view
1120 
1121  }//loop leading-prong associated planes
1122 
1123  double efrac_6plane=-1.;
1124  double p1_Eplane[20]={0.};//energy deposit upto 20 planes
1125 
1126  int nplanes=p1PlaneEnergy.size();
1127  for( int ip=0; ip<nplanes; ++ip ){
1128  for( int j=0; j<20; ++j ){
1129  if( ip<=j ) p1_Eplane[j]+=p1PlaneEnergy[ip];
1130  }
1131 
1132  if( ip<nplanes-6 ){
1133  double e_frac6=(p1PlaneEnergy[ip]+p1PlaneEnergy[ip+1]+p1PlaneEnergy[ip+2]+p1PlaneEnergy[ip+3]+p1PlaneEnergy[ip+4]+p1PlaneEnergy[ip+5])/ProngEnergy[Index1];
1134  if( efrac_6plane<e_frac6 ) efrac_6plane=e_frac6;
1135  }
1136  else if( nplanes<=6 ) efrac_6plane=1.;
1137 
1138  }//end loop for p1PlaneEnergy
1139 
1140  double efrac_Eplane[20]={0.};
1141  for( int ip=0; ip<20; ++ip ){
1142  efrac_Eplane[ip] = p1_Eplane[ip]/ProngEnergy[Index1];
1143  }
1144 
1145  double efrac_plane2=0.;
1146  double efrac_plane3=0.;
1147  double efrac_plane4=0.;
1148 
1149  if( p1PlaneEnergy.size()>1 && ProngEnergy[Index1]>0. )
1150  efrac_plane2=p1PlaneEnergy[1]/ProngEnergy[Index1];
1151  if( p1PlaneEnergy.size()>2 && ProngEnergy[Index1]>0. )
1152  efrac_plane3=p1PlaneEnergy[2]/ProngEnergy[Index1];
1153  if( p1PlaneEnergy.size()>3 && ProngEnergy[Index1]>0. )
1154  efrac_plane4=p1PlaneEnergy[3]/ProngEnergy[Index1];
1155 
1156  emaxfrac_6p=efrac_6plane;
1157  efrac_10p=efrac_Eplane[9];
1158  efrac_p2=efrac_plane2;
1159  efrac_p3=efrac_plane3;
1160  efrac_p4=efrac_plane4;
1161  efrac_2sig=totalE_road20/ProngEnergy[Index1];
1162 
1163  }//leading-prong energy > 0
1164 
1165  //for sub-leading prong
1166  if( nprongs>1 ){
1167  for( int i=0;i<nprongs; ++i ){
1168  if( i==Index1 ) continue;
1169  if( MaxEnergy2<ProngEnergy[i] ){
1170  MaxEnergy2=ProngEnergy[i];
1171  Index2=i;
1172  }
1173  }//sub-leading
1174  double Eden=ProngEnergy[Index1]+ProngEnergy[Index2];
1175  double Enum=ProngEnergy[Index1]-ProngEnergy[Index2];
1176  if( Eden>0.0 ) ratio=Enum/Eden;
1177  dphi=fabs(ProngPhi[Index1]-ProngPhi[Index2]);
1178  if( dphi>2.*TMath::Pi() ) dphi=dphi-2.*TMath::Pi();
1179  if( dphi>TMath::Pi() ) dphi=2.*TMath::Pi()-dphi;
1180 
1181  //fraction of mip hits for sub-leading prong
1182  if( prongIndex.size()>0 ){
1183  double nmips_p2=0.;
1184  double ntotal_p2=0.;
1185  for( unsigned int ic=0; ic<prongIndex.size(); ++ic ){
1186  if( prongIndex[ic] != Index2 ) continue;
1187  ntotal_p2 += 1.;
1188  if( pcellPECorr[ic]>100. && pcellPECorr[ic]<240. )
1189  nmips_p2 += 1.;
1190  }//
1191  if( nmips_p2>0. ) p2Fmip = nmips_p2/ntotal_p2;
1192  }//
1193 
1194  }//has at least two prongs
1195 
1196  Ebalance=ratio;
1197  Dphi=dphi;
1198  }//has at least one prong
1199 
1200  //LID
1201  art::FindOneP<slid::EventLID> fslLID(slices, evt, flidLabel);
1202  cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(fslLID.at(sliceIdx));
1203  art::Ptr<slid::EventLID> elid = rlid.ref();
1204  if( elid ){
1205  lid = elid->Value();
1206  }
1207 
1208  }//has only one vertex
1209  }//has a valid reco vertex
1210 
1211  _otree->Fill();
1212  }//end loop of slice
1213 
1214  return;
1215  }//end analyze
1216 
1217  void Xeff::endJob()
1218  {
1219  art::ServiceHandle<art::TFileService> tfs;
1220  TH1* hPOT = tfs->make<TH1F>("hTotalPOT", ";; POT", 1, 0, 1);
1221  hPOT->Fill(.5, fTotalPOT);
1222  TH1* hSPILL = tfs->make<TH1F>("hTotalSpill",";; SPILL", 1, 0, 1);
1223  hSPILL->Fill(.5, fTotalSpill);
1224  }
1225 
1226 }//end namespace
1227 
1228 namespace ncs{
1229 
1230  DEFINE_ART_MODULE(Xeff);
1231 
1232 }
1233 
ncs::Xeff::fProngLabel
std::string fProngLabel
label for kalman tracks
Definition: Xeff_module.cc:95
simb::MCParticle::E
double E(const int i=0) const
Definition: MCParticle.h:232
xml_schema::id
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
ncs::Xeff::Nmip
float Nmip
Definition: Xeff_module.cc:181
rb::CellHit::TNS
float TNS() const
Definition: CellHit.h:46
ncs::Xeff::efrac_p3
float efrac_p3
Definition: Xeff_module.cc:211
ncs::Xeff::prong1cellCell
int prong1cellCell[500]
Definition: Xeff_module.cc:223
art::Event::subRun
SubRunNumber_t subRun() const
Definition: Event.h:72
BackTracker.h
back track the reconstruction to the simulation
rb::Cluster::W
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
simb::MCParticle::PdgCode
int PdgCode() const
Definition: MCParticle.h:211
simb::MCNeutrino::CCNC
int CCNC() const
Definition: MCNeutrino.h:148
simb::MCNeutrino::QSqr
double QSqr() const
Definition: MCNeutrino.h:157
simb::MCNeutrino::Theta
double Theta() const
angle between incoming and outgoing leptons, in radians
Definition: MCNeutrino.cxx:63
rb::Cluster::NCell
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
simb::MCTruth::GetNeutrino
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:76
ncs::Xeff::KtrackNplanes
float KtrackNplanes
Definition: Xeff_module.cc:168
simb::MCParticle::Py
double Py(const int i=0) const
Definition: MCParticle.h:230
ncs::Xeff::KtrackTheta
float KtrackTheta
Definition: Xeff_module.cc:169
ncs::Xeff::prong1cellPECorr
float prong1cellPECorr[500]
Definition: Xeff_module.cc:224
ncs::Xeff::prong1cellX
float prong1cellX[500]
Definition: Xeff_module.cc:217
ncs::Xeff::KtrackLength
float KtrackLength
Definition: Xeff_module.cc:163
simb::MCFlux::ftpx
double ftpx
Definition: MCFlux.h:79
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
ncs::Xeff::p1StopZ
float p1StopZ
Definition: Xeff_module.cc:198
ncs::Xeff::mesonPx
float mesonPx
Definition: Xeff_module.cc:142
ncs::Xeff::spillPOT
double spillPOT
Definition: Xeff_module.cc:108
ncs::Xeff::p1StartY
float p1StartY
Definition: Xeff_module.cc:193
art::SubRun::subRun
SubRunNumber_t subRun() const
Definition: SubRun.h:44
nprongs
int nprongs
Definition: FillTrainingTrees.h:26
cheat::BackTracker::SliceToNeutrinoInteractions
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...
Definition: BackTracker_service.cc:1062
ncs::Xeff::qsqr
float qsqr
Definition: Xeff_module.cc:134
ncs::Xeff::fSliceLabel
std::string fSliceLabel
Definition: Xeff_module.cc:93
ncs::Xeff::~Xeff
virtual ~Xeff()
Definition: Xeff_module.cc:253
ncs::Xeff::p1StopX
float p1StopX
Definition: Xeff_module.cc:196
mf::LogInfo
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
Definition: MessageLogger.h:214
rb::CellHit::Plane
unsigned short Plane() const
Definition: CellHit.h:39
ncs::Xeff::p1Theta
float p1Theta
Definition: Xeff_module.cc:191
rb::CellHit::View
geo::View_t View() const
Definition: CellHit.h:41
simb::MCFlux::ftptype
int ftptype
Definition: MCFlux.h:82
ncs::Xeff::fKTrackLabel
std::string fKTrackLabel
label for slices
Definition: Xeff_module.cc:94
ncs::Xeff::analyze
void analyze(const art::Event &evt)
Definition: Xeff_module.cc:438
p
const char * p
Definition: xmltok.h:285
simb::MCNeutrino::Nu
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:146
ncs::Xeff::trueTheta
float trueTheta
Definition: Xeff_module.cc:123
ncs::Xeff::efrac_2sig
float efrac_2sig
Definition: Xeff_module.cc:213
geo::kX
Vertical planes which measure X.
Definition: PlaneGeo.h:28
simb::MCNeutrino::Pt
double Pt() const
transverse momentum of interaction, in GeV/c
Definition: MCNeutrino.cxx:74
simb::MCParticle::Px
double Px(const int i=0) const
Definition: MCParticle.h:229
ncs::Xeff::prong1cellY
float prong1cellY[500]
Definition: Xeff_module.cc:218
ncs::Xeff::KtrackEnergy_ratio
float KtrackEnergy_ratio
Definition: Xeff_module.cc:167
calib::Calibrator::MakeRecoHit
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
Definition: Calibrator_service.cc:575
simb::MCNeutrino::HitNuc
int HitNuc() const
Definition: MCNeutrino.h:152
rb::Cluster
A collection of associated CellHits.
Definition: Cluster.h:47
ncs::Xeff::fCellHitLabel
std::string fCellHitLabel
label for vertex
Definition: Xeff_module.cc:100
abs
void abs(TH1 *hist)
Definition: wrong_sign_uncertainty.C:716
ncs::Xeff::lepPx
float lepPx
Definition: Xeff_module.cc:137
ncs::Xeff::lepPy
float lepPy
Definition: Xeff_module.cc:138
ip
TString ip
Definition: loadincs.C:5
ncs::Xeff::emaxfrac_6p
float emaxfrac_6p
Definition: Xeff_module.cc:208
ncs::Xeff::beginJob
void beginJob()
Definition: Xeff_module.cc:273
simb::kCCQE
charged current quasi-elastic
Definition: MCNeutrino.h:96
PandAna.reco_validation.prod5_pid_validation.ratio
def ratio(spec1, spec2, nbins, pot, binrange=(0, 1))
Definition: prod5_pid_validation.py:554
ncs::Xeff::p1StartX
float p1StartX
Definition: Xeff_module.cc:192
ncs::Xeff::prong1cellNcell
int prong1cellNcell
Definition: Xeff_module.cc:216
flux
Loaders::FluxType flux
Definition: plot_prediction_systs.C:48
geo::kY
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
ncs::Xeff::p1Phi
float p1Phi
Definition: Xeff_module.cc:190
MCFlux.h
object containing MC flux information
ncs::Xeff::prong1cellPlane
int prong1cellPlane[500]
Definition: Xeff_module.cc:222
ncs::Xeff::Dphi
float Dphi
Definition: Xeff_module.cc:187
rb::RecoHit
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
ncs::Xeff::lepE
float lepE
Definition: Xeff_module.cc:140
febshutoff_auto.start
start
Definition: febshutoff_auto.py:53
ncs::Xeff::reconfigure
void reconfigure(const fhicl::ParameterSet &p)
Definition: Xeff_module.cc:258
ncs::Xeff::p1EnergyYView
float p1EnergyYView
Definition: Xeff_module.cc:200
ncs::Xeff::fInstLabel2D
std::string fInstLabel2D
instance label for prongs 3D
Definition: Xeff_module.cc:97
ncs::Xeff::prong1cellW
float prong1cellW[500]
Definition: Xeff_module.cc:225
ncs::Xeff::fMCFluxLabel
std::string fMCFluxLabel
instance label for prongs 2D
Definition: Xeff_module.cc:98
rb::CellHit::Cell
unsigned short Cell() const
Definition: CellHit.h:40
ncs::Xeff::Ebalance
float Ebalance
Definition: Xeff_module.cc:186
ncs::Xeff::hasOneVertex
int hasOneVertex
Definition: Xeff_module.cc:151
ncs::Xeff::efrac_10p
float efrac_10p
Definition: Xeff_module.cc:209
simb::MCNeutrino::InteractionType
int InteractionType() const
Definition: MCNeutrino.h:150
ncs::Xeff::p1StartZ
float p1StartZ
Definition: Xeff_module.cc:194
ncs::DEFINE_ART_MODULE
DEFINE_ART_MODULE(ROCKMRE)
ncs::Xeff::KtrackStartY
float KtrackStartY
Definition: Xeff_module.cc:172
check_time_usage.float
float
Definition: check_time_usage.py:7
simb::MCNeutrino::Lepton
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:147
simb::MCNeutrino::W
double W() const
Definition: MCNeutrino.h:154
ncs::Xeff::KtrackStartZ
float KtrackStartZ
Definition: Xeff_module.cc:173
ncs::Xeff::Eslice
float Eslice
Definition: Xeff_module.cc:159
breakpoints::beginJob
void beginJob()
art::PtrVector::push_back
void push_back(Ptr< U > const &p)
Definition: PtrVector.h:441
ncs::Xeff::endSubRun
void endSubRun(const art::SubRun &sr)
Definition: Xeff_module.cc:412
CLHEP::L
static constexpr double L
Definition: SystemOfUnits.h:104
simb::MCNeutrino::Y
double Y() const
Definition: MCNeutrino.h:156
ncs::Xeff::p1Fmip
float p1Fmip
Definition: Xeff_module.cc:205
ncs::Xeff::KtrackNcells
int KtrackNcells
Definition: Xeff_module.cc:164
ncs::Xeff::p1EnergyXView
float p1EnergyXView
Definition: Xeff_module.cc:199
ncs::Xeff::mesonPy
float mesonPy
Definition: Xeff_module.cc:143
ncs::Xeff::sl_Ndeaddcm
int sl_Ndeaddcm
Definition: Xeff_module.cc:118
ncs::Xeff::KtrackStartX
float KtrackStartX
Definition: Xeff_module.cc:171
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
fhicl::ParameterSet::get
T get(std::string const &key) const
Definition: ParameterSet.h:231
ncs::Xeff::fTotalSpill
int fTotalSpill
Definition: Xeff_module.cc:84
evt
int evt
Definition: FillTrainingTrees.h:22
ncs::Xeff::KtrackPhi
float KtrackPhi
Definition: Xeff_module.cc:170
pot
#define pot
Definition: template_basic.C:49
ncs::Xeff::HadY
float HadY
Definition: Xeff_module.cc:132
ncs::Xeff::KtrackStopY
float KtrackStopY
Definition: Xeff_module.cc:175
sr
caf::StandardRecord * sr
Definition: FillTrainingTrees.h:18
simb::MCNeutrino::X
double X() const
Definition: MCNeutrino.h:155
calib::j
const double j
Definition: BetheBloch.cxx:29
ncs::Xeff::recoVy
float recoVy
Definition: Xeff_module.cc:153
simb::MCFlux::ftpz
double ftpz
Definition: MCFlux.h:81
Vertex.h
Vertex location in position and time.
rb
Perform a "2 point" Hough transform on a collection of hits.
Definition: dEdxCalculator.h:20
sim
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:16
ncs::Xeff::truePt
float truePt
Definition: Xeff_module.cc:121
ncs::Xeff::MakeRecoHit
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
Definition: Xeff_module.cc:428
ncs::Xeff::p1Length
float p1Length
Definition: Xeff_module.cc:195
simb
This class describes a particle created in the detector Monte Carlo simulation.
Definition: TrackEva_module.cc:41
om::cout
OStream cout
Definition: OStream.cxx:6
nplanes
int nplanes
Definition: geom.C:145
ncs::Xeff::fVertexLabel
std::string fVertexLabel
label for generator information
Definition: Xeff_module.cc:99
ncs::Xeff::OscWt
float OscWt
Definition: Xeff_module.cc:135
Event.h
art::EventID::event
EventNumber_t event() const
Definition: EventID.h:116
rb::Cluster::Cell
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
ncs::Xeff::Fmip
float Fmip
Definition: Xeff_module.cc:182
ncs::Xeff::trueE
float trueE
Definition: Xeff_module.cc:122
simb::MCParticle::Vx
double Vx(const int i=0) const
Definition: MCParticle.h:220
art::TFileDirectory::make
T * make(ARGS...args) const
Definition: TFileDirectory.h:82
rb::RecoHit::GeV
float GeV() const
Definition: RecoHit.cxx:69
ncs::Xeff::p1NCellsXView
int p1NCellsXView
Definition: Xeff_module.cc:202
ncs::Xeff::HadX
float HadX
Definition: Xeff_module.cc:131
rawdata::RawDigit::ADC
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
ncs::Xeff::sl_spillPOT
double sl_spillPOT
Definition: Xeff_module.cc:117
HoughResult.h
Data resulting from a Hough transform on the cell hit positions.
ncs::Xeff::p1Energy
float p1Energy
Definition: Xeff_module.cc:189
art::Ptr::get
T const * get() const
Definition: Ptr.h:321
ncs::Xeff::_otree
TTree * _otree
Definition: Xeff_module.cc:111
dir
TDirectory * dir
Definition: macro.C:5
simb::MCFlux::fntype
int fntype
Definition: MCFlux.h:51
art::DataViewImpl::getByLabel
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
ncs::Xeff::HadW
float HadW
Definition: Xeff_module.cc:133
ncs::Xeff::p1NCellsYView
int p1NCellsYView
Definition: Xeff_module.cc:203
sumdata::POTSum::goodspills
int goodspills
Definition: POTSum.h:31
ncs::Xeff::mesonPz
float mesonPz
Definition: Xeff_module.cc:144
ncs::Xeff::lepPz
float lepPz
Definition: Xeff_module.cc:139
hit
Definition: structs.h:12
geom
void geom(int which=0)
Definition: geom.C:163
simb::MCParticle::Pz
double Pz(const int i=0) const
Definition: MCParticle.h:231
rb::Cluster::ExtentPlane
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
ncs::Xeff::prong1cellE
float prong1cellE[500]
Definition: Xeff_module.cc:220
simb::MCParticle::Vz
double Vz(const int i=0) const
Definition: MCParticle.h:222
ncs::Xeff::endJob
void endJob()
Definition: Xeff_module.cc:1217
ncs::Xeff::prong1cellZ
float prong1cellZ[500]
Definition: Xeff_module.cc:219
ncs::Xeff::fInstLabel3D
std::string fInstLabel3D
label for prong
Definition: Xeff_module.cc:96
ncs::Xeff::efrac_p2
float efrac_p2
Definition: Xeff_module.cc:210
ncs::Xeff::p2Fmip
float p2Fmip
Definition: Xeff_module.cc:207
simb::MCFlux::ftpy
double ftpy
Definition: MCFlux.h:80
ncs::Xeff::sl_Nslices
int sl_Nslices
Definition: Xeff_module.cc:116
rb::Cluster::IsNoise
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
ncs::Xeff::fTotalPOT
double fTotalPOT
Definition: Xeff_module.cc:83
sumdata::POTSum::totgoodpot
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
ncs::Xeff::sliceIndex
int sliceIndex
Definition: Xeff_module.cc:160
ncs::Xeff::recoVx
float recoVx
Definition: Xeff_module.cc:152
art::Event::run
RunNumber_t run() const
Definition: Event.h:77
ncs::Xeff::prong1cellADC
int prong1cellADC[500]
Definition: Xeff_module.cc:221
ncs::Xeff::KtrackStopZ
float KtrackStopZ
Definition: Xeff_module.cc:176
ncs::Xeff::KtrackNcells_ratio
float KtrackNcells_ratio
Definition: Xeff_module.cc:166
ncs::Xeff::_btree
TTree * _btree
Definition: Xeff_module.cc:85
simb::MCNeutrino
Event generator information.
Definition: MCNeutrino.h:18
w
Float_t w
Definition: plot.C:20
simb::MCNeutrino::Mode
int Mode() const
Definition: MCNeutrino.h:149
ncs::Xeff::KtrackEnergy
float KtrackEnergy
Definition: Xeff_module.cc:165
ncs::Xeff::recoVz
float recoVz
Definition: Xeff_module.cc:154
ncs::Xeff::sl_subrunnum
int sl_subrunnum
Definition: Xeff_module.cc:114
ncs::Xeff::efrac_p4
float efrac_p4
Definition: Xeff_module.cc:212
art::Event::id
EventID id() const
Definition: Event.h:56
simb::MCParticle::Vy
double Vy(const int i=0) const
Definition: MCParticle.h:221
Geometry.h
Encapsulate the geometry of one entire detector (near, far, ndos)
ncs::Xeff::KtrackStopX
float KtrackStopX
Definition: Xeff_module.cc:174
ncs::Xeff::flidLabel
std::string flidLabel
label for cell hits
Definition: Xeff_module.cc:101
art::Handle::failedToGet
bool failedToGet() const
Definition: Handle.h:196
cet::maybe_ref
Definition: maybe_ref.h:35
ncs::Xeff::p1StopY
float p1StopY
Definition: Xeff_module.cc:197
string
enum BeamMode string
Generated on Mon Mar 1 2021 14:13:48 for NOvA by   doxygen 1.8.11