79 std::vector<float> *MC_Q2 = 0;
80 std::vector<float> *MC_W = 0;
81 std::vector<float> *MC_Y = 0;
82 std::vector<float> *MC_X = 0;
83 std::vector<float> *MC_Theta = 0;
84 std::vector<float> *MC_T = 0;
85 std::vector<float> *MCVertX = 0;
86 std::vector<float> *MCVertY = 0;
87 std::vector<float> *MCVertZ = 0;
88 std::vector<float> *MCNuPx = 0;
89 std::vector<float> *MCNuPy = 0;
90 std::vector<float> *MCNuPz = 0;
91 std::vector<int> *NType = 0;
92 std::vector<int> *CCNC = 0;
93 std::vector<int> *Mode = 0;
94 std::vector<int> *Gint=0;
95 std::vector<int> *TgtPDG=0;
96 std::vector<int> *
GT_T=0;
97 std::vector<int> *InterT=0;
98 std::vector<float> *Weight=0;
100 std::vector<int> *
PDG = 0;
101 std::vector<int> *Mother=0;
102 std::vector<int> *PDGMother=0;
103 std::vector<int> *MCPTrkID = 0;
104 std::vector<float> *MCPTime = 0;
105 std::vector<float> *MCPStartX = 0;
106 std::vector<float> *MCPStartY = 0;
107 std::vector<float> *MCPStartZ = 0;
108 std::vector<float> *MCPStartPX = 0;
109 std::vector<float> *MCPStartPY = 0;
110 std::vector<float> *MCPStartPZ = 0;
111 std::vector<std::string> *MCPProc = 0;
112 std::vector<std::string> *MCPEndProc = 0;
113 std::vector<float> *MCPEndX = 0;
114 std::vector<float> *MCPEndY = 0;
115 std::vector<float> *MCPEndZ = 0;
116 std::vector<float> *TrajMCPX = 0;
117 std::vector<float> *TrajMCPY = 0;
118 std::vector<float> *TrajMCPZ = 0;
119 std::vector<int> *TrajMCPTrajIndex = 0;
122 _inttree->SetBranchStatus(
"Event", 1);
123 _inttree->SetBranchStatus(
"SubRun", 1);
124 _inttree->SetBranchStatus(
"Run", 1);
126 _inttree->SetBranchStatus(
"NType", 1);
127 _inttree->SetBranchStatus(
"CCNC", 1);
128 _inttree->SetBranchStatus(
"MC_Q2", 1);
129 _inttree->SetBranchStatus(
"MC_W", 1);
130 _inttree->SetBranchStatus(
"MC_Y", 1);
131 _inttree->SetBranchStatus(
"MC_X", 1);
132 _inttree->SetBranchStatus(
"MC_Theta", 1);
133 _inttree->SetBranchStatus(
"MC_T", 1);
134 _inttree->SetBranchStatus(
"Mode", 1);
135 _inttree->SetBranchStatus(
"Gint", 1);
136 _inttree->SetBranchStatus(
"TgtPDG", 1);
137 _inttree->SetBranchStatus(
"GT_T", 1);
138 _inttree->SetBranchStatus(
"MCVertX", 1);
139 _inttree->SetBranchStatus(
"MCVertY", 1);
140 _inttree->SetBranchStatus(
"MCVertZ", 1);
141 _inttree->SetBranchStatus(
"MCNuPx", 1);
142 _inttree->SetBranchStatus(
"MCNuPy", 1);
143 _inttree->SetBranchStatus(
"MCNuPz", 1);
144 _inttree->SetBranchStatus(
"InterT", 1);
145 _inttree->SetBranchStatus(
"Weight", 1);
148 _inttree->SetBranchStatus(
"PDG", 1);
149 _inttree->SetBranchStatus(
"MCPTrkID", 1);
150 _inttree->SetBranchStatus(
"MCPTime", 1);
151 _inttree->SetBranchStatus(
"MCPStartX", 1);
152 _inttree->SetBranchStatus(
"MCPStartY", 1);
153 _inttree->SetBranchStatus(
"MCPStartZ", 1);
154 _inttree->SetBranchStatus(
"MCPEndX", 1);
155 _inttree->SetBranchStatus(
"MCPEndY", 1);
156 _inttree->SetBranchStatus(
"MCPEndZ", 1);
157 _inttree->SetBranchStatus(
"Mother", 1);
158 _inttree->SetBranchStatus(
"PDGMother", 1);
159 _inttree->SetBranchStatus(
"MCPStartPX", 1);
160 _inttree->SetBranchStatus(
"MCPStartPY", 1);
161 _inttree->SetBranchStatus(
"MCPStartPZ", 1);
162 _inttree->SetBranchStatus(
"MCPProc", 1);
163 _inttree->SetBranchStatus(
"MCPEndProc", 1);
164 _inttree->SetBranchStatus(
"TrajMCPX", 1);
165 _inttree->SetBranchStatus(
"TrajMCPY", 1);
166 _inttree->SetBranchStatus(
"TrajMCPZ", 1);
167 _inttree->SetBranchStatus(
"TrajMCPTrajIndex", 1);
171 _inttree->SetBranchAddress(
"Event", &Event);
172 _inttree->SetBranchAddress(
"SubRun", &SubRun);
173 _inttree->SetBranchAddress(
"Run", &Run);
176 _inttree->SetBranchAddress(
"NType", &NType);
177 _inttree->SetBranchAddress(
"CCNC", &CCNC);
178 _inttree->SetBranchAddress(
"MC_Q2", &MC_Q2);
179 _inttree->SetBranchAddress(
"MC_W", &MC_W);
180 _inttree->SetBranchAddress(
"MC_Y", &MC_Y);
181 _inttree->SetBranchAddress(
"MC_X", &MC_X);
182 _inttree->SetBranchAddress(
"MC_Theta", &MC_Theta);
183 _inttree->SetBranchAddress(
"MC_T", &MC_T);
184 _inttree->SetBranchAddress(
"Mode", &Mode);
185 _inttree->SetBranchAddress(
"Gint", &Gint);
186 _inttree->SetBranchAddress(
"TgtPDG", &TgtPDG);
187 _inttree->SetBranchAddress(
"GT_T", >_T);
188 _inttree->SetBranchAddress(
"MCVertX", &MCVertX);
189 _inttree->SetBranchAddress(
"MCVertY", &MCVertY);
190 _inttree->SetBranchAddress(
"MCVertZ", &MCVertZ);
191 _inttree->SetBranchAddress(
"MCNuPx", &MCNuPx);
192 _inttree->SetBranchAddress(
"MCNuPy", &MCNuPy);
193 _inttree->SetBranchAddress(
"MCNuPz", &MCNuPz);
194 _inttree->SetBranchAddress(
"InterT", &InterT);
195 _inttree->SetBranchAddress(
"Weight", &Weight);
198 _inttree->SetBranchAddress(
"PDG", &PDG);
199 _inttree->SetBranchAddress(
"MCPTrkID", &MCPTrkID);
200 _inttree->SetBranchAddress(
"MCPTime", &MCPTime);
201 _inttree->SetBranchAddress(
"MCPStartX", &MCPStartX);
202 _inttree->SetBranchAddress(
"MCPStartY", &MCPStartY);
203 _inttree->SetBranchAddress(
"MCPStartZ", &MCPStartZ);
204 _inttree->SetBranchAddress(
"MCPEndX", &MCPEndX);
205 _inttree->SetBranchAddress(
"MCPEndY", &MCPEndY);
206 _inttree->SetBranchAddress(
"MCPEndZ", &MCPEndZ);
207 _inttree->SetBranchAddress(
"Mother", &Mother);
208 _inttree->SetBranchAddress(
"PDGMother", &PDGMother);
209 _inttree->SetBranchAddress(
"MCPStartPX", &MCPStartPX);
210 _inttree->SetBranchAddress(
"MCPStartPY", &MCPStartPY);
211 _inttree->SetBranchAddress(
"MCPStartPZ", &MCPStartPZ);
212 _inttree->SetBranchAddress(
"MCPProc", &MCPProc);
213 _inttree->SetBranchAddress(
"MCPEndProc", &MCPEndProc);
214 _inttree->SetBranchAddress(
"TrajMCPX", &TrajMCPX);
215 _inttree->SetBranchAddress(
"TrajMCPY", &TrajMCPY);
216 _inttree->SetBranchAddress(
"TrajMCPZ", &TrajMCPZ);
217 _inttree->SetBranchAddress(
"TrajMCPTrajIndex", &TrajMCPTrajIndex);
227 std::cout <<
"Treating Evt " << Event <<
std::endl;
233 for(
size_t i = 0; i < NType->size(); i++)
235 _CCNC.push_back(CCNC->at(i));
236 _NType.push_back(NType->at(i));
237 _MC_Q2.push_back(MC_Q2->at(i));
238 _MC_W.push_back(MC_W->at(i));
239 _MC_Y.push_back(MC_Y->at(i));
240 _MC_X.push_back(MC_X->at(i));
242 _Mode.push_back(Mode->at(i));
243 _MC_T.push_back(MC_T->at(i));
244 _InterT.push_back(InterT->at(i));
248 _MCNuPx.push_back(MCNuPx->at(i));
249 _MCNuPy.push_back(MCNuPy->at(i));
250 _MCNuPz.push_back(MCNuPz->at(i));
253 for(
size_t i = 0; i < Gint->size(); i++)
255 _Gint.push_back(Gint->at(i));
256 _TgtPDG.push_back(TgtPDG->at(i));
257 _GT_T.push_back(GT_T->at(i));
258 _Weight.push_back(Weight->at(i));
261 std::map<int, int> mcp_kid_mother_map;
263 for(
size_t i = 0; i < MCPStartPX->size(); i++ )
267 int this_kid_trkid = MCPTrkID->at(i);
268 int this_kid_mother = Mother->at(i);
272 if(this_kid_mother == 0) {
278 for(
size_t j = 0; j < MCPStartPX->size(); j++ )
280 int mother_trk_id = MCPTrkID->at(j);
281 int mother = Mother->at(j);
284 if(mother_trk_id == this_kid_mother)
292 mcp_kid_mother_map[this_kid_trkid] = this_kid_mother;
297 this_kid_mother = mother;
305 std::vector<int> IndexesToKeep;
308 for(
size_t i = 0; i < MCPStartPX->size(); i++)
310 TVector3 spoint(MCPStartX->at(i), MCPStartY->at(i), MCPStartZ->at(i));
311 TVector3 epoint(MCPEndX->at(i), MCPEndY->at(i), MCPEndZ->at(i));
316 if( mcp_kid_mother_map.find(MCPTrkID->at(i)) != mcp_kid_mother_map.end() )
326 std::cout <<
"Keeping D0 or V0" <<
std::endl;
327 std::cout <<
"Index " << i <<
" TrkID " << MCPTrkID->at(i) <<
" pdg " << PDG->at(i) <<
" mother pdg " << PDGMother->at(i);
328 std::cout <<
" process " << process <<
" endprocess " << endprocess <<
std::endl;
329 std::cout <<
" Start point X: " << spoint.X() <<
" Y: " << spoint.Y() <<
" Z: " << spoint.Z() <<
std::endl;
330 std::cout <<
" End point X: " << epoint.X() <<
" Y: " << epoint.Y() <<
" Z: " << epoint.Z() <<
std::endl;
334 IndexesToKeep.push_back(i);
339 std::cout <<
"Keeping Backscatter" <<
std::endl;
340 std::cout <<
"Index " << i <<
" TrkID " << MCPTrkID->at(i) <<
" pdg " << PDG->at(i) <<
" mother pdg " << PDGMother->at(i);
341 std::cout <<
" process " << process <<
" endprocess " << endprocess <<
std::endl;
342 std::cout <<
" Start point X: " << spoint.X() <<
" Y: " << spoint.Y() <<
" Z: " << spoint.Z() <<
std::endl;
343 std::cout <<
" End point X: " << epoint.X() <<
" Y: " << epoint.Y() <<
" Z: " << epoint.Z() <<
std::endl;
347 IndexesToKeep.push_back(i);
352 std::cout <<
"Keeping Bremsstrahlung" <<
std::endl;
353 std::cout <<
"Index " << i <<
" TrkID " << MCPTrkID->at(i) <<
" pdg " << PDG->at(i) <<
" mother pdg " << PDGMother->at(i);
354 std::cout <<
" process " << process <<
" endprocess " << endprocess <<
std::endl;
355 std::cout <<
" Start point X: " << spoint.X() <<
" Y: " << spoint.Y() <<
" Z: " << spoint.Z() <<
std::endl;
356 std::cout <<
" End point X: " << epoint.X() <<
" Y: " << epoint.Y() <<
" Z: " << epoint.Z() <<
std::endl;
360 IndexesToKeep.push_back(i);
368 if(MCPProc->at(i) ==
"primary"){
370 std::cout <<
"Keeping Primary particle" <<
std::endl;
371 std::cout <<
"Index " << i <<
" TrkID " << MCPTrkID->at(i) <<
" pdg " << PDG->at(i) <<
" mother pdg " << PDGMother->at(i);
372 std::cout <<
" process " << process <<
" endprocess " << endprocess <<
std::endl;
373 std::cout <<
" Start point X: " << spoint.X() <<
" Y: " << spoint.Y() <<
" Z: " << spoint.Z() <<
std::endl;
374 std::cout <<
" End point X: " << epoint.X() <<
" Y: " << epoint.Y() <<
" Z: " << epoint.Z() <<
std::endl;
378 IndexesToKeep.push_back(i);
383 for(
size_t i = 0; i < IndexesToKeep.size(); i++)
385 _PDG.push_back(PDG->at(i));
399 _Mother.push_back(Mother->at(i));
402 for(
size_t itraj = 0; itraj < TrajMCPX->size(); itraj++){
404 if(TrajMCPTrajIndex->at(itraj) == (
int) i){
405 _TrajMCPX.push_back(TrajMCPX->at(itraj));
406 _TrajMCPY.push_back(TrajMCPY->at(itraj));
407 _TrajMCPZ.push_back(TrajMCPZ->at(itraj));
bool isBremsstrahlung(const TVector3 &spoint, const int &pdg, const int &motherpdg)
std::vector< std::string > _MCPProc
std::vector< float > _MCPStartX
std::vector< float > _MCPEndZ
std::vector< float > _MCPStartZ
std::vector< float > _Weight
std::vector< float > _MCPStartPZ
std::vector< float > _MCNuPy
std::vector< float > _TrajMCPY
std::vector< int > _MCPTrkID
std::vector< float > _MCPTime
std::vector< float > _MCPStartY
bool hasDecayedInCalo(const TVector3 &point)
std::vector< float > _MCVertY
std::vector< float > _TrajMCPX
std::vector< int > _TrajMCPTrajIndex
std::vector< float > _TrajMCPZ
std::vector< float > _MCPEndX
std::vector< float > _MC_Y
std::vector< float > _MCPEndY
bool PointInTPC(const TVector3 &point)
std::vector< int > daughtersToKeep
std::vector< float > _MC_T
std::vector< float > _MCPStartPX
std::vector< int > _InterT
std::vector< std::string > _MCPEndProc
bool isBackscatter(const TVector3 &spoint, const TVector3 &epoint)
std::vector< int > _NType
std::vector< float > _MC_W
std::vector< int > _TgtPDG
std::vector< float > _MCVertX
std::vector< int > _Mother
std::vector< float > _MC_X
std::vector< float > _MC_Theta
std::vector< float > _MCNuPz
std::vector< int > _PDGMother
std::vector< float > _MC_Q2
std::vector< float > _MCVertZ
std::vector< float > _MCNuPx
std::vector< float > _MCPStartPY
QTextStream & endl(QTextStream &s)