TestBeamCommissioning_module.cc
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////
2 /// \file BeamlineDummyReco_module.cc
3 /// \module analyzer
4 /// \brief Makes commissioning plots for test beam to validate the
5 /// integrated (beamline+detector) system.
6 /// \authors Abhilash (Syracuse)
7 /// Mike Wallbank (University of Cincinnati) <wallbank@fnal.gov>
8 /// \date February 2019
9 ////////////////////////////////////////////////////////////////////////////
10 
11 // framework
19 #include "fhiclcpp/ParameterSet.h"
22 
23 // nova
24 #include "RawData/RawBeamline.h"
25 #include "RawData/RawDigit.h"
27 #include "BeamlineRecoBase/ToF.h"
28 #include "RecoBaseHit/CellHit.h"
31 
32 //Root includes
33 
34 #include "TH2.h"
35 
36 // -----------------------------------------------------------------------
37 namespace testbeam {
38 
39  // -----------------------------------------------------------------------
41 
42  public:
43 
45  // The compiler-generated destructor is fine for non-base
46  // classes without bare pointers or other resource use.
47 
48  // Required functions.
49  void analyze(art::Event const & e) override;
50 
51  // Selected optional functions.
52  void beginJob() override;
53  void endJob() override;
54  void reconfigure(const fhicl::ParameterSet& pset);
55 
56  private:
57 
58  TH1D* fPE;
59  TH1D* freconP;
60  TH1D* ftrueP;
61  TH1D* ftrueE;
62  TH2D* fMomentumHits; //< momentum [GeV]vs hits in detector
63  TH2D* fMomentumPE;
64  TH2D* ftruePhits;
65  TH2D* ftrueEhits;
66  TH2D* ftrueRecoP;
67  TH1D* fEndX;
68  TH1D* fEndY;
69  TH1D* fEndZ;
70  TH1D* fendplane;
71  TH2D* ftruePtof;
74 
75 
76  };
77  //------------------------------------------------------------------------
78 
81 
82  fPE = tfs->make<TH1D>("PE ","PE",100,0,15000);
83  freconP=tfs->make<TH1D>("RecoMomentum","RecoMomentum",100,700,1300);
84  ftrueP = tfs->make<TH1D>("trueP ","trueP",100,700,1300);
85  ftrueE=tfs->make<TH1D>("trueE","trueE",100,0,2000);
86 
87  fMomentumHits = tfs->make<TH2D>("reco momentum vs hits",";momentum vs hits;",500,700,1300,100,0,100);
88  fMomentumPE =tfs->make<TH2D>("recomomentum vs summed PE",";recomomentum vs summed PE;",500,700,1300,100,0,12000);
89  ftruePhits =tfs->make<TH2D>("true momentum vs hits ",";true momentum vs hits;",500,700,1300,100,0,100);
90  ftrueEhits =tfs->make<TH2D>("true energy vs hits",";true energy vs hits;",500,700,1800,100,0,100);
91  ftrueRecoP=tfs->make<TH2D>("true vs reco momentum",";true vs reco momentum;",500,700,1400,500,700,1400);
92  fEndX=tfs->make<TH1D>("End Position X",";End Position X;",100,0,10);
93  fEndY=tfs->make<TH1D>("End Position Y",";End Position Y;",100,0,10);
94  fEndZ=tfs->make<TH1D>("End Position Z",";End Position Z;",100,-10,0);
95  fendplane=tfs->make<TH1D>("End plane",";End plane;",65,0,65);
96  ftruePtof=tfs->make<TH2D>("True momentum vs tof","True momentum vs tof",500,700,1400,100,0,100);
97 
98  fMomentumHits->GetXaxis()->SetTitle("reco momentum (MeV)");
99  fMomentumHits->GetYaxis()->SetTitle("hits (RawDigitsize)");
100  fMomentumPE->GetXaxis()->SetTitle("Momentum (GeV)");
101  fMomentumPE->GetYaxis()->SetTitle("Total PE");
102  ftruePhits->GetXaxis()->SetTitle("True momentum (MeV)");
103  ftruePhits->GetYaxis()->SetTitle("hits (RawDigitsize)");
104 
105 
106 
107  }
108 
109 
110  // ---------------------------------------c--------------------------------
112  this->reconfigure(p);
113  }
114 
115  // -----------------------------------------------------------------------
117  fRawDigitLabel = pset.get<std::string>("RawDigitLabel");
118  fWCTrackLabel = pset.get<std::string>("WCTrackLabel");
119  }
120 
121  // -----------------------------------------------------------------------
123  {
125  std::vector<art::Ptr<simb::MCTruth> > mcTruths;
126  if (evt.getByLabel("generator", mcTruthHandle))
127  art::fill_ptr_vector(mcTruths, mcTruthHandle);
128 
129  std::cout << "There are " << mcTruths.size() << " mc truths in event " << evt.event() << std::endl;
130  int fparticles=0;
131  int fevent =0;
132  double totalPE=0;
133  double frecoP=0;
134  double truemomentum=0;
135  double trueenergy=0;
136  // int pdgid=0;
137  double endx=0;
138  double endy=0;
139  double endz=0;
140  int lastplane=0;
141  fevent=evt.event();
142  std::cout<<"event number "<<fevent<<std::endl;
143  double ftoftime=0;
144 
145  // print more stuff
146  std::cout << "There is " << mcTruths.size() << " MCTruth simulation-level objects in this event" << std::endl;
147 
148 // if (mcTruths.size()) {
149 // std::cout << "There are " << mcTruths[0]->NParticles() << " particles in this simulation" << std::endl;
150 // if (mcTruths[0]->NParticles()) {
151 // const simb::MCParticle& particle = mcTruths[0]->GetParticle(0);
152 // //std::cout << "Particle momentum " << particle.Momentum().Vect().Mag() << std::endl;
153 // std::cout << "Particle momentum " << particle.P() << std::endl;
154 // }
155 
156  for (std::vector<art::Ptr<simb::MCTruth> >::const_iterator mcTruthIt = mcTruths.begin(); mcTruthIt != mcTruths.end(); ++mcTruthIt)
157  {
158  std::cout << " There are " << (*mcTruthIt)->NParticles() << " particles in this MCTruth." << std::endl;
159  fparticles=(*mcTruthIt)->NParticles();
160  for(int j=0;j<(*mcTruthIt)->NParticles();j++)
161  {
162  simb::MCParticle particle=(*mcTruthIt)->simb::MCTruth::GetParticle(j);
163  std::cout<<"The particle "<<j<<" has momentum "<<particle.Pz()*1000<<std::endl;
164  std::cout<<"The particle "<<j<<" has energy "<<particle.E()*1000<<std::endl;
165  std::cout<<"The particle "<<j<<" has pdg code "<<particle.PdgCode()<<std::endl;
166  truemomentum=particle.Pz()*1000;
167  trueenergy=particle.E()*1000;
168  // pdgid=particle.PdgCode();
169  endx=particle.EndX();
170  endy=particle.EndY();
171  endz=particle.EndZ();
172  std::cout<<"The end z is "<<endz<<std::endl;
173  }
174 
175  }
176 
177 // art::Handle<std::vector<sim::Particle> > mcparticleHandle;
178 // std::vector<art::Ptr<sim::Particle> > mcparticle;
179 // if (evt.getByLabel("geantgen", mcparticleHandle))
180 // art::fill_ptr_vector(mcparticle, mcparticleHandle);
181 //
182 // for (std::vector<art::Ptr<simb::MCParticle> >::const_iterator mcparticleIt = mcparticle.begin(); mcparticleIt != mcparticle.end(); ++mcparticleIt)
183 // {
184 // //std::cout << " There are " << (*cellhitsIt)->PE() << " PE in this cellhit" << std::endl;
185 // truemomentum=(*mcparticleIt)->P();
186 // trueenergy=(*mcparticleIt)->E();
187 // std::cout<<"true momentum is "<<truemomentum<<std::endl;
188 // std::cout<<"true energy is "<<trueenergy<<std::endl;
189 // // fparticles=(*mcTruthIt)->NParticles();
190 // }
191 
192 
193  if (fparticles==1 ) {
194 
195  //truemomentum=mcparticle.P();
196  //trueenergy=mcparticle.E();
197  std::cout<<"true momentum is "<<truemomentum<<std::endl;
198  std::cout<<"true energy is "<<trueenergy<<std::endl;
199  //get PE from event
200 
202  std::vector<art::Ptr<rb::CellHit>> cellhits;
203  if(evt.getByLabel("calhit",cellhithandle))
204  {
205  art::fill_ptr_vector(cellhits,cellhithandle);
206  }
207 
208  //std::cout<<"there are "<<cellhits.PE()<<" photoelectrons"<<std::endl;
209  //std::cout<<"cell hit size "<< cellhithandle->size()<<std::endl;
210 
211  for (std::vector<art::Ptr<rb::CellHit> >::const_iterator cellhitsIt = cellhits.begin(); cellhitsIt != cellhits.end(); ++cellhitsIt)
212  {
213  // std::cout << " There are " << (*cellhitsIt)->PE() << " PE in this cellhit" << std::endl;
214  totalPE+=(*cellhitsIt)->PE();
215 // fparticles=(*mcTruthIt)->NParticles();
216  lastplane=(*cellhitsIt)->Plane();
217  std::cout<<" Current plane is "<<(*cellhitsIt)->Plane()<<std::endl;
218  }
219  fPE->Fill(totalPE);
220  //fPE->Fill(cellhits.PE());
221 
222 
224  std::vector<art::Ptr<brb::ToF>> toftime;
225  if(evt.getByLabel("beamlinesiminput",tofhandle))
226  {
227  art::fill_ptr_vector(toftime,tofhandle);
228  }
229 
230  for(std::vector<art::Ptr<brb::ToF>>::const_iterator tofIt=toftime.begin();tofIt !=toftime.end();++tofIt)
231  {
232  ftoftime=(*tofIt)->Time();
233  std::cout<<"tof time is "<<ftoftime<<std::endl;
234  }
235 
236 
237  // Get WC tracks from the event
239  std::vector<art::Ptr<brb::WCTrack> > wcTracks;
240  if (evt.getByLabel(fWCTrackLabel, wcTrackHandle))
241  art::fill_ptr_vector(wcTracks, wcTrackHandle);
242 
243 
244  // loop over the tracks
245  for (std::vector<art::Ptr<brb::WCTrack> >::const_iterator wcTrackIt = wcTracks.begin();wcTrackIt != wcTracks.end(); ++wcTrackIt)
246  {
247 // std::cout << "WC track has " << (*wcTrackIt)->NHits() << " and momentum " << (*wcTrackIt)->Momentum() << std::endl;
248  frecoP=(*wcTrackIt)->Momentum();
249  }
250 
251 
252 
253 
254  // Get the RawDigits from the event
256  evt.getByLabel(fRawDigitLabel, rawHits);
257  // How many are there in this event?
258  std::cout << "Event " << evt.event() << " has " << rawHits->size() << " raw digits." << std::endl;
259  // Now I have a rawdata::RawDigit vector and I can do whatever I want with it.
260  // let's loop over them all and get the ADC.
261  ftrueRecoP->Fill(truemomentum,frecoP);
262  ftrueP->Fill(truemomentum);
263  ftrueE->Fill(trueenergy);
264  fMomentumHits->Fill(frecoP,rawHits->size());
265  fMomentumPE->Fill(frecoP,totalPE);
266  freconP->Fill(frecoP);
267  ftruePhits->Fill(truemomentum,rawHits->size());
268  ftrueEhits->Fill(trueenergy,rawHits->size());
269  fEndX->Fill(endx);
270  fEndY->Fill(endy);
271  fEndZ->Fill(endz);
272  fendplane->Fill(lastplane);
273  ftruePtof->Fill(truemomentum,ftoftime);
274  }
275 
276 // // Get raw digits out of the event
277 // art::Handle<std::vector<rawdata::RawDigit> > rawDigitHandle;
278 // std::vector<art::Ptr<rawdata::RawDigit> > rawDigits;
279 // if (evt.getByLabel(fRawDigitLabel, rawDigitHandle))
280 // art::fill_ptr_vector(rawDigits, rawDigitHandle);
281 
282  // Print number of rawdigits
283  //std::cout << "Event " << evt.event() << " has " << rawHits->size() << " raw digits." << std::endl;
284 
285  }
286 
287 
288  // -----------------------------------------------------------------------
290  {
291  // Implementation of optional member function here.
292  }
293 
294 }
295 
296 // -----------------------------------------------------------------------
double E(const int i=0) const
Definition: MCParticle.h:232
int PdgCode() const
Definition: MCParticle.h:211
double EndZ() const
Definition: MCParticle.h:227
void analyze(art::Event const &e) override
const char * p
Definition: xmltok.h:285
DEFINE_ART_MODULE(TestTMapFile)
Particle class.
double EndY() const
Definition: MCParticle.h:226
T get(std::string const &key) const
Definition: ParameterSet.h:231
Encapsulation of reconstructed Time-of-Flight (ToF) information. Part of beamline reconstruction for ...
int evt
Encapsulation of reconstructed Wire Chamber track. Part of beamline reconstruction for NOvA test beam...
const double j
Definition: BetheBloch.cxx:29
EventNumber_t event() const
Definition: Event.h:67
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
TestBeamCommissioning(fhicl::ParameterSet const &p)
OStream cout
Definition: OStream.cxx:6
T * make(ARGS...args) const
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
void reconfigure(const fhicl::ParameterSet &pset)
double Pz(const int i=0) const
Definition: MCParticle.h:231
Raw data definitions for beamline data used in NOvA test beam experiment.
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
double EndX() const
Definition: MCParticle.h:225
Float_t e
Definition: plot.C:35
enum BeamMode string