ROCKMRE_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"
8 
9 #include "Geometry/Geometry.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 
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
48 #include "fhiclcpp/ParameterSet.h"
54 
55 #include <cmath>
56 #include <vector>
57 #include <string>
58 #include <vector>
59 
60 //include TMVA
61 #include "TMVA/Tools.h"
62 #include "TMVA/Reader.h"
63 #include "TMVA/MethodCuts.h"
64 
65 class TVector3;
66 class TH1F;
67 class TTree;
68 
69 namespace rb{class Cluster; class Track;}
70 namespace sim{class Particle;}
71 namespace simb{class MCFlux; class MCTruth; class MCNeutrino;}
72 
73 namespace ncs {
74  class ROCKMRE : public art::EDAnalyzer {
75  public:
76  explicit ROCKMRE(fhicl::ParameterSet const& pset);
77  virtual ~ROCKMRE();
78 
79  void beginJob();
80  void analyze(const art::Event& evt);
81  void reconfigure(const fhicl::ParameterSet& p);
82  void endSubRun(const art::SubRun& sr);
83  void endJob();
84  double fTotalPOT;
86  TTree *_btree;
87 
88  private:
89 
90  const rb::RecoHit MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
91  double const& w);
92 
93  //folders to get objects from
94  std::string fSliceLabel; ///label for slices
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 fDataSpillLabel; //label for data spill
102  std::string flidLabel; //label for event LID
103  int isMC; //mode for MC
104 
105  //save event info to ntuples
106  int runnum;
108  int evtnum;
109  int Nslices;
110  double spillPOT;
112 
113  int Ndeaddcm;
114 
115  TTree *_otree;
116 
117  //true info for MC
118  double trueE;
119  int ccnc;
120  int PDG;
121  int origPDG;
122  int intType;
123  double Vx;
124  double Vy;
125  double Vz;
126 
127  //TTree *_evttree;
133  double sl_spillPOT;
135 
136  /*** Vertex ***/
137  float recoVx;//reco vertex X
138  float recoVy;//reco vertex Y
139  float recoVz;//reco vertex Z
140 
141  /*** slicer information ***/
142  int Ncells;//# of associate cells
143  int Nplanes;//# of extended planes
144  float recoE;//reco energy
145  float PlaneEnergy[192];//set the maximal plane # to 192
146  int MinPlane;//min plane index
147  int PlaneNcells[192];//# of cells in each plane
149 
150  /*** mip information ***/
151  float Nmip;//# of mip cells
152  float Fmip;//# of mip cells over Ncells
153 
154  /*** cell information ***/
155  int slcellNcell;//# of good cells
156  float slcellX[500];
157  float slcellY[500];
158  float slcellZ[500];
159  float slcellE[500];
160  int slcellADC[500];
161  int slcellPlane[500];
162  int slcellCell[500];
163  float slcellPECorr[500];
164  float slcellW[500];
165 
166  /*** prong information ***/
167  int Nprongs;
168  float Ebalance;
169  float Dphi;
170 
171  float prongEnergy[20];
172  float prongPhi[20];
173  float prongTheta[20];
174  float prongStartX[20];
175  float prongStartY[20];
176  float prongStartZ[20];
177  float prongLength[20];
178  float prongStopX[20];
179  float prongStopY[20];
180  float prongStopZ[20];
181  float prongEnergyXView[20];//energy for XView cells
182  float prongEnergyYView[20];//energy for YView cells
183  int prongNCells[20];//# of cells
184  int prongNCellsXView[20];//# of cells for XView
185  int prongNCellsYView[20];//# of cells for YView
186  int prongNplanes[20];
187 
188  //leading prong
190  float prong1Fmip;
191  float prong1cellX[500];
192  float prong1cellY[500];
193  float prong1cellZ[500];
194  float prong1cellE[500];
195  int prong1cellADC[500];
196  int prong1cellPlane[500];
197  int prong1cellCell[500];
198  float prong1cellPECorr[500];
199  float prong1cellW[500];
200 
201  //2D prongs
203  int prong2DInXView[20];//prong in X View or not
204  float prongEnergy2D[20];//energy
205  int prongNCells2D[20];//# of cells
206  float prongPhi2D[20];
207  float prongTheta2D[20];
208  float prongStartX2D[20];
209  float prongStartY2D[20];
210  float prongStartZ2D[20];
211  float prongLength2D[20];
212  float prongStopX2D[20];
213  float prongStopY2D[20];
214  float prongStopZ2D[20];
215  float Ebalance2D;
216 
217  float lid;
218  };
219 }
220 
221 namespace ncs{
222 
223  //......................................................................
224  ROCKMRE::ROCKMRE(fhicl::ParameterSet const& pset)
225  : EDAnalyzer(pset)
226  {
227  mf::LogInfo("ROCKMRE")<<__PRETTY_FUNCTION__<<"\n";
228  reconfigure(pset);
229  }
230 
231  //......................................................................
233  {
234  }
235 
236  //......................................................................
238  {
239  fSliceLabel = pset.get< std::string >("SliceLabel");
240  fProngLabel = pset.get< std::string >("ProngLabel");
241  fInstLabel3D = pset.get< std::string >("InstLabel3D");
242  fInstLabel2D = pset.get< std::string >("InstLabel2D");
243  fMCFluxLabel = pset.get< std::string >("MCFluxLabel");
244  fCellHitLabel = pset.get< std::string >("CellHitLabel");
245  fVertexLabel = pset.get< std::string >("VertexLabel");
246  flidLabel = pset.get< std::string >("lidLabel");
247 
248  isMC = pset.get<int>("MCevts");
249  fDataSpillLabel = pset.get<std::string> ("DataSpillLabel");
250  }
251 
252  //......................................................................
254  {
256 
257  fTotalPOT=0.;
258  fTotalSpill=0;
259 
260  _otree = tfs->make<TTree>("ROCKMRE","ROCKMRE particle");
261 
262  _otree->Branch("trueE",&trueE,"trueE/F");
263  _otree->Branch("ccnc",&ccnc,"ccnc/I");
264  _otree->Branch("PDG",&PDG,"PDG/I");
265  _otree->Branch("origPDG",&origPDG,"origPDG/I");
266  _otree->Branch("intType",&intType,"intType/I");
267  _otree->Branch("Vx",&Vx,"Vx/F");
268  _otree->Branch("Vy",&Vy,"Vy/F");
269  _otree->Branch("Vz",&Vz,"Vz/F");
270 
271  _otree->Branch("sl_runnum",&sl_runnum,"sl_runnum/I");
272  _otree->Branch("sl_subrunnum",&sl_subrunnum,"sl_subrunnum/I");
273  _otree->Branch("sl_evtnum",&sl_evtnum,"sl_evtnum/I");
274  _otree->Branch("sl_Nslices",&sl_Nslices,"sl_Nslices/I");
275  _otree->Branch("sl_goodSpill",&sl_goodSpill,"sl_goodSpill/I");
276  _otree->Branch("sl_spillPOT",&sl_spillPOT,"sl_spillPOT/D");
277 
278  _otree->Branch("sl_Ndeaddcm",&sl_Ndeaddcm,"sl_Ndeaddcm/I");
279 
280  _otree->Branch("recoVx",&recoVx,"recoVx/F");
281  _otree->Branch("recoVy",&recoVy,"recoVy/F");
282  _otree->Branch("recoVz",&recoVz,"recoVz/F");
283 
284  _otree->Branch("slcellNcell",&slcellNcell,"slcellNcell/I");
285  _otree->Branch("slcellX",&slcellX,"slcellX[slcellNcell]/F");
286  _otree->Branch("slcellY",&slcellY,"slcellY[slcellNcell]/F");
287  _otree->Branch("slcellZ",&slcellZ,"slcellZ[slcellNcell]/F");
288  _otree->Branch("slcellE",&slcellE,"slcellE[slcellNcell]/F");
289  _otree->Branch("slcellADC",&slcellADC,"slcellADC[slcellNcell]/I");
290  _otree->Branch("slcellPlane",&slcellPlane,"slcellPlane[slcellNcell]/I");
291  _otree->Branch("slcellCell",&slcellCell,"slcellCell[slcellNcell]/I");
292  _otree->Branch("slcellW",&slcellW,"slcellW[slcellNcell]/F");
293  _otree->Branch("slcellPECorr",&slcellPECorr,"slcellPECorr[slcellNcell]/F");
294 
295  _otree->Branch("Ncells",&Ncells,"Ncells/I");
296  _otree->Branch("Nplanes",&Nplanes,"Nplanes/I");
297  _otree->Branch("recoE",&recoE,"recoE/F");
298  _otree->Branch("MinPlane",&MinPlane,"MinPlane/I");
299  _otree->Branch("sliceIndex",&sliceIndex,"sliceIndex/I");
300 
301  _otree->Branch("Nmip",&Nmip,"Nmip/F");
302  _otree->Branch("Fmip",&Fmip,"Fmip/F");
303 
304  _otree->Branch("Nprongs",&Nprongs,"Nprongs/I");
305  _otree->Branch("Ebalance",&Ebalance,"Ebalance/F");
306  _otree->Branch("Dphi",&Dphi,"Dphi/F");
307 
308  _otree->Branch("prongLength",&prongLength,"prongLength[Nprongs]/F");
309  _otree->Branch("prongStopX",&prongStopX,"prongStopX[Nprongs]/F");
310  _otree->Branch("prongStopY",&prongStopY,"prongStopY[Nprongs]/F");
311  _otree->Branch("prongStopZ",&prongStopZ,"prongStopZ[Nprongs]/F");
312 
313  _otree->Branch("prongEnergy",&prongEnergy,"prongEnergy[Nprongs]/F");
314  _otree->Branch("prongPhi",&prongPhi,"prongPhi[Nprongs]/F");
315  _otree->Branch("prongTheta",&prongTheta,"prongTheta[Nprongs]/F");
316  _otree->Branch("prongStartX",&prongStartX,"prongStartX[Nprongs]/F");
317  _otree->Branch("prongStartY",&prongStartY,"prongStartY[Nprongs]/F");
318  _otree->Branch("prongStartZ",&prongStartZ,"prongStartZ[Nprongs]/F");
319  _otree->Branch("prongEnergyXView",&prongEnergyXView,"prongEnergyXView[Nprongs]/F");
320  _otree->Branch("prongEnergyYView",&prongEnergyYView,"prongEnergyYView[Nprongs]/F");
321  _otree->Branch("prongNCells",&prongNCells,"prongNCells[Nprongs]/I");
322  _otree->Branch("prongNCellsXView",&prongNCellsXView,"prongNCellsXView[Nprongs]/I");
323  _otree->Branch("prongNCellsYView",&prongNCellsYView,"prongNCellsYView[Nprongs]/I");
324  _otree->Branch("prongNplanes",&prongNplanes,"prongNplanes[Nprongs]/I");
325 
326  _otree->Branch("prong1cellNcell",&prong1cellNcell,"prong1cellNcell/I");
327  _otree->Branch("prong1Fmip",&prong1Fmip,"prong1Fmip/F");
328  _otree->Branch("prong1cellX",&prong1cellX,"prong1cellX[prong1cellNcell]/F");
329  _otree->Branch("prong1cellY",&prong1cellY,"prong1cellY[prong1cellNcell]/F");
330  _otree->Branch("prong1cellZ",&prong1cellZ,"prong1cellZ[prong1cellNcell]/F");
331  _otree->Branch("prong1cellE",&prong1cellE,"prong1cellE[prong1cellNcell]/F");
332  _otree->Branch("prong1cellADC",&prong1cellADC,"prong1cellADC[prong1cellNcell]/I");
333  _otree->Branch("prong1cellPlane",&prong1cellPlane,"prong1cellPlane[prong1cellNcell]/I");
334  _otree->Branch("prong1cellCell",&prong1cellCell,"prong1cellCell[prong1cellNcell]/I");
335  _otree->Branch("prong1cellW",&prong1cellW,"prong1cellW[prong1cellNcell]/F");
336  _otree->Branch("prong1cellPECorr",&prong1cellPECorr,"prong1cellPECorr[prong1cellNcell]/F");
337 
338  _otree->Branch("Nprongs2D",&Nprongs2D,"Nprongs2D/I");
339  _otree->Branch("prong2DInXView",&prong2DInXView,"prong2DInXView[Nprongs2D]/I");
340  _otree->Branch("prongEnergy2D",&prongEnergy2D,"prongEnergy2D[Nprongs2D]/F");
341  _otree->Branch("prongNCells2D",&prongNCells2D,"prongNCells2D[Nprongs2D]/I");
342  _otree->Branch("prongLength2D",&prongLength2D,"prongLength2D[Nprongs2D]/F");
343  _otree->Branch("prongStopX2D",&prongStopX2D,"prongStopX2D[Nprongs2D]/F");
344  _otree->Branch("prongStopY2D",&prongStopY2D,"prongStopY2D[Nprongs2D]/F");
345  _otree->Branch("prongStopZ2D",&prongStopZ2D,"prongStopZ2D[Nprongs2D]/F");
346  _otree->Branch("prongPhi2D",&prongPhi2D,"prongPhi2D[Nprongs2D]/F");
347  _otree->Branch("prongTheta2D",&prongTheta2D,"prongTheta2D[Nprongs2D]/F");
348  _otree->Branch("prongStartX2D",&prongStartX2D,"prongStartX2D[Nprongs2D]/F");
349  _otree->Branch("prongStartY2D",&prongStartY2D,"prongStartY2D[Nprongs2D]/F");
350  _otree->Branch("prongStartZ2D",&prongStartZ2D,"prongStartZ2D[Nprongs2D]/F");
351  _otree->Branch("Ebalance2D",&Ebalance2D,"Ebalance2D/F");
352 
353  _otree->Branch("lid",&lid,"lid/F");
354 
355  _btree=new TTree("beamInfo","beam info");
356  _btree->Branch("runnum",&runnum,"runnum/I");
357  _btree->Branch("subrunnum",&subrunnum,"subrunnum/I");
358  _btree->Branch("evtnum",&evtnum,"evtnum/I");
359  _btree->Branch("Nslices",&Nslices,"Nslices/I");
360  _btree->Branch("spillPOT",&spillPOT,"spillPOT/D");
361  _btree->Branch("goodSpill",&goodSpill,"goodSpill/I");
362 
363  _btree->Branch("Ndeaddcm",&Ndeaddcm,"Ndeaddcm/I");
364 
365  }
366 
367  //......................................................................
369  {
371  if( isMC ) sr.getByLabel(fMCFluxLabel, pot);
372  else sr.getByLabel(fDataSpillLabel, pot);
373 
374  if( !pot.failedToGet() ){
375  fTotalPOT += pot->totgoodpot;
376  fTotalSpill += pot->goodspills;
377 
378  std::cout << "POT in subrun " << sr.subRun()
379  << " = " << pot->totgoodpot <<" with "<<fTotalSpill<<" spills"<< std::endl;
380  }
381  }
382 
383  //......................................................................
385  double const& w)
386  {
388  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
389 
390  return rhit;
391  }
392 
393  //......................................................................
395  {
396 
397  lid=-1.;
398 
399  sl_runnum=-1;
400  sl_subrunnum=-1;
401  sl_evtnum=-1;
402  sl_Nslices=-1;
403  sl_goodSpill=-1;
404  sl_spillPOT=-1.;
405  sl_Ndeaddcm=-1;
406 
407  runnum=-1;
408  subrunnum=-1;
409  evtnum=-1;
410  Nslices=-1;
411  spillPOT=0.;
412  goodSpill=-1;
413  Ndeaddcm=-1;
414 
415  recoVx=-999.;
416  recoVy=-999.;
417  recoVz=-999.;
418 
419  Ncells=0;
420  Nplanes=0;
421  recoE=-1.;
422  MinPlane=-1;
423  sliceIndex=-1;
424  for( int i=0; i<192; ++i ){
425  PlaneEnergy[i] = 0.;
426  PlaneNcells[i] = 0;
427  }
428 
429  Nmip=-1.;
430  Fmip=-1.;
431 
432  Nprongs=0;
433  Ebalance=0.;
434  Dphi=-1.;
435 
436  for( int i=0; i<20; ++i ){
437  prongEnergy[i]=0.;
438  prongPhi[i]=-1.;
439  prongTheta[i]=-1.;
440  prongStartX[i]=-999.;
441  prongStartY[i]=-999.;
442  prongStartZ[i]=-999.;
443  prongStopX[i]=-999.;
444  prongStopY[i]=-999.;
445  prongStopZ[i]=-999.;
446  prongLength[i]=-1.;
447  prongEnergyXView[i]=0.;
448  prongEnergyYView[i]=0.;
449  prongNCells[i]=0;
450  prongNCellsXView[i]=0;
451  prongNCellsYView[i]=0;
452  prongNplanes[i]=-1;
453 
454  prong2DInXView[i]=0;
455  prongEnergy2D[i]=0.;
456  prongNCells2D[i]=0;
457  prongLength2D[i]=-1.;
458  prongPhi2D[i]=-1.;
459  prongTheta2D[i]=-1.;
460  prongStartX2D[i]=-999.;
461  prongStartY2D[i]=-999.;
462  prongStartZ2D[i]=-999.;
463  prongStopX2D[i]=-999.;
464  prongStopY2D[i]=-999.;
465  prongStopZ2D[i]=-999.;
466  }
467  Nprongs2D=0;
468  Ebalance2D=0.;
469 
470  slcellNcell=0;
471  prong1cellNcell=0;
472  prong1Fmip=-1.;
473  for( int i=0; i<500; ++i ){
474  slcellX[i]=-999.;
475  slcellY[i]=-999.;
476  slcellZ[i]=-999.;
477  slcellE[i]=-999.;
478  slcellADC[i]=0;
479  slcellPlane[i]=-1;
480  slcellCell[i]=-1;
481  slcellW[i]=-1.;
482  slcellPECorr[i]=-1.;
483 
484  prong1cellX[i]=-999.;
485  prong1cellY[i]=-999.;
486  prong1cellZ[i]=-999.;
487  prong1cellE[i]=-999.;
488  prong1cellADC[i]=0;
489  prong1cellPlane[i]=-1;
490  prong1cellCell[i]=-1;
491  prong1cellW[i]=-1.;
492  prong1cellPECorr[i]=-1.;
493 
494  }
495 
498 
499  //get data spill info
501  if( isMC ) evt.getByLabel(fMCFluxLabel, spillPot);
502  else evt.getByLabel(fDataSpillLabel, spillPot);
503 
504  if( spillPot.failedToGet() ) return;
505  spillPOT = spillPot->spillpot;
506  if( isMC ) spillPOT=spillPOT/1.0e+12;
507  if( spillPot->goodbeam ) goodSpill=1;
508  //std::cout<<"beam info: "<<goodSpill<<" "<<spillPOT<<std::endl;
509 
510  //event info
511  runnum = evt.run();
512  subrunnum = evt.subRun();
513  evtnum = evt.id().event();
514 
515  //get all hits info
517  evt.getByLabel("calhit", chits0);
518 
519  std::vector<int> sq_plane;
520  sq_plane.clear();
521  std::vector<int> sq_xview;
522  sq_xview.clear();
523  std::vector<float> sq_tns;
524  sq_tns.clear();
525  std::vector<int> sq_cell;
526  sq_cell.clear();
527 
528  int sq_ncells=chits0->size();
529  for(int i = 0; i < sq_ncells; ++i){
530  art::Ptr< rb::CellHit> hit(chits0, i);
531  double hittime=hit->TNS()/1000.;
532  sq_plane.push_back(hit->Plane());
533  sq_tns.push_back(hittime);
534  sq_cell.push_back(hit->Cell());
535  if( (hit->View())==geo::kX ) sq_xview.push_back(1);
536  else sq_xview.push_back(0);
537  }//loop all hits
538 
539  //find dead dcms
540  int Ndbdcm_intime[14]={0};
541  int Ndbdcm_outtime[14]={0};
542 
543  for( int ic=0; ic<sq_ncells; ++ic ){
544  if( sq_plane[ic]<64 ){
545  if( sq_xview[ic]==1 ){
546  if( sq_cell[ic]>63 ){
547  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[0]++;
548  else Ndbdcm_outtime[0]++;
549  }//dcm 1
550  else{
551  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[1]++;
552  else Ndbdcm_outtime[1]++;
553  }//dcm 2
554  }//XZ view
555  else {
556  if( sq_cell[ic]>31 ){
557  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[2]++;
558  else Ndbdcm_outtime[2]++;
559  }//dcm 3
560  else {
561  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[3]++;
562  else Ndbdcm_outtime[3]++;
563  }//dcm 4
564  }//YZ view
565  }//DB1
566  else if( sq_plane[ic]<128 ){
567  if( sq_xview[ic]==1 ){
568  if( sq_cell[ic]>63 ){
569  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[4]++;
570  else Ndbdcm_outtime[4]++;
571  }//dcm 1
572  else{
573  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[5]++;
574  else Ndbdcm_outtime[5]++;
575  }//dcm 2
576  }//XZ view
577  else {
578  if( sq_cell[ic]>31 ){
579  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[6]++;
580  else Ndbdcm_outtime[6]++;
581  }//dcm 3
582  else {
583  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[7]++;
584  else Ndbdcm_outtime[7]++;
585  }//dcm 4
586  }//YZ view
587  }//DB2
588  else if( sq_plane[ic]<192 ){
589  if( sq_xview[ic]==1 ){
590  if( sq_cell[ic]>63 ){
591  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[8]++;
592  else Ndbdcm_outtime[8]++;
593  }//dcm 1
594  else{
595  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[9]++;
596  else Ndbdcm_outtime[9]++;
597  }//dcm 2
598  }//XZ view
599  else {
600  if( sq_cell[ic]>31 ){
601  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[10]++;
602  else Ndbdcm_outtime[10]++;
603  }//dcm 3
604  else {
605  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[11]++;
606  else Ndbdcm_outtime[11]++;
607  }//dcm 4
608  }//YZ view
609  }//DB3
610  else {
611  if( sq_xview[ic]==1 ){
612  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[12]++;
613  else Ndbdcm_outtime[12]++;
614  }//XZ view
615  else {
616  if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[13]++;
617  else Ndbdcm_outtime[13]++;
618  }//YZ view
619  }//MC
620  }//loop all cells
621 
622  int ndeaddcms=0;
623  std::vector<std::string> dcmlocation;
624  dcmlocation.clear();
625  for( int id=0; id<14; ++id ){
626  if( Ndbdcm_outtime[id]==0 && Ndbdcm_intime[id]==0 ){
627  ndeaddcms++;
628  }
629  }//loop all 14 dcms associated cells
630 
631  Ndeaddcm=ndeaddcms;
632 
633  //get reco slicer info
635  evt.getByLabel(fSliceLabel, slices);
636  art::FindManyP<rb::Vertex> fmVertex(slices, evt, fVertexLabel);
637 
638  int nslice=0;
639  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
640 
641  const rb::Cluster& slice = (*slices)[sliceIdx];
642  if(slice.IsNoise()) continue;
643  ++nslice;
644  }//loop slices
645 
646  Nslices=nslice;
647 
648 
649  //std::cout<<"# of missing dcms: "<<ndeaddcms<<std::endl;
650  //std::cout<<"# of slices: "<<nslice<<std::endl;
651 
652  _btree->Fill();
653 
654  for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
655 
656  const rb::Cluster& slice = (*slices)[sliceIdx];
657  if(slice.IsNoise()) continue;
658 
659  if (!(fmVertex.isValid())) continue;
660  std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
661 
662  if (vert.size() != 1) continue;
663 
664  //std::cout<<"having vertex "<<std::endl;
665 
666  double nCells = slice.NCell();
667  int lengthPlanes = slice.ExtentPlane();
668  if( slice.ExtentPlane()>192 ) lengthPlanes=192;
669 
670  double recoEnergy=0.;
671  double nmip=0.;
672  double Ntotal=0.;
673  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
674  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
675  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, slice.W(chit.get())));
676 
677  if( !rhit.IsCalibrated() ) continue;
678  recoEnergy += rhit.GeV();
679 
680  if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
681  Ntotal++;
682  }
683 
684  sl_runnum = evt.run();
685  sl_subrunnum = evt.subRun();
686  sl_evtnum = evt.id().event();
687  sl_Nslices=nslice;
688 
690  if( !isMC ){
693  }
694 
695  //std::cout<<"# of cells "<<nCells<<std::endl;
696  if( nCells>300 )
697  std::cout<<"xbbu "<<nCells<<" "<<evt.run()<<" "<<evt.subRun()<<" "<<evt.id().event()<<" "<<sliceIdx<<std::endl;
698  if( nCells<10 || nCells>500 ) continue;
699 
700 /*
701  if( isMC ){
702  art::Handle<std::vector<rb::CellHit> > chits0;
703  evt.getByLabel(fCellHitLabel, chits0);
704  art::ServiceHandle<cheat::BackTracker> bt;
705  art::PtrVector<rb::CellHit> allhits;
706  for(size_t h = 0; h < chits0->size(); ++h)
707  allhits.push_back(art::Ptr<rb::CellHit>(chits0, h));
708  const art::Ptr<rb::Cluster> allslice(slices,sliceIdx);
709  std::vector<cheat::NeutrinoEffPur> iacts = bt->SliceToNeutrinoInteractions(allslice, allhits);
710 
711  if( iacts.empty() ){
712  std::cout<<"no matching between slice and nu interactions!"<<std::endl;
713  return;
714  }
715  art::Ptr<simb::MCTruth>& truth = iacts[0].neutrinoInt;
716  art::PtrVector<simb::MCTruth> pv;
717  pv.push_back(truth);
718  art::FindManyP<simb::MCFlux> fmf(pv, evt, fMCFluxLabel);
719  std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
720  simb::MCFlux flux;
721  if(fluxs.size() == 1){
722  flux = *fluxs[0];
723  }
724  else{
725  std::cout<<"There is no assoicated neutrino intereaction !"<<std::endl;
726  }
727  const simb::MCNeutrino& nu = truth->GetNeutrino();
728 
729  trueE = nu.Nu().E();
730  ccnc = nu.CCNC();
731  PDG = nu.Nu().PdgCode();
732  origPDG = flux.fntype;
733  intType = nu.InteractionType();
734  Vx=nu.Nu().Vx();
735  Vy=nu.Nu().Vy();
736  Vz=nu.Nu().Vz();
737  }//for MC only
738 */
739 
740  sliceIndex = sliceIdx;
741  recoE = recoEnergy;
742  Ncells = nCells;
743  Nplanes = lengthPlanes;
744  MinPlane = slice.MinPlane();
745  Nmip=nmip;
746  Fmip=nmip/Ntotal;
747 
748  recoVx=vert[0]->GetX();
749  recoVy=vert[0]->GetY();
750  recoVz=vert[0]->GetZ();
751 
752  art::FindOneP<slid::EventLID> fslLID(slices, evt, flidLabel);
753  cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(fslLID.at(sliceIdx));
754  art::Ptr<slid::EventLID> elid = rlid.ref();
755  if( elid ){
756  lid = elid->Value();
757  }
758 
759  //reco prongs
761  std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
762 
763  //std::cout<<"# of prongs "<<SelectedProng.size()<<std::endl;
764  if( SelectedProng.size()<1 ) continue;
765  int nprongs=SelectedProng.size();
766  if( nprongs>20 ) nprongs=20;
768  std::vector<double> ProngEnergy;
769  std::vector<double> ProngEnergyXView;
770  std::vector<double> ProngEnergyYView;
771  std::vector<double> ProngPhi;
772 
773 
774  for(int ip=0; ip<nprongs; ++ip ){
775  //reco direction
776  TVector3 dir1 = SelectedProng[ip]->Dir();
777  //recoEta[ip]=dir1.Eta();
778  prongPhi[ip]=dir1.Phi();
779  prongTheta[ip]=dir1.Theta();
780 
781  TVector3 start1 = SelectedProng[ip]->Start();
782  prongStartX[ip] = start1.X();
783  prongStartY[ip] = start1.Y();
784  prongStartZ[ip] = start1.Z();
785 
786  prongLength[ip] = SelectedProng[ip]->TotalLength();
787 
788  TVector3 Pstop = start1 + (SelectedProng[ip]->TotalLength())*dir1;
789  prongStopX[ip] = Pstop.X();
790  prongStopY[ip] = Pstop.Y();
791  prongStopZ[ip] = Pstop.Z();
792 
793  prongNplanes[ip] = SelectedProng[ip]->ExtentPlane();
794 
795  //reco momentum
796  double Eprong=0.;
797  double EprongX=0.;
798  double EprongY=0.;
799 
800  int npcells=0;;
801  int nxpcells=0;
802  int nypcells=0;
803 
804  for( unsigned int icell=0; icell<SelectedProng[ip]->NCell(); ++ icell){
805  const art::Ptr<rb::CellHit>& chit = SelectedProng[ip]->Cell(icell);
806  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[ip]->W(chit.get())));
807 
808  if( !rhit.IsCalibrated() ) continue;
809  Eprong += rhit.GeV();
810  npcells++;
811 
812  if( chit->View() == geo::kX ){
813  EprongX += rhit.GeV();
814  ++nxpcells;
815  }
816  if( chit->View() == geo::kY ){
817  EprongY += rhit.GeV();
818  ++nypcells;
819  }
820 
821  }//loop prong associated cells
822  ProngEnergy.push_back(Eprong);
823  ProngEnergyXView.push_back(EprongX);
824  ProngEnergyYView.push_back(EprongY);
825  ProngPhi.push_back(dir1.Phi());
826  prongEnergy[ip]=Eprong;
827  prongEnergyXView[ip]=EprongX;
828  prongEnergyYView[ip]=EprongY;
829  prongNCells[ip]=npcells;
830  prongNCellsXView[ip]=nxpcells;
831  prongNCellsYView[ip]=nypcells;
832  }//loop prongs
833 
834  int Index1=-1;
835  double MaxEnergy1=-999.;
836  if( nprongs>0 ){
837  for( int i=0;i<nprongs; ++i ){
838  if( MaxEnergy1<ProngEnergy[i] ){
839  MaxEnergy1=ProngEnergy[i];
840  Index1=i;
841  }
842  }//leading
843  }
844 
845  if( Index1>=0 ){
846 
847  int p1_ncells=0;
848  int p1_mipcells=0;
849  std::vector<float> p1cellX;
850  std::vector<float> p1cellY;
851  std::vector<float> p1cellZ;
852  std::vector<float> p1cellE;
853  std::vector<int> p1cellPlane;
854  std::vector<int> p1cellCell;
855  std::vector<int> p1cellADC;
856  std::vector<float> p1cellPE;
857  std::vector<float> p1cellPECorr;
858  std::vector<int> p1cellXView;
859  std::vector<float> p1cellTNS;
860  std::vector<float> p1cellW;
861 
862  for( unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
863  const art::Ptr<rb::CellHit>& chit = SelectedProng[Index1]->Cell(icell);
864  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[Index1]->W(chit.get())));
865  double wx=SelectedProng[Index1]->W(chit.get());
866  double wy=wx;
867  if( !rhit.IsCalibrated() ) continue;
868  ++p1_ncells;
869  if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
870 
871  if( chit->View()==geo::kX ){
872  p1cellW.push_back(wx);
873  }//XZ view
874  else{
875  p1cellW.push_back(wy);
876  }//YZ view
877 
878  p1cellX.push_back(rhit.X());
879  p1cellY.push_back(rhit.Y());
880  p1cellZ.push_back(rhit.Z());
881  p1cellE.push_back(rhit.GeV());
882 
883  p1cellADC.push_back(chit->ADC());
884  p1cellPlane.push_back(chit->Plane());
885  p1cellCell.push_back(chit->Cell());
886  p1cellPECorr.push_back(rhit.PECorr());
887 
888  }//loop energetic prong associated cells
889  prong1cellNcell=p1_ncells;
890  if(p1_ncells>0 ) prong1Fmip = (float)p1_mipcells/p1_ncells;
891  for( int ic=0; ic<p1_ncells; ++ic ){
892  prong1cellX[ic]=p1cellX[ic];
893  prong1cellY[ic]=p1cellY[ic];
894  prong1cellZ[ic]=p1cellZ[ic];
895  prong1cellE[ic]=p1cellE[ic];
896 
897  prong1cellADC[ic]=p1cellADC[ic];
898  prong1cellPlane[ic]=p1cellPlane[ic];
899  prong1cellCell[ic]=p1cellCell[ic];
900  prong1cellPECorr[ic]=p1cellPECorr[ic];
901  prong1cellW[ic]=p1cellW[ic];
902  }//loop selected cells
903  }//for most energecti prong
904 
905  //slice associated cell information
906  int ncells_slice=0.;
907  std::vector<float> slicecellX;
908  std::vector<float> slicecellY;
909  std::vector<float> slicecellZ;
910  std::vector<float> slicecellE;
911  std::vector<int> slicecellPlane;
912  std::vector<int> slicecellCell;
913  std::vector<int> slicecellADC;
914  std::vector<float> slicecellPE;
915  std::vector<float> slicecellPECorr;
916  std::vector<int> slicecellXView;
917  std::vector<float> slicecellTNS;
918  std::vector<float> slicecellW;
919  std::vector<int> slicecellNoise;
920 
921  for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
922  const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
923  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, slice.W(chit.get())));
924 
925 
926  if( !rhit.IsCalibrated() ) continue;
927 
928  ncells_slice++;
929 
930  slicecellW.push_back(slice.W(chit.get()));
931 
932  slicecellX.push_back(rhit.X());
933  slicecellY.push_back(rhit.Y());
934  slicecellZ.push_back(rhit.Z());
935  slicecellE.push_back(rhit.GeV());
936 
937  slicecellADC.push_back(chit->ADC());
938  slicecellPlane.push_back(chit->Plane());
939  slicecellCell.push_back(chit->Cell());
940  slicecellPECorr.push_back(rhit.PECorr());
941 
942  }
943 
944  slcellNcell=ncells_slice;
945 
946  for( int ic=0; ic<ncells_slice; ++ic ){
947 
948  slcellX[ic]=slicecellX[ic];
949  slcellY[ic]=slicecellY[ic];
950  slcellZ[ic]=slicecellZ[ic];
951  slcellE[ic]=slicecellE[ic];
952 
953  slcellADC[ic]=slicecellADC[ic];
954  slcellPlane[ic]=slicecellPlane[ic];
955  slcellCell[ic]=slicecellCell[ic];
956  slcellPECorr[ic]=slicecellPECorr[ic];
957  slcellW[ic]=slicecellW[ic];
958 
959  }//loop selected cells
960 
961  //study 2D prongs
963  std::vector<art::Ptr<rb::Prong>> SelectedProng2D = fmProng2D.at(0);
964 
965  std::vector<int> Prong2DInXView;
966  std::vector<double> ProngEnergy2D;
967  int nprongs2D=SelectedProng2D.size();
968  if( nprongs2D>20 ) nprongs2D=20;
969  Nprongs2D = nprongs2D;
970  if( nprongs2D>0 ){
971  for(int ip=0; ip<nprongs2D; ++ip ){
972  prongNCells2D[ip]=SelectedProng2D[ip]->NCell();
973  int prongview=SelectedProng2D[ip]->View();
974 
975  prong2DInXView[ip]=prongview;//0 for XZ view
976  prongLength2D[ip]=SelectedProng2D[ip]->TotalLength();
977 
978  TVector3 dir1 = SelectedProng2D[ip]->Dir();
979  prongPhi2D[ip]=dir1.Phi();
980  prongTheta2D[ip]=dir1.Theta();
981 
982  TVector3 start1 = SelectedProng2D[ip]->Start();
983  prongStartX2D[ip] = start1.X();
984  prongStartY2D[ip] = start1.Y();
985  prongStartZ2D[ip] = start1.Z();
986 
987  TVector3 Pstop = start1 + (SelectedProng2D[ip]->TotalLength())*dir1;
988  prongStopX2D[ip] = Pstop.X();
989  prongStopY2D[ip] = Pstop.Y();
990  prongStopZ2D[ip] = Pstop.Z();
991 
992  Prong2DInXView.push_back(prongview);
993  //reco momentum
994  double Eprong=0.;
995  for( unsigned int icell=0; icell<SelectedProng2D[ip]->NCell(); ++ icell){
996  const art::Ptr<rb::CellHit>& chit = SelectedProng2D[ip]->Cell(icell);
997  const TVector3 mean = SelectedProng2D[ip]->MeanXYZ();
998  const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng2D[ip]->W(chit.get())));
999 
1000  if( !rhit.IsCalibrated() ) continue;
1001  Eprong += rhit.GeV();
1002  }
1003 
1004  prongEnergy2D[ip]=Eprong;
1005  ProngEnergy2D.push_back(Eprong);
1006  }
1007 
1008  int index3D=-1;
1009  double max3D=-99.;
1010  for( int i=0; i<nprongs; ++i ){
1011  if( max3D<ProngEnergy[i] ){
1012  max3D=ProngEnergy[i];
1013  index3D=i;
1014  }
1015  }//loop 3D prongs
1016  int index2D=-1;
1017  double max2D=-99.;
1018  for( int i=0;i<nprongs2D; ++i ){
1019  if( max2D<ProngEnergy2D[i] ){
1020  max2D=ProngEnergy2D[i];
1021  index2D=i;
1022  }
1023  }//loop 2D prongs
1024  double ebalance2D=1.;
1025  if( nprongs>0 && nprongs2D>0 ){
1026  double Eden=0.;
1027  double Enum=0.;
1028  if( Prong2DInXView[index2D]==0 ){//XZ view
1029  Eden=ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
1030  Enum=ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
1031  }
1032  else {//YZ view
1033  Eden=ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
1034  Enum=ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
1035  }
1036  if( Eden>0. ) ebalance2D=fabs(Enum)/Eden;
1037  }
1038  Ebalance2D=ebalance2D;
1039  }//having 2D prongs
1040 
1041  _otree->Fill();
1042 
1043  //_evttree->Fill();
1044  }//end loop of slice
1045 
1046  //_evttree->Fill();
1047  return;
1048  }//end analyze
1049 
1051  {
1053  TH1* hPOT = tfs->make<TH1F>("hTotalPOT", ";; POT", 1, 0, 1);
1054  hPOT->Fill(.5, fTotalPOT);
1055  TH1* hSPILL = tfs->make<TH1F>("hTotalSpill",";; SPILL", 1, 0, 1);
1056  hSPILL->Fill(.5, fTotalSpill);
1057  }
1058 
1059 }//end namespace
1060 
1061 namespace ncs{
1062 
1064 
1065 }
1066 
float prongStopY2D[20]
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
float TNS() const
Definition: CellHit.h:46
float prong1cellE[500]
void analyze(const art::Event &evt)
SubRunNumber_t subRun() const
Definition: Event.h:72
back track the reconstruction to the simulation
int prongNCells2D[20]
virtual double W(const rb::CellHit *chit) const
Estimate the unmeasured coordinate of chit.
Definition: Cluster.cxx:121
unsigned int NCell(geo::View_t view) const
Number of cells in view view.
Definition: Cluster.cxx:134
float prongPhi[20]
float slcellZ[500]
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
SubRunNumber_t subRun() const
Definition: SubRun.h:44
int slcellPlane[500]
int nprongs
int prong1cellCell[500]
float prong1cellPECorr[500]
float PlaneEnergy[192]
int prong1cellADC[500]
MaybeLogger_< ELseverityLevel::ELsev_info, false > LogInfo
unsigned short Plane() const
Definition: CellHit.h:39
float prongEnergyYView[20]
float prongTheta2D[20]
geo::View_t View() const
Definition: CellHit.h:41
const char * p
Definition: xmltok.h:285
std::string flidLabel
float prong1cellZ[500]
float slcellE[500]
Vertical planes which measure X.
Definition: PlaneGeo.h:28
rb::RecoHit MakeRecoHit(rb::CellHit const &cellhit, double w)
float prong1cellW[500]
A collection of associated CellHits.
Definition: Cluster.h:47
virtual ~ROCKMRE()
float prongStartZ2D[20]
TString ip
Definition: loadincs.C:5
std::string fInstLabel3D
label for prong
float prong1cellY[500]
float prongPhi2D[20]
std::string fProngLabel
label for slices
float prongStartZ[20]
int prongNplanes[20]
const rb::RecoHit MakeRecoHit(art::Ptr< rb::CellHit > const &chit, double const &w)
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
float slcellX[500]
object containing MC flux information
int prong1cellPlane[500]
Calibrated quantities relying on position in the orthogonal view. To generate a rb::CellHit from a rb...
Definition: RecoHit.h:19
unsigned short Cell() const
Definition: CellHit.h:40
float prongEnergyXView[20]
int prong2DInXView[20]
DEFINE_ART_MODULE(ROCKMRE)
float prongStopX[20]
void beginJob()
float prongStartX[20]
std::string fCellHitLabel
label for vertex
float prongEnergy2D[20]
float prongStopZ2D[20]
float prongStopX2D[20]
float prongStopZ[20]
T get(std::string const &key) const
Definition: ParameterSet.h:231
float slcellPECorr[500]
int evt
#define pot
caf::StandardRecord * sr
Vertex location in position and time.
Perform a "2 point" Hough transform on a collection of hits.
void reconfigure(const fhicl::ParameterSet &p)
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:15
std::string fSliceLabel
This class describes a particle created in the detector Monte Carlo simulation.
int prongNCellsXView[20]
OStream cout
Definition: OStream.cxx:6
std::string fVertexLabel
label for generator information
unsigned int MinPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.cxx:462
float slcellW[500]
EventNumber_t event() const
Definition: EventID.h:116
art::Ptr< rb::CellHit > Cell(geo::View_t view, unsigned int viewIdx) const
Get the ith cell from view view.
Definition: Cluster.cxx:145
float prong1cellX[500]
int prongNCells[20]
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
float GeV() const
Definition: RecoHit.cxx:69
float prongStartX2D[20]
int16_t ADC(uint32_t i) const
Definition: RawDigit.cxx:58
Data resulting from a Hough transform on the cell hit positions.
T const * get() const
Definition: Ptr.h:321
float prongTheta[20]
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
int goodspills
Definition: POTSum.h:31
Definition: structs.h:12
float prongEnergy[20]
void geom(int which=0)
Definition: geom.C:163
int slcellCell[500]
unsigned int ExtentPlane(geo::View_t view=geo::kXorY) const
Definition: Cluster.h:250
float prongStartY[20]
float slcellY[500]
int slcellADC[500]
std::string fDataSpillLabel
label for cell hits
float prongLength2D[20]
bool IsNoise() const
Is the noise flag set?
Definition: Cluster.h:163
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
float prongStartY2D[20]
RunNumber_t run() const
Definition: Event.h:77
std::string fInstLabel2D
instance label for prongs 3D
Float_t w
Definition: plot.C:20
float prongLength[20]
void endSubRun(const art::SubRun &sr)
double spillpot
POT for spill normalized by 10^12.
Definition: SpillData.h:26
EventID id() const
Definition: Event.h:56
float prongStopY[20]
Encapsulate the geometry of one entire detector (near, far, ndos)
int prongNCellsYView[20]
bool failedToGet() const
Definition: Handle.h:196
std::string fMCFluxLabel
instance label for prongs 2D
int PlaneNcells[192]