LArSoft  v06_85_00
Liquid Argon Software toolkit - http://larsoft.org/
MCTrackCollectionAnaAlg.cxx
Go to the documentation of this file.
1 
11 #include "TTree.h"
12 #include <numeric>
13 
15 
17 {
18  fTree = tree;
19  fFillTree = fill;
20 
21  fTree->Branch("mct_run",&fRun,"mct_run/i");
22  fTree->Branch("mct_event",&fEvent,"mct_event/i");
23 
24  fTree->Branch("mct_nmctracks",&fNMCTracks,"mct_nmctracks/i");
25  fTree->Branch("mct_dp",&fDParticle,"mct_dp/i");
26  fTree->Branch("mct_dpfrac",&fDParticleFraction,"mct_dpfrac/F");
27 
28  fTree->Branch("mct_dporigin",&fDParticleOrigin,"mct_dporigin/I");
29  fTree->Branch("mct_dppdg",&fDParticlePdgCode,"mct_dppdg/I");
30  fTree->Branch("mct_dptrackid",&fDParticleTrackId,"mct_dptrackid/i");
31  fTree->Branch("mct_dpmpdg",&fDParticleMotherPdgCode,"mct_dpmpdg/I");
32  fTree->Branch("mct_dpmtrackid",&fDParticleMotherTrackId,"mct_dpmtrackid/i");
33  fTree->Branch("mct_dpapdg",&fDParticleAncestorPdgCode,"mct_dpapdg/I");
34  fTree->Branch("mct_dpatrackid",&fDParticleAncestorTrackId,"mct_dpatrackid/i");
35 
36  fTree->Branch("mct_dpstartx",&fDParticleStartX,"mct_dpstartx/F");
37  fTree->Branch("mct_dpstarty",&fDParticleStartY,"mct_dpstarty/F");
38  fTree->Branch("mct_dpstartz",&fDParticleStartZ,"mct_dpstartz/F");
39  fTree->Branch("mct_dpstarte",&fDParticleStartE,"mct_dpstarte/F");
40  fTree->Branch("mct_dpendx",&fDParticleEndX,"mct_dpendx/F");
41  fTree->Branch("mct_dpendy",&fDParticleEndY,"mct_dpendy/F");
42  fTree->Branch("mct_dpendz",&fDParticleEndZ,"mct_dpendz/F");
43  fTree->Branch("mct_dpende",&fDParticleEndE,"mct_dpende/F");
44 
45  fTree->Branch("mct_coly",&fCollectionY,"mct_coly/F");
46  fTree->Branch("mct_colz",&fCollectionZ,"mct_colz/F");
47  fTree->Branch("mct_colx",&fCollectionX,"mct_colx/F");
48  fTree->Branch("mct_colrmsy",&fCollectionRMSY,"mct_colrmsy/F");
49  fTree->Branch("mct_colrmsz",&fCollectionRMSZ,"mct_colrmsz/F");
50  fTree->Branch("mct_colrmsx",&fCollectionRMSX,"mct_colrmsx/F");
51  fTree->Branch("mct_colE",&fCollectionEnergy,"mct_colE/F");
52 
53  fTree->Branch("mct_miny",&fMinY,"mct_miny/F");
54  fTree->Branch("mct_maxy",&fMaxY,"mct_maxy/F");
55  fTree->Branch("mct_minz",&fMinZ,"mct_minz/F");
56  fTree->Branch("mct_maxz",&fMaxZ,"mct_maxz/F");
57  fTree->Branch("mct_minx",&fMinX,"mct_minx/F");
58  fTree->Branch("mct_maxx",&fMaxX,"mct_maxx/F");
59 
60 }
61 
62 void sim::MCTrackCollectionAnaAlg::FillTree(unsigned int run, unsigned int event,
63  const std::vector<sim::MCTrack>& mctVec)
64 {
65  fRun = run;
66  fEvent = event;
67 
68  fNMCTracks = mctVec.size();
70 
71  fMinY = 99999;
72  fMinZ = 99999;
73  fMinX = 99999;
74  fMaxY = -99999;
75  fMaxZ = -99999;
76  fMaxX = -99999;
77 
78  if(mctVec.size()==0){
79  fTree->Fill();
80  return;
81  }
82 
83 
84  std::vector<double> y_vals;
85  std::vector<double> z_vals;
86  std::vector<double> x_vals;
87  std::vector<double> stepL_vals;
88 
89  std::vector<double> mct_length(mctVec.size(),0);
90  //int d_index = -1; // unused
91  for(size_t i_p=0; i_p<mctVec.size(); i_p++){
92 
93 
94  fCollectionEnergy += mctVec[i_p].Start().E();
95 
96  //y_vals.reserve( y_vals.size() + mctVec[i_p].size()-1 );
97  //z_vals.reserve( z_vals.size() + mctVec[i_p].size()-1 );
98  //x_vals.reserve( x_vals.size() + mctVec[i_p].size()-1 );
99 
100  for(size_t i_s=1; i_s < mctVec[i_p].size(); i_s++){
101 
102  TVector3 const& vec1 = mctVec[i_p][i_s-1].Position().Vect();
103  TVector3 const& vec2 = mctVec[i_p][i_s].Position().Vect();
104  double stepL = (vec2-vec1).Mag();
105 
106  double thisy = 0.5*(vec1.Y() + vec2.Y());
107  double thisz = 0.5*(vec1.Z() + vec2.Z());
108  double thisx = 0.5*(vec1.X() + vec2.X());
109 
110  y_vals.push_back(thisy);
111  z_vals.push_back(thisz);
112  x_vals.push_back(thisx);
113  stepL_vals.push_back(stepL);
114 
115  mct_length[i_p] += stepL;
116 
117  if(thisy > fMaxY) fMaxY = thisy;
118  if(thisy < fMinY) fMinY = thisy;
119  if(thisz > fMaxZ) fMaxZ = thisz;
120  if(thisz < fMinZ) fMinZ = thisz;
121  if(thisx > fMaxX) fMaxX = thisx;
122  if(thisx < fMinX) fMinX = thisx;
123 
124  }//end loop over steps
125  }//end loop over tracks
126 
127  double totalL = std::accumulate(mct_length.begin(),mct_length.end(),0.0);
128 
129  fDParticle = std::distance(mct_length.begin(),std::max_element(mct_length.begin(),mct_length.end()));
130  fDParticleFraction = mct_length.at(fDParticle) / totalL;
132 
133  double sumy=0,sumz=0,sumx=0;
134  for(size_t i_step=0; i_step<stepL_vals.size(); i_step++){
135  sumy += stepL_vals[i_step]*y_vals[i_step];
136  sumz += stepL_vals[i_step]*z_vals[i_step];
137  sumx += stepL_vals[i_step]*x_vals[i_step];
138  }
139 
140  fCollectionY = sumy/totalL;
141  fCollectionZ = sumz/totalL;
142  fCollectionX = sumx/totalL;
143 
144  double sumy2=0,sumz2=0,sumx2=0;
145  for(size_t i_step=0; i_step<stepL_vals.size(); i_step++){
146  sumy2 += stepL_vals[i_step]*(y_vals[i_step]-fCollectionY)*(y_vals[i_step]-fCollectionY);
147  sumz2 += stepL_vals[i_step]*(z_vals[i_step]-fCollectionZ)*(z_vals[i_step]-fCollectionZ);
148  sumx2 += stepL_vals[i_step]*(x_vals[i_step]-fCollectionX)*(x_vals[i_step]-fCollectionX);
149  }
150 
151  fCollectionRMSY = std::sqrt(sumy2/totalL);
152  fCollectionRMSZ = std::sqrt(sumz2/totalL);
153  fCollectionRMSX = std::sqrt(sumx2/totalL);
154 
155  if(fFillTree) fTree->Fill();
156 }
157 
159 {
160  size_t nsteps = mctrack.size();
161 
162  fDParticleOrigin = mctrack.Origin();
163  fDParticlePdgCode = mctrack.PdgCode();
164  fDParticleTrackId = mctrack.TrackID();
165  fDParticleStartY = mctrack[0].Y();
166  fDParticleStartZ = mctrack[0].Z();
167  fDParticleStartX = mctrack[0].X();
168  fDParticleStartE = mctrack[0].E();
169  fDParticleEndY = mctrack[nsteps-1].Y();
170  fDParticleEndZ = mctrack[nsteps-1].Z();
171  fDParticleEndX = mctrack[nsteps-1].X();
172  fDParticleEndE = mctrack[nsteps-1].E();
177 }
simb::Origin_t Origin() const
Definition: MCTrack.h:40
void SetOutputTree(TTree *, bool fill=true)
unsigned int AncestorTrackID() const
Definition: MCTrack.h:56
void FillTree(unsigned int, unsigned int, const std::vector< sim::MCTrack > &)
int AncestorPdgCode() const
Definition: MCTrack.h:55
unsigned int MotherTrackID() const
Definition: MCTrack.h:50
void fill(const art::PtrVector< recob::Hit > &hits, int only_plane)
Class def header for mctrack data container.
int PdgCode() const
Definition: MCTrack.h:41
int MotherPdgCode() const
Definition: MCTrack.h:49
unsigned int TrackID() const
Definition: MCTrack.h:42
void FillDominantParticleInfo(const sim::MCTrack &)
Event finding and building.