FillTruth.cxx
Go to the documentation of this file.
1 #include <regex>
2 
3 #include "CAFMaker/FillTruth.h"
4 
5 #include "CAFMaker/Utils.h"
6 
7 #include "Geometry/Geometry.h"
8 #include "GeometryObjects/Geo.h"
10 #include "MCReweight/FluxWeights.h"
11 #include "G4Rwgt/G4WeightTable.h"
12 #include "RecoBase/CellHit.h"
13 #include "RecoBase/Cluster.h"
15 #include "Simulation/Particle.h"
19 
20 #include "TMatrixD.h"
21 #include "TVectorD.h"
22 
24 
25 namespace caf
26 {
27  //......................................................................
28  SRParticleTruth FillParticleTruth(const std::vector<rb::Cluster>& sliceList,
30  std::vector<cheat::TrackIDE>& allTracks,
31  std::vector<cheat::TrackIDE>& sliceTracks,
32  std::vector<cheat::TrackIDE>& allTracksBirks,
33  std::vector<cheat::TrackIDE>& sliceTracksBirks,
34  int sliceIdx)
35  {
36  SRParticleTruth truth;
38  cheat::ParticleEffPur puff = bt->ClusterToParticle(sliceList,
39  PtrVecToVec(hits),
40  sliceIdx);
41  truth.pdg = puff.pdg;
42  truth.eff = puff.efficiency;
43  truth.pur = puff.purity;
44  if(puff.trackID > -1 ){
45  const sim::Particle* part = bt->TrackIDToParticle(puff.trackID);
46  truth.trkID = puff.trackID;
47  truth.motherpdg = bt->TrackIDToMotherParticle(puff.trackID)->PdgCode();
48  truth.motherp = bt->TrackIDToMotherParticle(puff.trackID)->Momentum(0);
49  truth.p = part->Momentum(0);
50  truth.start = part->Position(0).Vect();
51 
52  for (int i=0; i<part->NumberDaughters();++i){
53 
54  truth.daughterlist.push_back(bt->TrackIDToParticle(part->Daughter(i))->PdgCode());
55  const sim::Particle * p = bt->TrackIDToParticle(part->Daughter(i));
56  truth.daughterVisEnergies.push_back(FindDaughterVisE(*p,allTracks));
57 
58  }
59 
60  int mId = part->Mother();
61  while (mId != 0){
62  truth.motherlist.push_back(bt->TrackIDToParticle(mId)->PdgCode());
63  mId = bt->TrackIDToParticle(mId)->Mother();
64  }
65 
66  int trackId = part->TrackId();
67  // Find visible energy
68  double visE = 0;
69  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
70  it != allTracks.end(); ++it)
71  {
72  if(trackId == it->trackID)
73  {
74  visE += it->energy;
75  break;
76  }
77  }
78  // Find visible energy in slice
79  double inSliceVisE = 0;
80  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracks.begin();
81  it != sliceTracks.end(); ++it)
82  {
83  if(trackId == it->trackID)
84  {
85  inSliceVisE += it->energy;
86  break;
87  }
88  }
89 
90  //now do visible energy with birks suppression
91  // Find visible energy
92  double visEBirks = 0;
93  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
94  it != allTracksBirks.end(); ++it)
95  {
96  if(trackId == it->trackID)
97  {
98  visEBirks += it->energy;
99  break;
100  }
101  }
102  // Find visible energy in slice
103  double inSliceVisEBirks = 0;
104  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
105  it != sliceTracksBirks.end(); ++it)
106  {
107  if(trackId == it->trackID)
108  {
109  inSliceVisEBirks += it->energy;
110  break;
111  }
112  }
113  //find visible energy of daughters
114  double daughterVisE = 0;
115  double daughterInSliceVisE = 0;
116  double daughterVisEBirks = 0;
117  double daughterInSliceVisEBirks = 0;
118  for (int i=0; i<part->NumberDaughters();++i){
119  const sim::Particle * p = bt->TrackIDToParticle(part->Daughter(i));
120  daughterVisE += FindDaughterVisE(*p, allTracks);
121  daughterInSliceVisE += FindDaughterVisE(*p, sliceTracks);
122  daughterVisEBirks += FindDaughterVisE(*p, allTracksBirks);
123  daughterInSliceVisEBirks += FindDaughterVisE(*p, sliceTracksBirks);
124  }
125 
126  truth.visE = visE;
127  truth.visEinslc = inSliceVisE;
128  truth.daughterVisE = daughterVisE;
129  truth.daughterVisEinslc = daughterInSliceVisE;
130  truth.visEBirks = visEBirks;
131  truth.visEinslcBirks = inSliceVisEBirks;
132  truth.daughterVisEBirks = daughterVisEBirks;
133  truth.daughterVisEinslcBirks = daughterInSliceVisEBirks;
134 
135 
136 //find the energy of each process that contributed to the prong to decide which was
137 //the dominant . Also keep indices to identify the particle and its energy that likely
138 //caused it. Only saving the particle of the MOST ENERGETIC fls hit.
139 
140  int tempFlsIdx = 0;
141  double tempBestE = 0;
142 
143  double bestE = -1;
144  int bestHitIdx = -1;
145  int bestFlsIdx = -1;
146 
147 
148  double inelasticProton = 0;
149  double inelasticPhoton = 0 ;
150  double elasticProton = 0;
151  std::map<std::string,double> processMap;
152  const sim::ParticleNavigator& nav = bt->ParticleNavigator();
153 
154  for (unsigned int idx = 0; idx < hits.size(); idx++){
155  const std::vector<sim::FLSHit> flshits = bt->HitToFLSHit(hits[idx]);
156 
157  for(unsigned int idx_fls = 0; idx_fls < flshits.size(); idx_fls++){
158  tempFlsIdx = idx_fls;
159  int flshit_trackid = flshits[idx_fls].GetTrackID();
160 
161  const sim::ParticleHistory particleHistory(&nav, flshit_trackid);
162  const sim::Particle *partHist = particleHistory[particleHistory.size()-1];
163  const sim::Particle *partHistFirst = particleHistory[0];
164 
165  if(idx_fls == 0 && partHistFirst->PdgCode() ==2112){
166  truth.primNeutronE.push_back(partHistFirst->E());
167  }
168 
169  processMap.insert(std::pair<std::string, double>("Primary", 0));
170  processMap.insert(std::pair<std::string, double>("hadElastic", 0));
171  processMap.insert(std::pair<std::string, double>("Decay", 0));
172  processMap.insert(std::pair<std::string, double>("pi-Inelastic", 0));
173  processMap.insert(std::pair<std::string, double>("protonInleastic", 0));
174  processMap.insert(std::pair<std::string, double>("neutronInelastic", 0));
175  processMap.insert(std::pair<std::string, double>("other", 0));
176 
177  tempBestE = flshits[idx_fls].GetEdep();
178  std::map<std::string, double>::iterator it;
179  it = processMap.find(partHist->Process());
180 
181  if(partHistFirst->PdgCode() == 2112){
182  if(partHist->Process() == "neutronInelastic"){
183  if(partHist->PdgCode() ==2212){
184  inelasticProton += flshits[idx_fls].GetEdep();
185  }
186  if(partHist->PdgCode() ==22){
187  inelasticPhoton += flshits[idx_fls].GetEdep();
188  }
189 
190  }
191  if(partHist->Process() == "hadElastic"){
192  if(partHist->PdgCode() == 2212){
193  elasticProton += flshits[idx_fls].GetEdep();
194  }
195  }
196 
197  }
198  if(it != processMap.end()){
199  processMap[partHist->Process()] += flshits[idx_fls].GetEdep();
200  }
201  else{
202  processMap["other"] += flshits[idx_fls].GetEdep();
203 
204  }
205 
206  if(tempBestE > bestE){
207  bestE = tempBestE;
208  bestFlsIdx = tempFlsIdx;
209  bestHitIdx = idx;
210  }
211  }//for flsh hits
212  } // for hits
213 
214 
215  truth.primNeutronProcessE.push_back(inelasticProton);
216  truth.primNeutronProcessE.push_back(inelasticPhoton);
217  truth.primNeutronProcessE.push_back(elasticProton);
218 
219  if (bestHitIdx == -1) {
220  // If there is no best hit index, then who knows what happened in the event...
221  truth.processMax = kModeUnknown;
222  truth.processParticleE = -5;
223  } else {
224  // Should in principle always go down this loop, but above is for protection.
225  const std::vector<sim::FLSHit> flshits = bt->HitToFLSHit (hits[bestHitIdx]);
226  int flshit_trackid = flshits[bestFlsIdx].GetTrackID();
227  const sim::ParticleHistory particleHist(&nav, flshit_trackid);
228  const sim::Particle *partHist = particleHist[particleHist.size()-1];
229 
230 
231  std::string processMax = "Unknown";
232  double currentMax = 0;
233 
234  for(auto its = processMap.cbegin(); its != processMap.cend(); ++its){
235  if(its->second > currentMax){
236  currentMax = its->second;
237  processMax = its->first;
238  }
239  }
240  if ( processMax == "Primary" ) { truth.processMax = kPrimary; }
241  else if ( processMax == "hadElastic" ) { truth.processMax = kHadElastic; }
242  else if ( processMax == "Decay" ) { truth.processMax = kDecay; }
243  else if ( processMax == "pi-Inelastic" ) { truth.processMax = kPionInelastic; }
244  else if ( processMax == "protonInleastic" ) { truth.processMax = kProtonInelastic; }
245  else if ( processMax == "neutronInelastic") { truth.processMax = kNeutronInelastic; }
246  else if ( processMax == "other" ) { truth.processMax = kOther; }
247  else { truth.processMax = kModeUnknown; }
248  truth.processParticleE = partHist->E();
249  } // if bestHitIdx == -1
250 
251  }//if puff
252 
253  return truth;
254 }
255 
256 
257  //.....................................................................
259  std::vector<cheat::TrackIDE>& allTracks,
260  std::vector<cheat::TrackIDE>& sliceTracks,
261  std::vector<cheat::TrackIDE>& allTracksBirks,
262  std::vector<cheat::TrackIDE>& sliceTracksBirks,
263  std::vector<SRTrueParticle>* vec,
264  const std::vector<sim::TrueEnergy>& TrueEnergies)
265  {
267  int trackId = part.TrackId();
268  // Find visible energy
269  double visE = 0;
270  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
271  it != allTracks.end(); ++it)
272  {
273  if(trackId == it->trackID)
274  {
275  visE += it->energy;
276  break;
277  }
278 
279  }
280  // Find visible energy in slice
281  double inSliceVisE = 0;
282  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracks.begin();
283  it != sliceTracks.end(); ++it)
284  {
285  if(trackId == it->trackID)
286  {
287  inSliceVisE += it->energy;
288  break;
289  }
290  }
291  //now calculate visible energies with birks suppression
292  // Find visible energy
293  double visEBirks = 0;
294  for (std::vector<cheat::TrackIDE>::iterator it = allTracksBirks.begin();
295  it != allTracksBirks.end(); ++it)
296  {
297  if(trackId == it->trackID)
298  {
299  visEBirks += it->energy;
300  break;
301  }
302 
303  }
304  // Find visible energy in slice
305  double inSliceVisEBirks = 0;
306  for (std::vector<cheat::TrackIDE>::iterator it = sliceTracksBirks.begin();
307  it != sliceTracksBirks.end(); ++it)
308  {
309  if(trackId == it->trackID)
310  {
311  inSliceVisEBirks += it->energy;
312  break;
313  }
314  }
315 
316 
317  double daughterVisE = 0;
318  double daughterInSliceVisE = 0;
319  double daughterVisEBirks = 0;
320  double daughterInSliceVisEBirks = 0;
321 
322 
323 
324  for (int i=0; i<part.NumberDaughters();++i){
325  const sim::Particle * p = bt->TrackIDToParticle(part.Daughter(i));
326 
327  daughterVisE += FindDaughterVisE(*p,allTracks);
328  daughterInSliceVisE += FindDaughterVisE(*p,sliceTracks);
329  daughterVisEBirks += FindDaughterVisE(*p,allTracksBirks);
330  daughterInSliceVisEBirks += FindDaughterVisE(*p,sliceTracksBirks);
331  }
332 
333 
334  // Find the total escaping energy in slice.
335  float TotEscE = -5;
336  float EnterEn = -5;
337  for(const sim::TrueEnergy& trueE: TrueEnergies){
338  if ( trackId == trueE.TrackId() ) {
339  EnterEn = trueE.EnteringEnergy();
340  TotEscE = trueE.TotalEscEnergy();
341  break;
342  }
343  }
344 
345  // Make a new SRTrueParticle and get an accessor
346  vec->push_back(SRTrueParticle());
347  SRTrueParticle& srPart = vec->back();
348  // Set fields
349  srPart.pdg = part.PdgCode();
350  srPart.visE = visE;
351  srPart.visEinslc = inSliceVisE;
352  srPart.daughterVisE = daughterVisE;
353  srPart.daughterVisEinslc = daughterInSliceVisE;
354  srPart.visEBirks = visEBirks;
355  srPart.visEinslcBirks = inSliceVisEBirks;
356  srPart.daughterVisEBirks = daughterVisEBirks;
357  srPart.daughterVisEinslcBirks = daughterInSliceVisEBirks;
358  srPart.totEscE = TotEscE;
359  srPart.enteringE = EnterEn;
360  srPart.p = part.Momentum();
361  srPart.trkID = trackId;
362 
363  float maxProInel = 0;
364  float maxProInelS =0;
365  float maxProElS =0;
366  float maxProEl =0;
367  float maxPhoInel = 0;
368  float maxPhoInelS = 0;
369 
370  float maxProInelTrueE = 0;
371  float maxProElTrueE = 0;
372  float maxPhoInelTrueE =0;
373 
374  double inelasticProtonSum = 0;
375  double inelasticPhotonSum = 0 ;
376  double elasticProtonSum = 0;
377  double inelasticProtonInslcSum = 0;
378  double inelasticPhotonInslcSum = 0 ;
379  double elasticProtonInslcSum = 0;
380 
381 
382  double inelasticProtonMax = 0;
383  double inelasticPhotonMax = 0 ;
384  double elasticProtonMax = 0;
385  double inelasticProtonInslcMax = 0;
386  double inelasticPhotonInslcMax = 0 ;
387  double elasticProtonInslcMax = 0;
388 
389 
390  for (int i=0; i<part.NumberDaughters();++i){
391  const sim::Particle * p = bt->TrackIDToParticle(part.Daughter(i));
392  srPart.daughterlist.push_back(bt->TrackIDToParticle(part.Daughter(i))->PdgCode());
393  srPart.daughterEnergies.push_back(p->E());
394 
395 
396  if(p->Process() == "neutronInelastic"){
397  if(p->PdgCode() ==2212){
398 
399  inelasticProtonSum += FindDaughterVisENonRecur(*p,allTracks);
400  inelasticProtonInslcSum += FindDaughterVisENonRecur(*p,sliceTracks);
401 
402  if( FindDaughterVisENonRecur(*p,allTracks) > maxProInel){
403  inelasticProtonMax = FindDaughterVisENonRecur(*p,allTracks);
404  maxProInel = inelasticProtonMax;
405 
406  }
407  if( FindDaughterVisENonRecur(*p,sliceTracks) > maxProInelS){
408  inelasticProtonInslcMax = FindDaughterVisENonRecur(*p,sliceTracks);
409  maxProInelS = inelasticProtonInslcMax;
410  maxProInelTrueE = p->E();
411  }
412 
413 
414  }
415  if(p->PdgCode() ==22){
416 
417  inelasticPhotonSum += FindDaughterVisENonRecur(*p,allTracks);
418  inelasticPhotonInslcSum += FindDaughterVisENonRecur(*p,sliceTracks);
419 
420  if( FindDaughterVisENonRecur(*p,allTracks) > maxPhoInel){
421  inelasticPhotonMax = FindDaughterVisENonRecur(*p,allTracks);
422  maxPhoInel = inelasticPhotonMax;
423  }
424  if( FindDaughterVisENonRecur(*p,sliceTracks) > maxPhoInelS){
425  inelasticPhotonInslcMax = FindDaughterVisENonRecur(*p,sliceTracks);
426  maxPhoInelS = inelasticPhotonInslcMax;
427  maxPhoInelTrueE = p->E();
428  }
429 
430 
431  }
432 
433  }
434  if(p->Process() == "hadElastic"){
435 
436  if(p->PdgCode() == 2212){
437 
438  elasticProtonSum += FindDaughterVisENonRecur(*p,allTracks);
439  elasticProtonInslcSum += FindDaughterVisENonRecur(*p,sliceTracks);
440 
441  if( FindDaughterVisENonRecur(*p,sliceTracks) > maxProElS){
442  elasticProtonInslcMax = FindDaughterVisENonRecur(*p,sliceTracks);
443  maxProElS = elasticProtonInslcMax;
444  maxProElTrueE = p->E();
445  }
446  if( FindDaughterVisENonRecur(*p,allTracks) > maxProEl){
447  elasticProtonMax = FindDaughterVisENonRecur(*p,allTracks);
448  maxProEl = elasticProtonMax;
449  }
450 
451  }
452 
453  }
454 
455  }
456 
457  srPart.inelasticProtonSumVisE = inelasticProtonSum;
458  srPart.inelasticPhotonSumVisE = inelasticPhotonSum;
459  srPart.elasticProtonSumVisE = elasticProtonSum;
460 
461 
462  srPart.inelasticProtonSumVisEinslc = inelasticProtonInslcSum;
463  srPart.inelasticPhotonSumVisEinslc = inelasticPhotonInslcSum;
464  srPart.elasticProtonSumVisEinslc = elasticProtonInslcSum;
465 
466 
467  srPart.inelasticProtonMaxVisE = inelasticProtonMax;
468  srPart.inelasticPhotonMaxVisE = inelasticPhotonMax;
469  srPart.elasticProtonMaxVisE = elasticProtonMax;
470 
471 
472  srPart.inelasticProtonMaxVisEinslc = inelasticProtonInslcMax;
473  srPart.inelasticPhotonMaxVisEinslc = inelasticPhotonInslcMax;
474  srPart.elasticProtonMaxVisEinslc = elasticProtonInslcMax;
475 
476 
477  srPart.maxElasticProtonTrueE = maxProElTrueE;
478  srPart.maxInelasticProtonTrueE = maxProInelTrueE;
479  srPart.maxInelasticPhotonTrueE = maxPhoInelTrueE;
480 
481 
482 
483 
484  }
485 
486  //.....................................................................
488  std::vector<cheat::TrackIDE>& tracks)
489  {
491  int trackId = part.TrackId();
492  // Find visible energy
493  double visE = 0;
494  for (std::vector<cheat::TrackIDE>::iterator it = tracks.begin();
495  it != tracks.end(); ++it)
496  {
497  if(trackId == it->trackID)
498  {
499  visE += it->energy;
500  break;
501  }
502  }
503  for (int i=0; i<part.NumberDaughters();++i){
504  const sim::Particle * p = bt->TrackIDToParticle(part.Daughter(i));
505  visE += FindDaughterVisE(*p,tracks);
506  }
507 
508  return visE;
509  }
510  //.....................................................................
511 
513  std::vector<cheat::TrackIDE>& tracks)
514  {
516  int trackId = part.TrackId();
517  // Find visible energy
518  double visE = 0;
519  for (std::vector<cheat::TrackIDE>::iterator it = tracks.begin();
520  it != tracks.end(); ++it)
521  {
522  if(trackId == it->trackID)
523  {
524  visE += it->energy;
525  break;
526  }
527  }
528 
529 
530  return visE;
531  }
532 
533 
534  //.....................................................................
536  std::vector<cheat::TrackIDE>& allTracks,
537  std::vector<cheat::TrackIDE>& sliceTracks,
538  std::vector<cheat::TrackIDE>& allTracksBirks,
539  std::vector<cheat::TrackIDE>& sliceTracksBirks,
540  std::vector<SRCosmic>* vec)
541  {
544 
545  // Get information about true cosmic
546  vec->push_back(SRCosmic());
547  SRCosmic& srTruth = vec->back();
548 
549  const art::Ptr<simb::MCTruth> &tru = effPur->neutrinoInt;
550 
551  // BUT GetParticle(int i) requires looping over list of particles in
552  // single MCTruth object. So we have made an MCTruth object for each
553  // "cosmic" and GetParticle(0) corresponds to the only particle then
554  srTruth.pdg = tru->GetParticle(0).PdgCode();
555 
556  srTruth.E = tru->GetParticle(0).Momentum().E();
557  srTruth.visE = effPur->energyTotal;
558  srTruth.visEinslc = effPur->energySlice;
559  srTruth.eff = effPur->efficiency;
560  srTruth.pur = effPur->purity;
561  srTruth.nhitslc = effPur->nSliceHits;
562  srTruth.nhittot = effPur->nTotalHits;
563 
564  srTruth.time = tru->GetParticle(0).Position().T();
565  srTruth.vtx = tru->GetParticle(0).Position().Vect();
566  srTruth.p = tru->GetParticle(0).Momentum();
567 
568 
569  srTruth.azimuth = atan2(tru->GetParticle(0).Momentum().Pz(),
570  tru->GetParticle(0).Momentum().Px());
571  srTruth.zenith = -(tru->GetParticle(0).Momentum().Py())/tru->GetParticle(0).Momentum().P();
572 
573  std::vector<const sim::Particle*> particles = bt->MCTruthToParticles(tru);
574 
575 
576  const simb::MCTrajectory& traj = particles[0]->Trajectory();
577 
578 
579  bool entered = 0;
580  bool exited = 0;
581 
582  for(size_t pt = 0; pt < traj.size(); ++pt){
583  TVector3 point = traj.Position(pt).Vect();
584 
585  // Check to see if trajectory point is in the detector
586  bool inDet = fabs(point.X()) <= geom->DetHalfWidth() // Check X inside detector
587  && fabs(point.Y()) <= geom->DetHalfHeight() // then check Y
588  && point.Z() >= 0 && point.Z() <= geom->DetLength(); // then check Z
589 
590  if(inDet && !entered){
591  // Then we just entered the detector
592  entered = true;
593 
594  if(pt > 0){
595  // Put xyz from last point into an array for geo function call
596  TVector3 lastPoint = traj.Position(pt - 1).Vect();
597  double pointArray[3] = {point.X(), point.Y(), point.Z()};
598 
599  // Put xyz from direction vector into an array.
600  TVector3 dir = (lastPoint - point).Unit();
601  double dirArray[3] = {dir.X(), dir.Y(), dir.Z()};
602 
603  double enterArray[3] = {0,0,0};
604 
605  geo::ProjectToBoxEdge(pointArray, dirArray,
606  - geom->DetHalfWidth() , geom->DetHalfWidth(),
607  - geom->DetHalfHeight() , geom->DetHalfHeight(),
608  0 , geom->DetLength(),
609  enterArray);
610 
611  srTruth.enter.SetXYZ(enterArray[0], enterArray[1], enterArray[2]);
612  // No good way to interpolate between momentum at pt-1 and pt.
613  srTruth.penter = traj.Momentum(pt-1);
614  }
615  else{ // If this is the case, we started inside the detector.
616  srTruth.enter = point;
617  srTruth.penter = traj.Momentum(pt);
618  }
619  }
620 
621  if(!inDet && entered && !exited){
622  // Then we just exited the detector
623  exited= true;
624 
625  // Put xyz from last point into an array for geo function call
626  TVector3 lastPoint = traj.Position(pt - 1).Vect();
627  double lastPointArray[3] = {lastPoint.X(), lastPoint.Y(), lastPoint.Z()};
628 
629  // Put xyz from direction vector into an array.
630  TVector3 dir = (point - lastPoint).Unit();
631  double dirArray[3] = {dir.X(), dir.Y(), dir.Z()};
632 
633  double exitArray[3] = {0,0,0};
634 
635  geo::ProjectToBoxEdge(lastPointArray, dirArray,
636  - geom->DetHalfWidth() , geom->DetHalfWidth(),
637  - geom->DetHalfHeight() , geom->DetHalfHeight(),
638  0 , geom->DetLength(),
639  exitArray);
640 
641  srTruth.exit.SetXYZ(exitArray[0], exitArray[1], exitArray[2]);
642  }
643  }
644  srTruth.stop = traj.Position(traj.size() - 1).Vect();
645 
646  // Make sure the particle actually entered the detector
647  if(!entered){
648  // If the primary never entered the detector, set the
649  // coordinates to -5*10^9
650 
651  srTruth.enter.SetXYZ(-5e9, -5e9, -5e9);
652  srTruth.exit.SetXYZ(-5e9, -5e9, -5e9);
653  }
654 
655  // If the particle entered and did not exit, set the exit position to stop
656  if(entered && !exited){
657  srTruth.exit = srTruth.stop;
658  }
659  // Find the Michel electrons
660  FindAndAddMichels(particles, allTracks, &srTruth.michel);
661  }
662 
663  //......................................................................
664  void FindAndAddMichels(std::vector<const sim::Particle*> particles,
665  std::vector<cheat::TrackIDE>& allTracks,
666  std::vector<SRTrueMichelE>* michelVec)
667  {
669 
670  for(const sim::Particle* part: particles){
671  // Check to see if we're dealing with an electron, continue if not.
672  int pdg = part->PdgCode();
673  if (abs(pdg) != 11) continue;
674 
675  // Check to see if mother is a muon, continue if not
676  int motherId = part->Mother();
677  if(motherId == 0) continue;
678  const sim::Particle* mother = bt->TrackIDToParticle( motherId );
679  int pdgMother = 0;
680  if(mother) pdgMother = mother->PdgCode();
681  if(abs(pdgMother) != 13) continue;
682 
683  // By definition, Michel electrons only arise in muon decay.
684  // However, in GEANT4, negative muons only produce electrons via "Decay"
685  // if they decay in flight; otherwise, both decays and capture
686  // are handled via the "capture-at-rest" module.
687  // Since actual muon capture is a process containing no electrons,
688  // any electrons that arose from "muMinusCaptureAtRest" are Michels.
689  if(part->Process() != "Decay" && part->Process() != "muMinusCaptureAtRest") continue;
690 
691  // Add particle if it passed all criteria
692  AddTrueMichelEToVec(*part, *bt->TrackIDToParticle( motherId ),
693  allTracks, michelVec);
694  }
695  }
696 
697  //......................................................................
698  void AddTrueMichelEToVec(const sim::Particle& michel,
699  const sim::Particle& motherMuon,
700  std::vector<cheat::TrackIDE>& allTracks,
701  std::vector<SRTrueMichelE>* vec)
702  {
703  vec->push_back(SRTrueMichelE());
704  SRTrueMichelE& srMichel = vec->back();
705 
706  // Find the visible energy from the FLS hits
707  float visE = 0;
708  for (std::vector<cheat::TrackIDE>::iterator it = allTracks.begin();
709  it != allTracks.end(); ++it)
710  {
711  if(michel.TrackId() == it->trackID)
712  {
713  visE = it->energy;
714  break;
715  }
716  }
717 
718  srMichel.E = michel.E();
719  srMichel.visE = visE;
720  srMichel.time = michel.T();
721  srMichel.mustop = motherMuon.EndPosition();
722  srMichel.p = michel.Momentum();
723  }
724 
725  //......................................................................
726  void AddPreFSI( const art::Ptr<simb::MCTruth>& truth, SRNeutrino& nu )
727  {
728  // Add information about "pre-FSI hadrons", which are recoil hadrons inside the nucleus after
729  // resonance decays and hadronization and before FSI. This excludes hadrons from hadronization
730  // and resonance decays occuring outside the nucleus, as well as hadrons from coherent scattering
731  // and scattering on hydrogen.
732 
733  // Also included is the mapping between pre-FSI hadrons and their final state descendants.
734  // Final state particle information in the CAFs comes from sim::Particle objects from the
735  // BackTracker service and is stored in SRTrueParticle objects in SRNeutrino. Information of
736  // final state descendants of pre-FSI hadrons is stored in simb::MCParticle objects from
737  // simb::MCTruth. This situation requires matching the pre-FSI hadron final state descendant
738  // simb::MCParticle objects to the corresponding sim::Particle or SRTrueParticle objects by
739  // pdg code and 4-momentum.
740 
741  // First, retrieve the pre-FSI hadrons and their final state descendants
742 
743  std::map<const simb::MCParticle*,std::vector<const simb::MCParticle*>> prefsi_final_descendants;
744 
745  for ( int iprefsi = 0; iprefsi < truth->NParticles(); ++iprefsi )
746  {
747  const simb::MCParticle* prefsi_part = &truth->GetParticle( iprefsi );
748  // Require hadron undergoing intranuclear transport (status code = 14)
749  if ( prefsi_part->StatusCode() != 14 ) continue;
750  // Require hadron not a product of FSI
751  if ( truth->GetParticle( prefsi_part->Mother() ).StatusCode() == 14 ) continue;
752  // Find final state descendants
753  std::vector<const simb::MCParticle*> final_descendants;
754  for ( int ifinal = 0; ifinal < truth->NParticles(); ++ifinal )
755  {
756  const simb::MCParticle* final_part = &truth->GetParticle( ifinal );
757  // Require final state particle (particle (status code = 1)
758  if ( final_part->StatusCode() != 1 ) continue;
759  // Exclude at-rest final state particles produced by GENIE presumably for baryon accounting.
760  // Iddentical paricles with zero momentum are not unique and cannot be matched as described above.
761  if ( final_part->Momentum().Vect().Mag() == 0 ) continue;
762  int iancestor = final_part->Mother();
763  int search_count = 0;
764  while ( iancestor != iprefsi && iancestor > 0 && search_count < truth->NParticles() )
765  {
766  iancestor = truth->GetParticle( iancestor ).Mother();
767  search_count++;
768  }
769  if ( iancestor == iprefsi ) final_descendants.push_back( final_part );
770  }
771  prefsi_final_descendants[ prefsi_part ] = final_descendants;
772  }
773 
774  // Determine mapping between pre-FSI hadrons and their final state descendants in the nu.prim
775  // vector and add the pre-FSI hadron information to the record.
776 
777  for ( auto pfd : prefsi_final_descendants )
778  {
779  SRNuGenParticle nugen_part;
780  nugen_part.pdg = pfd.first->PdgCode();
781  nugen_part.status = pfd.first->StatusCode();
782  nugen_part.p = pfd.first->Momentum();
783 
784  const int prefsiID = (int)nu.prefsi.size();
785 
786  for ( unsigned int iprim = 0; iprim < nu.prim.size(); ++iprim )
787  {
788  for ( auto final_part : pfd.second )
789  {
790  double max_diff = 1e-4;
791  if ( final_part->PdgCode() == nu.prim[ iprim ].pdg &&
792  fabs( final_part->Px() - nu.prim[ iprim ].p.px ) < max_diff &&
793  fabs( final_part->Py() - nu.prim[ iprim ].p.py ) < max_diff &&
794  fabs( final_part->Pz() - nu.prim[ iprim ].p.pz ) < max_diff &&
795  fabs( final_part->E() - nu.prim[ iprim ].p.E ) < max_diff )
796  {
797  nugen_part.primList.push_back( iprim );
798  nu.prim[ iprim ].prefsiID = prefsiID;
799  }
800  }
801  }
802 
803  nu.prefsi.push_back( nugen_part );
804 
805  int nfinal = pfd.second.size();
806  int nprim = nu.prefsi[ prefsiID ].primList.size();
807  if ( nfinal != nprim )
808  {
809  std::cout << Form( "Mismatch: pre-FSI N final state descendants = %d, N primary matches = %d", nfinal, nprim ) << std::endl;
810  //if ( nfinal != nprim ) LOG_WARNING( "FillTruth" ) << Form( "Mismatch: pre-FSI N final state descendants = %d, N primary matches = %d", nfinal, nprim );
811  }
812  }
813 
814  }
815 
816  //......................................................................
818  {
819  // TODO: I think this is actually 809.something? And maybe computing
820  // it on an event-by-event basis could be interesting to someone?
821  if(det == novadaq::cnv::kFARDET) return 810;
822 
823  // TODO: ideally this would be in some more appropriate location.
824  // Perhaps something like Geometry::BeamCoordsToDetectorCoords()?
825  const double rotation[] = {
826  9.9990e-1, -8.2300e-4, -1.4088e-2,
827  3.0533e-6, 9.9831e-1, -5.8103e-2,
828  1.4112e-2, 5.8097e-2, 9.8821e-1
829  };
830 
831  const double translation[] = {200.5991,
832  6000.9565,
833  -99900.2012};
834 
835  const double origin[] = {nu.beam.v.X(),
836  nu.beam.v.Y(),
837  nu.beam.v.Z()};
838 
839  const TMatrixD rot(3, 3, rotation);
840  const TVectorD trans(3, translation);
841  const TVectorD orig(3, origin);
842 
843  const TVectorD hadVtx = rot*orig+trans;
844 
845  const TVector3 hadVtx3(hadVtx[0], hadVtx[1], hadVtx[2]);
846 
847  return (nu.vtx - hadVtx3).Mag()/100000;
848  }
849 
850  //......................................................................
851  std::vector<SRGenieWeights>
853  {
854  const unsigned int N = table.NShifts();
855 
856  std::vector<caf::SRGenieWeights> ret(N);
857 
858  for(unsigned int i = 0; i < N; ++i){
859  ret[i].minus2sigma = table.Minus2Sigma(i);
860  ret[i].minus1sigma = table.Minus1Sigma(i);
861  ret[i].plus1sigma = table.Plus1Sigma(i);
862  ret[i].plus2sigma = table.Plus2Sigma(i);
863  }
864 
865  return ret;
866  }
867 
868 
869  //......................................................................
871  {
872  SRFluxWeights frew;
873  frew.cv = flxwgts.GetTotalCentralValue();
874  frew.vuniv = flxwgts.GetTotalMultiUniverses();
875  frew.nvuniv = flxwgts.GetNumberOfUniverses();
876  return frew;
877  }
878 
879  //......................................................................
881  {
882  SRGeant4Weights g4rw;
883 
884  g4rw.piplus_cv = g4wgts.Piplus_CV();
885  g4rw.piminus_cv = g4wgts.Piminus_CV();
886  g4rw.proton_cv = g4wgts.Proton_CV();
887 
888  g4rw.piplus_univ = g4wgts.Piplus();
889  g4rw.piminus_univ = g4wgts.Piminus();
890  g4rw.proton_univ = g4wgts.Proton();
891 
892  // It's redundant filling 3 fields with the same value
893  // but CAFAna likes the size of "vec_name" to be called
894  // "nvec_name"
895  g4rw.npiplus_univ = g4wgts.nUniverses();
896  g4rw.npiminus_univ = g4wgts.nUniverses();
897  g4rw.nproton_univ = g4wgts.nUniverses();
898 
899 
900  return g4rw;
901  }
902 
903 
904  //......................................................................
906  {
907  switch (simbGeneratorEnum)
908  {
910  return kGENIE;
911 
913  return kGIBUU;
914 
915  default:
916  return kUnknownGenerator;
917 
918  }
919  }
920 
921  //......................................................................
922  std::vector<unsigned int> DecodeGeneratorVersion(const std::string &versionString)
923  {
924  std::vector<unsigned int> ret;
925  std::regex pattern("(\\d+)");
926  auto vals_begin = std::sregex_iterator(versionString.begin(), versionString.end(), pattern);
927  auto vals_end = std::sregex_iterator();
928 
929  for (auto it = vals_begin; it != vals_end; it++)
930  ret.push_back(std::stoi(it->str()));
931 
932  return ret;
933  }
934 
935 }
SRLorentzVector p
Momentum 4-vector.
Definition: SRTrueMichelE.h:27
double E(const int i=0) const
Definition: MCParticle.h:232
caf::generator_ CAFGeneratorEnum(simb::Generator_t simbGeneratorEnum)
Definition: FillTruth.cxx:905
double Plus1Sigma(int i) const
float daughterVisE
Visible Energy in detector for all daughters of this particle, all summed FLSHits that made CellHits ...
float daughterVisEinslcBirks
Visible Energy in detector for all daughters of this particle, slice summed FLSHits that made CellHit...
double TrueNeutrinoDistance(novadaq::cnv::DetId det, const SRNeutrino &nu)
Definition: FillTruth.cxx:817
SRVector3D vtx
Vertex position in detector coordinates [cm].
Definition: SRNeutrino.h:47
float maxInelasticProtonTrueE
Energy of the proton daughter going through the most energetic inelastic process. ...
double Minus2Sigma(int i) const
back track the reconstruction to the simulation
const TLorentzVector & Position(const int i=0) const
Definition: MCParticle.h:218
int PdgCode() const
Definition: MCParticle.h:211
float visE
Visible Energy in detector, all summed FLSHits that made CellHits [GeV].
Definition: SRTrueMichelE.h:22
float piplus_cv
Reweight for the piplus central value (cv)
Reweight information for geant4 systematic.
float inelasticPhotonMaxVisEinslc
Vis energy in slc coming from max inelastic process with photons in final state linked to primary...
std::vector< float > daughterVisEnergies
Energy of each particle contributing to the prong.
float FindDaughterVisE(const sim::Particle &part, std::vector< cheat::TrackIDE > &tracks)
Definition: FillTruth.cxx:487
set< int >::iterator it
SRVector3D enter
Cosmic entrance point in detector coordinates. [cm] When the primary doesn&#39;t enter the detector...
Definition: SRCosmic.h:40
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
const TLorentzVector & EndPosition() const
Definition: MCParticle.h:224
float daughterVisEinslc
Visible Energy in detector for all daughters of this particle, slice summed FLSHits that made CellHit...
SRBeam beam
Information about neutrino production.
Definition: SRNeutrino.h:93
void SetXYZ(float x, float y, float z)
Definition: SRVector3D.cxx:31
void AddPreFSI(const art::Ptr< simb::MCTruth > &truth, SRNeutrino &nu)
Definition: FillTruth.cxx:726
void AddParticleToVec(const sim::Particle &part, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, std::vector< SRTrueParticle > *vec, const std::vector< sim::TrueEnergy > &TrueEnergies)
Definition: FillTruth.cxx:258
const sim::Particle * TrackIDToMotherParticle(int const &id) const
float elasticProtonSumVisEinslc
Vis energy in slc coming from sum of elastic processes with protons in final state linked to primary...
const sim::ParticleNavigator & ParticleNavigator() const
Get a reference to the ParticleNavigator.
Definition: BackTracker.h:744
std::vector< float > piplus_univ
Reweight values for the piplus multi-universe.
SRLorentzVector p
True energy 4-vector of the best matched particle.
float elasticProtonSumVisE
Visible energy coming from sum of elastic processes with protons in final state linked to primary...
std::vector< float > proton_univ
Reweight values for the proton multi-universe.
float azimuth
Azimuth angle (w.r.t y-axis)
Definition: SRCosmic.h:38
int Mother() const
Definition: MCParticle.h:212
float inelasticPhotonSumVisEinslc
Vis energy in slc coming from sum of inelastic processes with photons in final state linked to primar...
const char * p
Definition: xmltok.h:285
std::vector< unsigned int > primList
Indices to primary (i.e. final state) descendants in SRNeutrino.
std::vector< SRTrueParticle > prim
Primary daughters, lepton comes first in vector.
Definition: SRNeutrino.h:84
float Y() const
Definition: SRVector3D.h:33
short pdg
pdg code
cheat::ParticleEffPur ClusterToParticle(const std::vector< const rb::Cluster * > &sliceList, const std::vector< const rb::CellHit * > &hits, const int &sliceIdx)
int trackID
Truth particle track ID.
Definition: BackTracker.h:66
SRLorentzVector p
True momentum [GeV].
Definition: SRCosmic.h:32
float visEinslc
Sum of FLS hits that made CellHits from this neutrino in this subevent [GeV].
Definition: SRCosmic.h:24
float Piplus_CV() const
Definition: G4WeightTable.h:27
double DetLength() const
int nTotalHits
Total number of hits the neutrino left in the detector.
Definition: BackTracker.h:54
Reweight information for flux systematic.
Definition: SRFluxWeights.h:12
void abs(TH1 *hist)
double energyTotal
Sum of FLS hits from the neutrino contributing to any hit in the event.
Definition: BackTracker.h:52
enum simb::_ev_generator Generator_t
generator used to produce event, if applicable
int StatusCode() const
Definition: MCParticle.h:210
SRLorentzVector p
True momentum [GeV].
float eff
True deposited energy collection efficiency for the best matched particle relative to the slice...
int NParticles() const
Definition: MCTruth.h:74
double purity
Putity of particle relative to track.
Definition: BackTracker.h:68
float piminus_cv
Reweight for the piminus central value (cv)
SRGeant4Weights Geant4Reweights(const g4rwgt::G4WeightTable &g4wgts)
Definition: FillTruth.cxx:880
float inelasticProtonMaxVisE
Visible energy coming from max inelastic process with protons in final state linked to primary...
std::string Process() const
Definition: MCParticle.h:214
std::vector< float > GetTotalMultiUniverses() const
Get a vector of total multiuniverse weights (every entry is the product of the HP and Beam Focusing)...
Definition: FluxWeights.cxx:58
float elasticProtonMaxVisEinslc
Vis energy in slc coming from max elastic process with protons in final state linked to primary...
Particle class.
float daughterVisEinslcBirks
Visible Energy in detector for all daughters of this particle, slice summed FLSHits that made CellHit...
float daughterVisEBirks
Visible Energy in detector for all daughters of this particle, all summed FLSHits that made CellHits ...
double Plus2Sigma(int i) const
SRFluxWeights FluxReweights(const fxwgt::FluxWeights &flxwgts)
Definition: FillTruth.cxx:870
double efficiency
Efficiency (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:48
generator_
Known generators of neutrino interactions (extend as other generators are used)
Definition: SREnums.h:28
int NumberDaughters() const
Definition: MCParticle.h:216
float FindDaughterVisENonRecur(const sim::Particle &part, std::vector< cheat::TrackIDE > &tracks)
Definition: FillTruth.cxx:512
std::vector< T > PtrVecToVec(const art::PtrVector< T > &xs)
Definition: Utils.h:10
std::vector< int > daughterlist
std::vector< float > Piplus() const
Definition: G4WeightTable.h:21
float pur
Slicer purity for this truth interaction.
Definition: SRCosmic.h:26
The truth information of reco objects within a slice.
SRVector3D vtx
Vertex position in detector coordinates [cm].
Definition: SRCosmic.h:34
float Piminus_CV() const
Definition: G4WeightTable.h:28
int TrackId() const
Definition: MCParticle.h:209
int Daughter(const int i) const
Definition: MCParticle.cxx:112
The SRNeutrino is a representation of neutrino interaction information.
Definition: SRNeutrino.h:21
float zenith
Zenith angle (w.r.t y-axis)
Definition: SRCosmic.h:39
SRVector3D stop
Cosmic end point in detector coordinates, regardless of whether it is inside or outside the detector...
Definition: SRCosmic.h:43
float E
True energy [GeV].
Definition: SRCosmic.h:22
float X() const
Definition: SRVector3D.h:32
float visEinslc
Visible Energy in detector, slice summed FLSHits that made CellHits [GeV].
SRParticleTruth FillParticleTruth(const std::vector< rb::Cluster > &sliceList, const art::PtrVector< rb::CellHit > &hits, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, int sliceIdx)
Function to fill the particle truth for a set of hits.
Definition: FillTruth.cxx:28
float daughterVisEBirks
Visible Energy in detector for all daughters of this particle, all summed FLSHits that made CellHits ...
TString part[npart]
Definition: Style.C:32
float visE
Visible Energy in detector, all summed FLSHits that made CellHits [GeV].
Far Detector at Ash River, MN.
std::vector< float > piminus_univ
Reweight values for the piminus multi-universe.
void hits()
Definition: readHits.C:15
unsigned int npiplus_univ
Number of piplus universes.
int PdgCode(Resonance_t res, int Q)
(resonance id, charge) -> PDG code
std::vector< float > vuniv
Reweight values for the multi-universe.
Definition: SRFluxWeights.h:18
float visEinslcBirks
Visible Energy in detector, slice summed FLSHits that made CellHits [GeV] with birks suppression...
void FindAndAddMichels(std::vector< const sim::Particle * > particles, std::vector< cheat::TrackIDE > &allTracks, std::vector< SRTrueMichelE > *michelVec)
Definition: FillTruth.cxx:664
std::vector< float > daughterEnergies
Vector containing energy of each daughter.
int pdg
PDG Code of the best matched truth particle.
float eff
Slicer efficiency for this truth interaction.
Definition: SRCosmic.h:25
float time
Time of first Michel electron trajectory point [GeV].
Definition: SRTrueMichelE.h:23
float daughterVisE
Visible Energy in detector for all daughters of this particle, all summed FLSHits that made CellHits ...
Store flux weigths for neutrino correction.
Definition: FluxWeights.h:15
float totEscE
The total escaping energy, from the particle and all of its daughters, using truth information [GeV]...
std::vector< float > Proton() const
Definition: G4WeightTable.h:23
std::vector< sim::FLSHit > HitToFLSHit(const rb::CellHit &hit) const
All the FLSHits that contributed to this hit, sorted from most to least light.
int motherpdg
PDG Code of the mother of the best matched truth particle.
Truth information for a Michel electron.
Definition: SRTrueMichelE.h:15
double efficiency
Efficiency of particle relative to track.
Definition: BackTracker.h:67
float enteringE
The kinetic energy the particle had when it first entered the detector, using truth information [GeV]...
double T(const int i=0) const
Definition: MCParticle.h:223
float pur
True deposited energy purity for the best matched particle.
Collect Geo headers and supply basic geometry functions.
Eigen::VectorXd vec
SRVector3D exit
Cosmic exit point in detector coordinates. [cm] When the primary doesn&#39;t enter the detector...
Definition: SRCosmic.h:41
int pdg
Pdg of particle.
Definition: BackTracker.h:69
double DetHalfHeight() const
std::vector< unsigned int > DecodeGeneratorVersion(const std::string &versionString)
Definition: FillTruth.cxx:922
float visEBirks
Visible Energy in detector, all summed FLSHits that made CellHits [GeV] with birks suppression...
TVector3 Unit() const
const TLorentzVector & Position(const size_type) const
The accessor methods described above.
float Z() const
Definition: SRVector3D.h:34
void AddTrueMichelEToVec(const sim::Particle &michel, const sim::Particle &motherMuon, std::vector< cheat::TrackIDE > &allTracks, std::vector< SRTrueMichelE > *vec)
Definition: FillTruth.cxx:698
float proton_cv
Reweight for the proton central value (cv)
double Minus1Sigma(int i) const
float maxInelasticPhotonTrueE
Energy of the photon daughter going through the most energetic inelastic process. ...
size_type size() const
Definition: PtrVector.h:308
float E
True energy of electron [GeV].
Definition: SRTrueMichelE.h:19
const sim::Particle * TrackIDToParticle(int const &id) const
Returns a pointer to the sim::Particle object corresponding to the given TrackID. ...
OStream cout
Definition: OStream.cxx:6
std::vector< float > primNeutronProcessE
const simb::MCParticle & GetParticle(int i) const
Definition: MCTruth.h:75
unsigned int nhitslc
Number of hits recorded in this slice by this truth interaction.
Definition: SRCosmic.h:27
std::vector< int > motherlist
float visE
Visible Energy in detector, all summed FLSHits that made CellHits [GeV].
float inelasticPhotonSumVisE
Visible energy coming from sum of inelastic processes with photons in final state linked to primary...
float visEinslcBirks
Visible Energy in detector, slice summed FLSHits that made CellHits [GeV] with birks suppression...
int status
Generator status code.
int nSliceHits
Number of hits from this neutrino in this slice.
Definition: BackTracker.h:53
float visEBirks
Visible Energy in detector, all summed FLSHits that made CellHits [GeV] with birks suppression...
gen_process_t processMax
The process conributing the most the prong.
double DetHalfWidth() const
float visE
Sum of FLS hits that made CellHits from this neutrino [GeV].
Definition: SRCosmic.h:23
unsigned int GetNumberOfUniverses() const
Get the total number of universes used.
Definition: FluxWeights.cxx:68
SRLorentzVector mustop
Stopping position of parent muon [cm].
Definition: SRTrueMichelE.h:26
float maxElasticProtonTrueE
Energy of the proton daughter going through the most energetic elastic process.
std::vector< SRTrueMichelE > michel
Vector of true Michel electrons.
Definition: SRCosmic.h:36
size_type size() const
Definition: MCTrajectory.h:165
float inelasticProtonSumVisEinslc
Vis energy in slc coming from sum of inelastic processes with protons in final state linked to primar...
TDirectory * dir
Definition: macro.C:5
SRNuGenParticle represents a particle within the simulation of a neutrino interaction.
Store +/-1,2sigma shifts for all GENIE reweighting systematics.
float inelasticPhotonMaxVisE
Visible energy coming from max inelastic process with photons in final state linked to primary...
unsigned int nvuniv
Number of universes.
Definition: SRFluxWeights.h:19
std::vector< int > daughterlist
float Proton_CV() const
Definition: G4WeightTable.h:29
double energySlice
Sum of FLS hits from the neutrino contributing to hits included in the slice.
Definition: BackTracker.h:51
const TLorentzVector & Momentum(const int i=0) const
Definition: MCParticle.h:219
void geom(int which=0)
Definition: geom.C:163
float elasticProtonMaxVisE
Visible energy coming from max elastic process with protons in final state linked to primary...
uint nUniverses() const
Definition: G4WeightTable.h:18
short pdg
pdg code
Definition: SRCosmic.h:19
std::vector< SRNuGenParticle > prefsi
Pre-FSI hadrons.
Definition: SRNeutrino.h:85
float processParticleE
Energy of the particle causing the process that contirbuted the most.
SRVector3D v
"vertex". Position of hadron/muon decay
Definition: SRBeam.h:46
float visEinslc
Visible Energy in detector, slice summed FLSHits that made CellHits [GeV].
float inelasticProtonMaxVisEinslc
Vis energy in slc coming from max inelastic process with protons in final state linked to primary...
float cv
Reweight for the central value (cv)
Definition: SRFluxWeights.h:17
void ProjectToBoxEdge(const double xyz[], const double dxyz[], double xlo, double xhi, double ylo, double yhi, double zlo, double zhi, double xyzout[])
Project along a direction from a particular starting point to the edge of a box.
Definition: Geo.cxx:38
unsigned int npiminus_univ
Number of piminus universes.
Truth information for cosmic rays.
Definition: SRCosmic.h:15
This module creates Common Analysis Files.
Definition: FileReducer.h:10
int trkID
GEANT trackId for particle.
Float_t e
Definition: plot.C:35
SRLorentzVector penter
True momentum when entering the detector [GeV].
Definition: SRCosmic.h:33
void AddCosmicTruthToVec(const cheat::NeutrinoEffPur *effPur, std::vector< cheat::TrackIDE > &allTracks, std::vector< cheat::TrackIDE > &sliceTracks, std::vector< cheat::TrackIDE > &allTracksBirks, std::vector< cheat::TrackIDE > &sliceTracksBirks, std::vector< SRCosmic > *vec)
Definition: FillTruth.cxx:535
unsigned int nproton_univ
Number of proton universes.
float daughterVisEinslc
Visible Energy in detector for all daughters of this particle, slice summed FLSHits that made CellHit...
SRLorentzVector motherp
True energy 4-vector of the mother particle.
int trkID
GEANT trackID.
const TLorentzVector & Momentum(const size_type) const
Encapsulate the geometry of one entire detector (near, far, ndos)
std::vector< float > Piminus() const
Definition: G4WeightTable.h:22
float time
interaction time.
Definition: SRCosmic.h:31
unsigned int nhittot
Total number of hits recorded for this truth interaction.
Definition: SRCosmic.h:28
int NParticles(int pdg_code, const TClonesArray *const particle_list)
unsigned int NShifts() const
float inelasticProtonSumVisE
Visible energy coming from sum of inelastic processes with protons in final state linked to primary...
TVector3 Vect() const
std::vector< SRGenieWeights > GenieReweightTable(const rwgt::GENIEReweightTable &table)
Definition: FillTruth.cxx:852
art::Ptr< simb::MCTruth > neutrinoInt
Truth information about neutrino interaction.
Definition: BackTracker.h:47
float GetTotalCentralValue() const
Get the total central value correction calculated as the product of the HP and Beam Focusing componen...
Definition: FluxWeights.cxx:46
SRLorentzVector p
Momentum 4-vector.
T atan2(T number)
Definition: d0nt_math.hpp:72
std::vector< float > primNeutronE
Energy of the primary neutron that is linked to the prong, if one exists.
std::vector< const sim::Particle * > MCTruthToParticles(art::Ptr< simb::MCTruth > const &mct) const
double purity
Purity (based on FLS energy) of neutrino interaction relative to slice.
Definition: BackTracker.h:49
enum BeamMode string