LArSoft  v07_13_02
Liquid Argon Software toolkit - http://larsoft.org/
TCShowerTemplateMaker_module.cc
Go to the documentation of this file.
1 // -------------------------------------------------
2 // makes shower profile templates
3 //
4 // Author: Rory Fitzpatrick (roryfitz@umich.edu)
5 // Created: 7/16/18
6 // -------------------------------------------------
7 
8 // Framework includes
16 #include "fhiclcpp/ParameterSet.h"
19 
37 
38 #include "TTree.h"
39 #include "TFile.h"
40 #include "TH1.h"
41 #include "TH3.h"
42 #include "TProfile.h"
43 #include "TProfile2D.h"
44 #include <TROOT.h>
45 #include <TStyle.h>
46 
47 namespace shower {
48 
50 
51  public:
52 
53  explicit TCShowerTemplateMaker(fhicl::ParameterSet const& pset);
54  virtual ~TCShowerTemplateMaker();
55 
56  void reconfigure(fhicl::ParameterSet const& pset);
57  void beginJob();
58  void analyze(const art::Event& evt);
59 
60  void showerProfile(std::vector< art::Ptr<recob::Hit> > showerhits, TVector3 shwvtx, TVector3 shwdir, double elep);
61  void showerProfileTrue(std::vector< art::Ptr<recob::Hit> > allhits, double elep);
63 
64  protected:
65 
66  // TTree* fTree;
70 
74 
78 
82 
86 
90 
94 
98 
102 
104  TH3F* fTransverse;
110 
111  // electrons
118  // photons
125  // other
132 
133  const int LBINS = 20;
134  const int LMIN = 0;
135  const int LMAX = 5;
136 
137  const int TBINS = 20;
138  const int TMIN = -5;
139  const int TMAX = 5;
140 
141  /*
142  const int EBINS = 10;
143  const double EMIN = 0.5;
144  const double EMAX = 10.5;
145  */
146 
147  const int EBINS = 20;
148  const double EMIN = 0.5;
149  const double EMAX = 20.5;
150 
151  const double X0 = 14;
152 
153  std::string fHitModuleLabel;
154  std::string fShowerModuleLabel;
155  std::string fGenieGenModuleLabel;
156  std::string fDigitModuleLabel;
157 
159 
160  }; // class TCShowerAnalysis
161 
162 } // shower
163 
164 // -------------------------------------------------
165 
167  EDAnalyzer(pset),
168  fHitModuleLabel (pset.get< std::string >("HitModuleLabel", "trajcluster" ) ),
169  fShowerModuleLabel (pset.get< std::string >("ShowerModuleLabel", "tcshower" ) ),
170  fGenieGenModuleLabel (pset.get< std::string >("GenieGenModuleLabel", "generator") ),
171  fDigitModuleLabel (pset.get< std::string >("DigitModuleLabel", "largeant") ),
172  fCalorimetryAlg (pset.get< fhicl::ParameterSet >("CalorimetryAlg") ) {
173  this->reconfigure(pset);
174 } // TCShowerTemplateMaker
175 
176 // -------------------------------------------------
177 
179 } // ~TCShowerTemplateMaker
180 
181 // -------------------------------------------------
182 
184 } // reconfigure
185 
186 // -------------------------------------------------
187 
189 
191 
192  // fTree = tfs->make<TTree>("tcshowerana", "tcshowerana");
193 
194  fShowerProfileSimLong = tfs->make<TProfile>("fShowerProfileSimLong", "longitudinal e- profile (true, simchannel);t;E (MeV)", LBINS, LMIN, LMAX);
195  fShowerProfileHitLong = tfs->make<TProfile>("fShowerProfileHitLong", "longitudinal e- profile (true, hit);t;E (MeV)", LBINS, LMIN, LMAX);
196  fShowerProfileRecoLong = tfs->make<TProfile>("fShowerProfileRecoLong", "longitudinal e- profile (reco);t;Q", LBINS, LMIN, LMAX);
197 
198  fShowerProfileSimLong2D = tfs->make<TProfile2D>("fShowerProfileSimLong2D", "longitudinal e- profile (true, simchannel);t;electron energy (MeV);E (MeV)", LBINS, LMIN, LMAX, EBINS, EMIN, EMAX);
199  fShowerProfileHitLong2D = tfs->make<TProfile2D>("fShowerProfileHitLong2D", "longitudinal e- profile (true, hit);t;electron energy (MeV);E (MeV)", LBINS, LMIN, LMAX, EBINS, EMIN, EMAX);
200  fShowerProfileRecoLong2D = tfs->make<TProfile2D>("fShowerProfileRecoLong2D", "longitudinal e- profile (reco);t;electron energy (MeV);Q", LBINS, LMIN, LMAX, EBINS, EMIN, EMAX);
201 
202  fShowerProfileSimTrans = tfs->make<TProfile>("fShowerProfileSimTrans", "transverse e- profile (true, simchannel);dist (cm);E (MeV)", TBINS, TMIN, TMAX);
203  fShowerProfileHitTrans = tfs->make<TProfile>("fShowerProfileHitTrans", "transverse e- profile (true, hit);dist (cm);E (MeV)", TBINS, TMIN, TMAX);
204  fShowerProfileRecoTrans = tfs->make<TProfile>("fShowerProfileRecoTrans", "transverse e- profile (reco);dist (cm);Q", TBINS, TMIN, TMAX);
205 
206  fShowerProfileSimTrans2D = tfs->make<TProfile2D>("fShowerProfileSimTrans2D", "transverse e- profile (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
207  fShowerProfileHitTrans2D = tfs->make<TProfile2D>("fShowerProfileHitTrans2D", "transverse e- profile (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
208  fShowerProfileRecoTrans2D = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D", "transverse e- profile (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
209 
210  fShowerProfileSimTrans2D_1 = tfs->make<TProfile2D>("fShowerProfileSimTrans2D_1", "transverse e- profile [0 <= t < 1] (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
211  fShowerProfileHitTrans2D_1 = tfs->make<TProfile2D>("fShowerProfileHitTrans2D_1", "transverse e- profile [0 <= t < 1] (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
212  fShowerProfileRecoTrans2D_1 = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_1", "transverse e- profile [0 <= t < 1] (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
213 
214  fShowerProfileSimTrans2D_2 = tfs->make<TProfile2D>("fShowerProfileSimTrans2D_2", "transverse e- profile [1 <= t < 2] (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
215  fShowerProfileHitTrans2D_2 = tfs->make<TProfile2D>("fShowerProfileHitTrans2D_2", "transverse e- profile [1 <= t < 2] (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
216  fShowerProfileRecoTrans2D_2 = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_2", "transverse e- profile [1 <= t < 2] (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
217 
218  fShowerProfileSimTrans2D_3 = tfs->make<TProfile2D>("fShowerProfileSimTrans2D_3", "transverse e- profile [2 <= t < 3] (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
219  fShowerProfileHitTrans2D_3 = tfs->make<TProfile2D>("fShowerProfileHitTrans2D_3", "transverse e- profile [2 <= t < 3] (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
220  fShowerProfileRecoTrans2D_3 = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_3", "transverse e- profile [2 <= t < 3] (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
221 
222  fShowerProfileSimTrans2D_4 = tfs->make<TProfile2D>("fShowerProfileSimTrans2D_4", "transverse e- profile [3 <= t < 4] (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
223  fShowerProfileHitTrans2D_4 = tfs->make<TProfile2D>("fShowerProfileHitTrans2D_4", "transverse e- profile [3 <= t < 4] (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
224  fShowerProfileRecoTrans2D_4 = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_4", "transverse e- profile [3 <= t < 4] (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
225 
226  fShowerProfileSimTrans2D_5 = tfs->make<TProfile2D>("fShowerProfileSimTrans2D_5", "transverse e- profile [4 <= t < 5] (true, simchannel);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
227  fShowerProfileHitTrans2D_5 = tfs->make<TProfile2D>("fShowerProfileHitTrans2D_5", "transverse e- profile [4 <= t < 5] (true, hit);t;electron energy (MeV);E (MeV)", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
228  fShowerProfileRecoTrans2D_5 = tfs->make<TProfile2D>("fShowerProfileRecoTrans2D_5", "transverse e- profile [4 <= t < 5] (reco);t;electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX);
229 
230  fLongitudinal = tfs->make<TH3F>("fLongitudinal", "longitudinal e- profile;t;electron energy (MeV);Q", LBINS, LMIN, LMAX, EBINS, EMIN, EMAX, 100, 0, 150000);
231  fTransverse = tfs->make<TH3F>("fTransverse", "transverse e- profile;dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 150000);
232  fTransverse_1 = tfs->make<TH3F>("fTransverse_1", "transverse e- profile [0 <= t < 1];dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 100000);
233  fTransverse_2 = tfs->make<TH3F>("fTransverse_2", "transverse e- profile [1 <= t < 2];dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 100000);
234  fTransverse_3 = tfs->make<TH3F>("fTransverse_3", "transverse e- profile [2 <= t < 3];dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 100000);
235  fTransverse_4 = tfs->make<TH3F>("fTransverse_4", "transverse e- profile [3 <= t < 4];dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 100000);
236  fTransverse_5 = tfs->make<TH3F>("fTransverse_5", "transverse e- profile [4 <= t < 5];dist (cm);electron energy (MeV);Q", TBINS, TMIN, TMAX, EBINS, EMIN, EMAX, 100, 0, 100000);
237 
238  // electrons
239  fLongitudinal_electron = tfs->make<TProfile>("fLongitudinal_electron", "longitudinal e- profile (reco);t;Q", LBINS, LMIN, LMAX);
240  fTransverse1_electron = tfs->make<TProfile>("fTransverse1_electron", "transverse e- profile [0 <= t < 1] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
241  fTransverse2_electron = tfs->make<TProfile>("fTransverse2_electron", "transverse e- profile [1 <= t < 2] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
242  fTransverse3_electron = tfs->make<TProfile>("fTransverse3_electron", "transverse e- profile [2 <= t < 3] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
243  fTransverse4_electron = tfs->make<TProfile>("fTransverse4_electron", "transverse e- profile [3 <= t < 4] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
244  fTransverse5_electron = tfs->make<TProfile>("fTransverse5_electron", "transverse e- profile [4 <= t < 5] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
245 
246  // photons
247  fLongitudinal_photon = tfs->make<TProfile>("fLongitudinal_photon", "longitudinal photon profile (reco);t;Q", LBINS, LMIN, LMAX);
248  fTransverse1_photon = tfs->make<TProfile>("fTransverse1_photon", "transverse photon profile [0 <= t < 1] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
249  fTransverse2_photon = tfs->make<TProfile>("fTransverse2_photon", "transverse photon profile [1 <= t < 2] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
250  fTransverse3_photon = tfs->make<TProfile>("fTransverse3_photon", "transverse photon profile [2 <= t < 3] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
251  fTransverse4_photon = tfs->make<TProfile>("fTransverse4_photon", "transverse photon profile [3 <= t < 4] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
252  fTransverse5_photon = tfs->make<TProfile>("fTransverse5_photon", "transverse photon profile [4 <= t < 5] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
253 
254  // other
255  fLongitudinal_other = tfs->make<TProfile>("fLongitudinal_other", "longitudinal other profile (reco);t;Q", LBINS, LMIN, LMAX);
256  fTransverse1_other = tfs->make<TProfile>("fTransverse1_other", "transverse other profile [0 <= t < 1] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
257  fTransverse2_other = tfs->make<TProfile>("fTransverse2_other", "transverse other profile [1 <= t < 2] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
258  fTransverse3_other = tfs->make<TProfile>("fTransverse3_other", "transverse other profile [2 <= t < 3] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
259  fTransverse4_other = tfs->make<TProfile>("fTransverse4_other", "transverse other profile [3 <= t < 4] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
260  fTransverse5_other = tfs->make<TProfile>("fTransverse5_other", "transverse other profile [4 <= t < 5] (reco);dist (cm);Q", TBINS, TMIN, TMAX);
261 
262 
263 } // beginJob
264 
265 // -------------------------------------------------
266 
268 
269  art::Handle< std::vector<recob::Hit> > hitListHandle;
270  std::vector<art::Ptr<recob::Hit> > hitlist;
271  if (evt.getByLabel(fHitModuleLabel,hitListHandle))
272  art::fill_ptr_vector(hitlist, hitListHandle);
273 
275  std::vector<art::Ptr<sim::SimChannel> > simchanlist;
276  if (evt.getByLabel(fDigitModuleLabel,scListHandle))
277  art::fill_ptr_vector(simchanlist, scListHandle);
278 
279  art::Handle< std::vector<recob::Shower> > showerListHandle;
280  std::vector<art::Ptr<recob::Shower> > showerlist;
281  if (evt.getByLabel(fShowerModuleLabel,showerListHandle))
282  art::fill_ptr_vector(showerlist, showerListHandle);
283 
284  art::Handle< std::vector<simb::MCTruth> > mctruthListHandle;
285  std::vector<art::Ptr<simb::MCTruth> > mclist;
286  if (evt.getByLabel(fGenieGenModuleLabel,mctruthListHandle))
287  art::fill_ptr_vector(mclist, mctruthListHandle);
288 
289  art::FindManyP<recob::Hit> shwfm(showerListHandle, evt, fShowerModuleLabel);
290 
291  if (mclist.size()) {
292  art::Ptr<simb::MCTruth> mctruth = mclist[0];
293  if (mctruth->NeutrinoSet()) {
294  if (std::abs(mctruth->GetNeutrino().Nu().PdgCode()) == 12 && mctruth->GetNeutrino().CCNC() == 0) {
295  double elep = mctruth->GetNeutrino().Lepton().E();
296  if (showerlist.size()) {
297  std::vector< art::Ptr<recob::Hit> > showerhits = shwfm.at(0);
298  showerProfile(showerhits, showerlist[0]->ShowerStart(), showerlist[0]->Direction(), elep);
299  }
300  showerProfileTrue(hitlist, elep);
301  showerProfileTrue(simchanlist, mctruth->GetNeutrino().Lepton());
302  }
303  // use this to look at background
304  /*
305  if (!(std::abs(mctruth->GetNeutrino().Nu().PdgCode()) == 12 && mctruth->GetNeutrino().CCNC() == 0)) {
306  if (showerlist.size()) {
307  std::vector< art::Ptr<recob::Hit> > showerhits = shwfm.at(0);
308  showerProfile(showerhits, showerlist[0]->ShowerStart(), showerlist[0]->Direction(), 0);
309  }
310  }
311  */
312 
313  }
314  }
315 
316  //fTree->Fill();
317 
318 } // analyze
319 
320 // -------------------------------------------------
321 
322 void shower::TCShowerTemplateMaker::showerProfile(std::vector< art::Ptr<recob::Hit> > showerhits, TVector3 shwvtx, TVector3 shwdir, double elep) {
323 
324  auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
326 
327  auto collectionPlane = geo::PlaneID(0, 0, 1);
328 
329  double shwVtxTime = detprop->ConvertXToTicks(shwvtx[0], collectionPlane);
330  double shwVtxWire = geom->WireCoordinate(shwvtx[1], shwvtx[2], collectionPlane);
331 
332  double shwTwoTime = detprop->ConvertXToTicks(shwvtx[0]+shwdir[0], collectionPlane);
333  double shwTwoWire = geom->WireCoordinate(shwvtx[1]+shwdir[1], shwvtx[2]+shwdir[2], collectionPlane);
334 
335  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
336  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
337 
338  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
339  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
340  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
341  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
342  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
343 
344  for (size_t i = 0; i < showerhits.size(); ++i) {
345  if (showerhits[i]->WireID().Plane != collectionPlane.Plane) continue;
346 
347  double wirePitch = geom->WirePitch(showerhits[i]->WireID());
348  double tickToDist = detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
349  tickToDist *= 1.e-3 * detprop->SamplingRate(); // 1e-3 is conversion of 1/us to 1/ns
350 
351  double xvtx = shwVtxTime * tickToDist;
352  double yvtx = shwVtxWire * wirePitch;
353 
354  double xtwo = shwTwoTime * tickToDist;
355  double ytwo = shwTwoWire * wirePitch;
356 
357  double xtwoorth = (ytwo - yvtx) + xvtx;
358  double ytwoorth = -(xtwo - xvtx) + yvtx;
359 
360  double xhit = showerhits[i]->PeakTime() * tickToDist;
361  double yhit = showerhits[i]->WireID().Wire * wirePitch;
362 
363  double ldist = std::abs((ytwoorth-yvtx)*xhit - (xtwoorth-xvtx)*yhit + xtwoorth*yvtx - ytwoorth*xvtx)/std::sqrt( pow((ytwoorth-yvtx), 2) + pow((xtwoorth-xvtx), 2) );
364  double tdist = ((ytwo-yvtx)*xhit - (xtwo-xvtx)*yhit + xtwo*yvtx - ytwo*xvtx)/std::sqrt( pow((ytwo-yvtx), 2) + pow((xtwo-xvtx), 2) );
365 
366  double to3D = 1. / sqrt( pow(xvtx-xtwo,2) + pow(yvtx-ytwo,2) ) ; // distance between two points in 3D space is one
367  ldist *= to3D;
368  tdist *= to3D;
369  double t = ldist/X0;
370 
371  double Q = showerhits[i]->Integral() * fCalorimetryAlg.LifetimeCorrection(showerhits[i]->PeakTime());
372 
373  ltemp->Fill(t, Q);
374  ttemp->Fill(tdist, Q);
375 
376  if (t < 1) ttemp_1->Fill(tdist, Q);
377  else if (t < 2) ttemp_2->Fill(tdist, Q);
378  else if (t < 3) ttemp_3->Fill(tdist, Q);
379  else if (t < 4) ttemp_4->Fill(tdist, Q);
380  else if (t < 5) ttemp_5->Fill(tdist, Q);
381 
382  } // loop through showerhits
383 
384  for (int i = 0; i < LBINS; ++i) {
385  if (ltemp->GetBinContent(i+1) == 0) continue;
386  fShowerProfileRecoLong->Fill(ltemp->GetBinCenter(i+1), ltemp->GetBinContent(i+1));
387  fShowerProfileRecoLong2D->Fill(ltemp->GetBinCenter(i+1), elep, ltemp->GetBinContent(i+1));
388  fLongitudinal->Fill(ltemp->GetBinCenter(i+1), elep, ltemp->GetBinContent(i+1));
389  }
390 
391  for (int i = 0; i < TBINS; ++i) {
392  if (ttemp->GetBinContent(i+1) == 0) continue;
393  fShowerProfileRecoTrans->Fill(ttemp->GetBinCenter(i+1), ttemp->GetBinContent(i+1));
394  fShowerProfileRecoTrans2D->Fill(ttemp->GetBinCenter(i+1), elep, ttemp->GetBinContent(i+1));
395  fTransverse->Fill(ttemp->GetBinCenter(i+1), elep, ttemp->GetBinContent(i+1));
396 
397  fShowerProfileRecoTrans2D_1->Fill(ttemp_1->GetBinCenter(i+1), elep, ttemp_1->GetBinContent(i+1));
398  fShowerProfileRecoTrans2D_2->Fill(ttemp_2->GetBinCenter(i+1), elep, ttemp_2->GetBinContent(i+1));
399  fShowerProfileRecoTrans2D_3->Fill(ttemp_3->GetBinCenter(i+1), elep, ttemp_3->GetBinContent(i+1));
400  fShowerProfileRecoTrans2D_4->Fill(ttemp_4->GetBinCenter(i+1), elep, ttemp_4->GetBinContent(i+1));
401  fShowerProfileRecoTrans2D_5->Fill(ttemp_5->GetBinCenter(i+1), elep, ttemp_5->GetBinContent(i+1));
402 
403  fTransverse_1->Fill(ttemp_1->GetBinCenter(i+1), elep, ttemp_1->GetBinContent(i+1));
404  fTransverse_2->Fill(ttemp_2->GetBinCenter(i+1), elep, ttemp_2->GetBinContent(i+1));
405  fTransverse_3->Fill(ttemp_3->GetBinCenter(i+1), elep, ttemp_3->GetBinContent(i+1));
406  fTransverse_4->Fill(ttemp_4->GetBinCenter(i+1), elep, ttemp_4->GetBinContent(i+1));
407  fTransverse_5->Fill(ttemp_5->GetBinCenter(i+1), elep, ttemp_5->GetBinContent(i+1));
408 
409  }
410 
411  return;
412 
413 } // showerProfile
414 
415 // -------------------------------------------------
416 
418 
419  auto const* detprop = lar::providerFrom<detinfo::DetectorPropertiesService>();
421  auto collectionPlane = geo::PlaneID(0, 0, 1);
424  std::map<int,double> trkID_E;
425 
426  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
427  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
428 
429  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
430  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
431  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
432  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
433  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
434 
435  double xvtx = -999;
436  double yvtx = -999;
437  double zvtx = -999;
438  double xtwo = -999;
439  double ytwo = -999;
440  double ztwo = -999;
441  double shwvtxT = -999;
442  double shwvtxW = -999;
443  double shwtwoT = -999;
444  double shwtwoW = -999;
445 
446  double shwvtxx = -999;
447  double shwvtxy = -999;
448  double shwtwox = -999;
449  double shwtwoy = -999;
450  double xtwoorth = -999;
451  double ytwoorth = -999;
452 
453  double wirePitch = -999;
454  double tickToDist = -999;
455 
456  bool foundParent = false;
457 
458  for (size_t i = 0; i < allhits.size(); ++i) {
459  if (allhits[i]->WireID().Plane != collectionPlane.Plane) continue;
460 
461  art::Ptr<recob::Hit> hit = allhits[i];
462  std::vector<sim::TrackIDE> trackIDs = btserv->HitToEveTrackIDEs(hit);
463 
464  for (size_t j = 0; j < trackIDs.size(); ++j) {
465  // only want energy associated with the electron and electron must have neutrino mother
466  if ( std::abs((piserv->TrackIdToParticle_P(trackIDs[j].trackID))->PdgCode()) != 11) continue;
467  if ( std::abs((piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Mother()) != 0) continue;
468 
469  if (!foundParent) {
470  xvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vx();
471  yvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vy();
472  zvtx = (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Vz();
473 
474  xtwo = xvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Px();
475  ytwo = yvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Py();
476  ztwo = zvtx + (piserv->TrackIdToParticle_P(trackIDs[j].trackID))->Pz();
477 
478  shwvtxT = detprop->ConvertXToTicks(xvtx, collectionPlane);
479  shwvtxW = geom->WireCoordinate(yvtx, zvtx, collectionPlane);
480 
481  shwtwoT = detprop->ConvertXToTicks(xtwo, collectionPlane);
482  shwtwoW = geom->WireCoordinate(ytwo, ztwo, collectionPlane);
483 
484  wirePitch = geom->WirePitch(hit->WireID());
485  tickToDist = detprop->DriftVelocity(detprop->Efield(),detprop->Temperature());
486  tickToDist *= 1.e-3 * detprop->SamplingRate(); // 1e-3 is conversion of 1/us to 1/ns
487 
488  shwvtxx = shwvtxT * tickToDist;
489  shwvtxy = shwvtxW * wirePitch;
490 
491  shwtwox = shwtwoT * tickToDist;
492  shwtwoy = shwtwoW * wirePitch;
493 
494  xtwoorth = (shwtwoy - shwvtxy) + shwvtxx;
495  ytwoorth = -(shwtwox - shwvtxx) + shwvtxy;
496 
497  foundParent = true;
498  }
499  double xhit = hit->PeakTime() * tickToDist;
500  double yhit = hit->WireID().Wire * wirePitch;
501 
502  double ldist = abs((ytwoorth-shwvtxy)*xhit - (xtwoorth-shwvtxx)*yhit + xtwoorth*shwvtxy - ytwoorth*shwvtxx)/sqrt( pow((ytwoorth-shwvtxy), 2) + pow((xtwoorth-shwvtxx), 2) );
503  double tdist = ((shwtwoy-shwvtxy)*xhit - (shwtwox-shwvtxx)*yhit + shwtwox*shwvtxy - shwtwoy*shwvtxx)/sqrt( pow((shwtwoy-shwvtxy), 2) + pow((shwtwox-shwvtxx), 2) );
504 
505  double to3D = sqrt( pow(xvtx-xtwo , 2) + pow(yvtx-ytwo , 2) + pow(zvtx-ztwo , 2) ) / sqrt( pow(shwvtxx-shwtwox,2) + pow(shwvtxy-shwtwoy,2) ) ; // distance between two points in 3D space is one
506  ldist *= to3D;
507  tdist *= to3D;
508  double t = ldist/X0;
509 
510  double energy = trackIDs[j].energy;
511 
512  ltemp->Fill(t, energy);
513  ttemp->Fill(tdist, energy);
514 
515  if (t < 1) ttemp_1->Fill(tdist, energy);
516  else if (t < 2) ttemp_2->Fill(tdist, energy);
517  else if (t < 3) ttemp_3->Fill(tdist, energy);
518  else if (t < 4) ttemp_4->Fill(tdist, energy);
519  else if (t < 5) ttemp_5->Fill(tdist, energy);
520 
521  } // loop through track IDE
522 
523  } // loop through all hits
524 
525  for (int i = 0; i < LBINS; ++i) {
526  if (ltemp->GetBinContent(i+1) == 0) continue;
527  fShowerProfileHitLong->Fill(ltemp->GetBinCenter(i+1), ltemp->GetBinContent(i+1));
528  fShowerProfileHitLong2D->Fill(ltemp->GetBinCenter(i+1), elep, ltemp->GetBinContent(i+1));
529  }
530 
531  for (int i = 0; i < TBINS; ++i) {
532  if (ttemp->GetBinContent(i+1) == 0) continue;
533  fShowerProfileHitTrans->Fill(ttemp->GetBinCenter(i+1), ttemp->GetBinContent(i+1));
534  fShowerProfileHitTrans2D->Fill(ttemp->GetBinCenter(i+1), elep, ttemp->GetBinContent(i+1));
535 
536  fShowerProfileHitTrans2D_1->Fill(ttemp_1->GetBinCenter(i+1), elep, ttemp_1->GetBinContent(i+1));
537  fShowerProfileHitTrans2D_2->Fill(ttemp_2->GetBinCenter(i+1), elep, ttemp_2->GetBinContent(i+1));
538  fShowerProfileHitTrans2D_3->Fill(ttemp_3->GetBinCenter(i+1), elep, ttemp_3->GetBinContent(i+1));
539  fShowerProfileHitTrans2D_4->Fill(ttemp_4->GetBinCenter(i+1), elep, ttemp_4->GetBinContent(i+1));
540  fShowerProfileHitTrans2D_5->Fill(ttemp_5->GetBinCenter(i+1), elep, ttemp_5->GetBinContent(i+1));
541  }
542 
543  return;
544 
545 } // showerProfileTrue
546 
547 // -------------------------------------------------
548 
552 
553  std::vector<sim::MCEnDep> alledep;
554 
555  TH1F* ltemp = new TH1F("ltemp", "ltemp", LBINS, LMIN, LMAX);
556  TH1F* ttemp = new TH1F("ttemp", "ttemp", TBINS, TMIN, TMAX);
557 
558  TH1F* ttemp_1 = new TH1F("ttemp_1", "ttemp_1", TBINS, TMIN, TMAX);
559  TH1F* ttemp_2 = new TH1F("ttemp_2", "ttemp_2", TBINS, TMIN, TMAX);
560  TH1F* ttemp_3 = new TH1F("ttemp_3", "ttemp_3", TBINS, TMIN, TMAX);
561  TH1F* ttemp_4 = new TH1F("ttemp_4", "ttemp_4", TBINS, TMIN, TMAX);
562  TH1F* ttemp_5 = new TH1F("ttemp_5", "ttemp_5", TBINS, TMIN, TMAX);
563 
564  // get electron energy depositions
565  for (size_t i = 0; i < allchan.size(); ++i) {
566  art::Ptr<sim::SimChannel> simchan = allchan[i];
567  if ( geom->View(simchan->Channel()) != geo::kV ) continue;
568  auto tdc_ide_map = simchan->TDCIDEMap();
569 
570  for(auto const& tdc_ide_pair : tdc_ide_map) {
571  for (auto const& ide : tdc_ide_pair.second) {
572  if (piserv->TrackIdToMotherParticle_P(ide.trackID) == NULL) continue;
573  if ( std::abs(piserv->TrackIdToMotherParticle_P(ide.trackID)->PdgCode()) != 11 ) continue;
574 
576  edep.SetVertex(ide.x,ide.y,ide.z);
577  edep.SetEnergy(ide.energy);
578  edep.SetTrackId(ide.trackID);
579 
580  alledep.push_back(edep);
581 
582  } // loop through ide
583  } // loop through tdc_ide
584 
585  } // loop through channels
586 
587  double x0 = electron.Vx();
588  double y0 = electron.Vy();
589  double z0 = electron.Vz();
590 
591  double x2 = electron.Px();
592  double y2 = electron.Py();
593  double z2 = electron.Pz();
594 
595  TVector3 v0(x2, y2, z2);
596  v0 = v0.Unit();
597 
598  for (size_t i = 0; i < alledep.size(); ++i) {
599  double x = (double)alledep[i].Vertex()[0];
600  double y = (double)alledep[i].Vertex()[1];
601  double z = (double)alledep[i].Vertex()[2];
602 
603  TVector3 v1(x-x0, y-y0, z-z0);
604 
605  double ldist = v0.Dot(v1);
606  double t = ldist/X0;
607  double tdist = (v0.Orthogonal()).Dot(v1);
608 
609  double energy = alledep[i].Energy();
610 
611  ltemp->Fill(t, energy);
612  ttemp->Fill(tdist, energy);
613 
614  if (t < 1) ttemp_1->Fill(tdist, energy);
615  else if (t < 2) ttemp_2->Fill(tdist, energy);
616  else if (t < 3) ttemp_3->Fill(tdist, energy);
617  else if (t < 4) ttemp_4->Fill(tdist, energy);
618  else if (t < 5) ttemp_5->Fill(tdist, energy);
619 
620  }
621 
622  for (int i = 0; i < LBINS; ++i) {
623  if (ltemp->GetBinContent(i+1) == 0) continue;
624  fShowerProfileSimLong->Fill(ltemp->GetBinCenter(i+1), ltemp->GetBinContent(i+1));
625  fShowerProfileSimLong2D->Fill(ltemp->GetBinCenter(i+1), electron.E(), ltemp->GetBinContent(i+1));
626  }
627 
628  for (int i = 0; i < TBINS; ++i) {
629  if (ttemp->GetBinContent(i+1) == 0) continue;
630  fShowerProfileSimTrans->Fill(ttemp->GetBinCenter(i+1), ttemp->GetBinContent(i+1));
631  fShowerProfileSimTrans2D->Fill(ttemp->GetBinCenter(i+1), electron.E(), ttemp->GetBinContent(i+1));
632 
633  fShowerProfileSimTrans2D_1->Fill(ttemp_1->GetBinCenter(i+1), electron.E(), ttemp_1->GetBinContent(i+1));
634  fShowerProfileSimTrans2D_2->Fill(ttemp_2->GetBinCenter(i+1), electron.E(), ttemp_2->GetBinContent(i+1));
635  fShowerProfileSimTrans2D_3->Fill(ttemp_3->GetBinCenter(i+1), electron.E(), ttemp_3->GetBinContent(i+1));
636  fShowerProfileSimTrans2D_4->Fill(ttemp_4->GetBinCenter(i+1), electron.E(), ttemp_4->GetBinContent(i+1));
637  fShowerProfileSimTrans2D_5->Fill(ttemp_5->GetBinCenter(i+1), electron.E(), ttemp_5->GetBinContent(i+1));
638  }
639 
640  return;
641 
642 } // showerProfileTrue
643 
644 // -------------------------------------------------
645 
Float_t x
Definition: compare.C:6
double E(const int i=0) const
Definition: MCParticle.h:237
const simb::MCParticle * TrackIdToParticle_P(int const &id)
TCShowerTemplateMaker(fhicl::ParameterSet const &pset)
int PdgCode() const
Definition: MCParticle.h:216
int CCNC() const
Definition: MCNeutrino.h:152
const simb::MCNeutrino & GetNeutrino() const
Definition: MCTruth.h:74
double Py(const int i=0) const
Definition: MCParticle.h:235
void SetEnergy(float e)
Definition: MCDataHolder.h:39
Planes which measure V.
Definition: geo_types.h:77
geo::WireID WireID() const
Initial tdc tick for hit.
Definition: Hit.h:234
Declaration of signal hit object.
Float_t y
Definition: compare.C:6
const simb::MCParticle & Nu() const
Definition: MCNeutrino.h:150
Double_t z
Definition: plot.C:279
The data type to uniquely identify a Plane.
Definition: geo_types.h:250
double Px(const int i=0) const
Definition: MCParticle.h:234
STL namespace.
void SetTrackId(unsigned int id)
Definition: MCDataHolder.h:41
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:313
object containing MC flux information
Float_t y2[n_points_geant4]
Definition: compare.C:26
auto vector(Vector const &v)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:265
const simb::MCParticle & Lepton() const
Definition: MCNeutrino.h:151
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:42
Provides recob::Track data product.
void showerProfile(std::vector< art::Ptr< recob::Hit > > showerhits, TVector3 shwvtx, TVector3 shwdir, double elep)
parameter set interface
double energy
Definition: plottest35.C:25
void SetVertex(float x, float y, float z)
Definition: MCDataHolder.h:32
EDAnalyzer(Table< Config > const &config)
Definition: EDAnalyzer.h:100
Declaration of cluster object.
Double_t edep
Definition: macro.C:13
View_t View(geo::PlaneID const &pid) const
Returns the view (wire orientation) on the channels of specified TPC plane.
Detector simulation of raw signals on wires.
raw::ChannelID_t Channel() const
Returns the readout channel this object describes.
Definition: SimChannel.h:332
double Vx(const int i=0) const
Definition: MCParticle.h:225
float PeakTime() const
Time of the signal peak, in tick units.
Definition: Hit.h:219
T * make(ARGS...args) const
Utility object to perform functions of association.
const std::vector< sim::TrackIDE > HitToEveTrackIDEs(recob::Hit const &hit)
bool getByLabel(std::string const &label, std::string const &productInstanceName, Handle< PROD > &result) const
Definition: DataViewImpl.h:344
double Pz(const int i=0) const
Definition: MCParticle.h:236
object containing MC truth information necessary for making RawDigits and doing back tracking ...
double Vz(const int i=0) const
Definition: MCParticle.h:227
TDCIDEs_t const & TDCIDEMap() const
Returns all the deposited energy information as stored.
Definition: SimChannel.h:331
Direction
Definition: AssnsIter.h:24
bool NeutrinoSet() const
Definition: MCTruth.h:75
TCEvent evt
Definition: DataStructs.cxx:5
void fill_ptr_vector(std::vector< Ptr< T >> &ptrs, H const &h)
Definition: Ptr.h:464
const simb::MCParticle * TrackIdToMotherParticle_P(int const &id)
Float_t x2[n_points_geant4]
Definition: compare.C:26
double LifetimeCorrection(double time, double T0=0) const
void showerProfileTrue(std::vector< art::Ptr< recob::Hit > > allhits, double elep)
double Vy(const int i=0) const
Definition: MCParticle.h:226
art framework interface to geometry description
void reconfigure(fhicl::ParameterSet const &pset)