RecoJMShowerAna_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 /// \brief Example to use RecoJMShower
3 /// for particle identification
4 /// \author Jianming Bian
5 ////////////////////////////////////////////////////////////////////////
6 
7 #include <cmath>
8 #include <iostream>
9 #include <fstream>
10 #include <string>
11 #include <vector>
12 
13 // ROOT includes
14 #include "TFile.h"
15 #include "TTree.h"
16 
17 // NOvA includes
18 #include "Calibrator/Calibrator.h"
19 #include "RecoJMShower/JMShower.h"
20 #include "SummaryData/POTSum.h"
21 
22 // ART includes
30 
31 
32 namespace jmshower {
33 
35  {
36  public:
37  explicit RecoJMShowerAna(fhicl::ParameterSet const& pset);
39  void analyze(const art::Event& evt);
40  void beginJob();
41  void beginSubRun(const art::SubRun& sr);
42  void reconfigure(const fhicl::ParameterSet& pset);
43 
44  std::string fPotLabel; ///< Module that produced the POTSum object
45 
46  private:
47 
48  TTree* fShower;
49  double m_showerRun;
51  double m_showerEvent;
52  double m_showerId;
53 
54 
60  double m_showerTrueNuP4[4];
61  double m_showerTrueNuVtx[3];
64 
68  double m_showerTrueVtx[3];
69  double m_showerTrueP4[4];
70  double m_showerTruePartHit[10];
72 
74  double m_showerNTrk;
79  double m_showerDir[3];
80  double m_showerStart[3];
81  double m_showerStop[3];
82  double m_showerPlaneEnergy[200];
83  double m_showerPlaneDedx[200];
88 
93  double m_showerDedxChisq[10];
94  double m_showerRadiusChisq[10];
95  double m_showerDedxNdof[10];
96  double m_showerDedxProb[10];
98  double m_showerDedxLLL[10];
99  double m_showerDedxLLL1[10];
100  double m_showerDedxLLL2[10];
101  double m_showerDedxLLL3[10];
102  double m_showerDedxLLL4[10];
103  double m_showerDedxLLL5[10];
104  double m_showerDedxLLL6[10];
105  double m_showerDedxLLT[10];
106  double m_showerPlaneRadius[200];
107  double m_showerADC;
110  //double m_showerEnergySvc;
113  double m_showerGap;
116 
117  };
118 }
119 
120 namespace jmshower{
121  //_________________________________________________________________________
122 
124  : EDAnalyzer(pset)
125  {
126  reconfigure(pset);
127  }
128 
129  //......................................................................
130 
132  {
133  fPotLabel =pset.get< std::string>("PotLabel");
134  }
135 
136  //......................................................................
137 
139  {
140  }
141 
142  //......................................................................
143 
145  {
146  std::cout<<"RecoJMShowerAna::beginJob()"<<std::endl;
148  fShower = tfs->make<TTree>("fShower","fShower");
149  fShower->Branch("showerRun",&m_showerRun);
150  fShower->Branch("showerSubRun",&m_showerSubRun);
151  fShower->Branch("showerEvent",&m_showerEvent);
152  fShower->Branch("showerTrueNuCCNC", &m_showerTrueNuCCNC);
153  fShower->Branch("showerTrueNuMode", &m_showerTrueNuMode);
154  fShower->Branch("showerTrueNuPdg", &m_showerTrueNuPdg);
155  fShower->Branch("showerTrueNuEnergy", &m_showerTrueNuEnergy);
156  fShower->Branch("showerTrueNuIsFid", &m_showerTrueNuIsFid);
157  fShower->Branch("showerTrueNuP4[4]", m_showerTrueNuP4);
158  fShower->Branch("showerTrueNuVtx[3]", m_showerTrueNuVtx);
159  fShower->Branch("showerEvtTrueEnergy[20]", m_showerEvtTrueEnergy);
160  fShower->Branch("showerEvtTruePdg[20]", m_showerEvtTruePdg);
161  fShower->Branch("showerId",&m_showerId);
162  fShower->Branch("showerNPlane", &m_showerNPlane);
163  fShower->Branch("showerNTrk", &m_showerNTrk);
164  fShower->Branch("showerNXcell", &m_showerNXcell);
165  fShower->Branch("showerNYcell", &m_showerNYcell);
166  fShower->Branch("showerDir[3]", m_showerDir);
167  fShower->Branch("showerStart[3]", m_showerStart);
168  fShower->Branch("showerStop[3]", m_showerStop);
169  fShower->Branch("showerADC", &m_showerADC);
170  fShower->Branch("showerEnergy0", &m_showerEnergy0);
171  fShower->Branch("showerEnergy", &m_showerEnergy);
172  fShower->Branch("showerNCell", &m_showerNCell);
173  fShower->Branch("showerRadius", &m_showerRadius);
174  fShower->Branch("showerGap", &m_showerGap);
175  fShower->Branch("showerIsFid", &m_showerIsFid);
176  fShower->Branch("showerNMIPPlane", &m_showerNMIPPlane);
177  fShower->Branch("showerContStartPlane", &m_showerContStartPlane);
178  fShower->Branch("showerRminmax", &m_showerRminmax);
179  fShower->Branch("showerPlaneEnergy[200]", m_showerPlaneEnergy);
180  fShower->Branch("showerPlaneDedx[200]", m_showerPlaneDedx);
181  fShower->Branch("showerDeltaPlaneDedx[200]", m_showerPlaneDeltaDedx);
182  fShower->Branch("showerDedxChisq[10]", m_showerDedxChisq);
183  fShower->Branch("showerRadiusChisq[10]", m_showerRadiusChisq);
184  fShower->Branch("showerDedxNdof[10]", m_showerDedxNdof);
185  fShower->Branch("showerDedxProb[10]", m_showerDedxProb);
186  fShower->Branch("showerTransDedxProb[10]", m_showerTransDedxProb);
187  fShower->Branch("showerDedxLLL[10]", m_showerDedxLLL);
188  fShower->Branch("showerDedxLLL1[10]", m_showerDedxLLL1);
189  fShower->Branch("showerDedxLLL2[10]", m_showerDedxLLL2);
190  fShower->Branch("showerDedxLLL3[10]", m_showerDedxLLL3);
191  fShower->Branch("showerDedxLLL4[10]", m_showerDedxLLL4);
192  fShower->Branch("showerDedxLLL5[10]", m_showerDedxLLL5);
193  fShower->Branch("showerDedxLLL6[10]", m_showerDedxLLL6);
194  fShower->Branch("showerDedxLLT[10]", m_showerDedxLLT);
195  fShower->Branch("showerPlaneRadius[200]", m_showerPlaneRadius);
196  fShower->Branch("showerPlaneDedxprobE[200]", m_showerPlaneDedxprobE);
197  fShower->Branch("showerPlaneDedxprobG[200]", m_showerPlaneDedxprobG);
198  fShower->Branch("showerPlaneDedxprobMu[200]", m_showerPlaneDedxprobMu);
199  fShower->Branch("showerPlaneDedxprobPi0[200]", m_showerPlaneDedxprobPi0);
200  fShower->Branch("showerTransCellDedx[20]", m_showerTransCellDedx);
201  fShower->Branch("showerTransCellProbE[20]", m_showerTransCellProbE);
202  fShower->Branch("showerTransCellProbPi0[20]", m_showerTransCellProbPi0);
203  fShower->Branch("showerTrueDang", &m_showerTrueDang);
204  fShower->Branch("showerTruePdg", &m_showerTruePdg);
205  fShower->Branch("showerTrueMother", &m_showerTrueMother);
206  fShower->Branch("showerTrueVtx[3]", m_showerTrueVtx);
207  fShower->Branch("showerTrueP4[4]", m_showerTrueP4);
208  fShower->Branch("showerTruePartHit[10]", m_showerTruePartHit);
209  fShower->Branch("showerTruePartEnergy[10]", m_showerTruePartEnergy);
210  }
211 
214  sr.getByLabel(fPotLabel,p);
215  std::cout<<"pot = "<<p->totgoodpot<<std::endl;
216  }
217  //......................................................................
219  {
220  //----------------------------------------------------
221  // Run, subrun and event information
222  //----------------------------------------------------
223  int run = evt.run();
224  int srun = evt.subRun();
225  int event = evt.id().event();
227  evt.getByLabel("recojmshower", showerhandle);
228  std::vector<art::Ptr<jmshower::JMShower> > svcshowercol;
229  for(unsigned int i = 0; i < showerhandle->size(); i++){
230  art::Ptr<jmshower::JMShower> p(showerhandle, i);
231  svcshowercol.push_back(p);
232  }
233 
234 
235 
236  //**********************************************************
237  // Fill showers information
238  //**********************************************************
239 
240  for(unsigned int i = 0; i < svcshowercol.size(); i++){
241  // General information
242  double showerRadius = svcshowercol[i]->Radius();
243  double showerEnergy = svcshowercol[i]->Energy();
244  double showerEnergy0 = svcshowercol[i]->DepositEnergy();
245  int showerIsFid = int(svcshowercol[i]->IsFiducial());
246  if(showerEnergy>0.001)showerRadius = showerRadius/showerEnergy;
247  m_showerNCell = svcshowercol[i]->NCell();
248  m_showerNPlane = svcshowercol[i]->ExtentPlane();
249  m_showerADC = svcshowercol[i]->TotalADC ();
251  m_showerEnergy0 = showerEnergy0;
252  m_showerDir[0] = (svcshowercol[i]->Dir())[0];
253  m_showerDir[1] = (svcshowercol[i]->Dir())[1];
254  m_showerDir[2] = (svcshowercol[i]->Dir())[2];
255  m_showerStart[0] = (svcshowercol[i]->Start())[0];
256  m_showerStart[1] = (svcshowercol[i]->Start())[1];
257  m_showerStart[2] = (svcshowercol[i]->Start())[2];
258  m_showerStop[0] = (svcshowercol[i]->Stop())[0];
259  m_showerStop[1] = (svcshowercol[i]->Stop())[1];
260  m_showerStop[2] = (svcshowercol[i]->Stop())[2];
261 
262 
263  m_showerRadius = showerRadius;
265  m_showerNMIPPlane = svcshowercol[i]->NMIPPlane();
266  // Fill longitudinal and transverse information
267 
268  for( int itype = int(jmshower::kNull)+1; itype != int(jmshower::kLast); ++itype)
269  { //type: 0-electron, 1-photon, 2-muon, 3-pion-zero, 4-inc-had, 5-proton, 6-neutron, 7-charged-pion;
270  m_showerDedxNdof[itype]=svcshowercol[i]->ExtentPlane();
271  m_showerDedxProb[itype] = svcshowercol[i]->DedxLongLL(itype);
272  m_showerTransDedxProb[itype] = svcshowercol[i]->DedxTransLL(itype);
273  m_showerDedxLLL[itype] = svcshowercol[i]->DedxLongLL(itype);
274  m_showerDedxLLT[itype] = svcshowercol[i]->DedxTransLL(itype);
275  }
276 
277  m_showerRun= run;
278  m_showerSubRun= srun;
279  m_showerEvent= event;
280  m_showerId= i;
281  m_showerNTrk= svcshowercol.size();
282  fShower->Fill();
283  }
284  }
285 }// end namespace jmshower
286 
287 namespace jmshower
288 {
290 }
SubRunNumber_t subRun() const
Definition: Event.h:72
void beginSubRun(const art::SubRun &sr)
const char * p
Definition: xmltok.h:285
DEFINE_ART_MODULE(TestTMapFile)
TODO.
Definition: FillPIDs.h:12
void analyze(const art::Event &evt)
bool IsFiducial(const TVector3 &nuVtx, const TVector3 &min, const TVector3 &max)
Definition: Flux.cxx:44
T get(std::string const &key) const
Definition: ParameterSet.h:231
double showerEnergy
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Definition: run.py:1
OStream cout
Definition: OStream.cxx:6
EventNumber_t event() const
Definition: EventID.h:116
void reconfigure(const fhicl::ParameterSet &pset)
RecoJMShowerAna(fhicl::ParameterSet const &pset)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double showerIsFid
double totgoodpot
normalized by 10^12 POT
Definition: POTSum.h:28
RunNumber_t run() const
Definition: Event.h:77
Definition: fwd.h:28
EventID id() const
Definition: Event.h:56
std::string fPotLabel
Module that produced the POTSum object.
static constexpr Double_t sr
Definition: Munits.h:164