6 #include <HDDM/hddm_r.hpp>
20 bool foundLambda =
false;
23 bool foundProton =
false;
24 if(debug) cout <<
"event: " << nevents <<
"----------------------------------------------------------------" << endl;
26 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
27 hddm_r::ReactionList reactions = re.getReactions();
29 hddm_r::ReactionList::iterator iter;
30 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
31 int ireaction = std::distance(reactions.begin(), iter);
33 hddm_r::VertexList vertices = iter->getVertices();
34 hddm_r::VertexList::iterator iter_vertex;
36 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
37 int ivertex = std::distance(vertices.begin(), iter_vertex);
39 hddm_r::ProductList products = iter_vertex->getProducts();
40 hddm_r::ProductList::iterator iter_product;
42 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
43 int iparticle = std::distance(products.begin(), iter_product);
45 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
46 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
47 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
49 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
55 idLambda = iter_product->getId();
61 if(iter_product->getParentId() == idLambda){
65 idProton = iter_product->getId();
83 else if(select_type==2){
85 bool foundLambda =
false;
88 bool foundProton =
false;
90 bool foundEta =
false;
95 if(debug) cout <<
"event: " << nevents <<
"----------------------------------------------------------------" << endl;
97 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
98 hddm_r::ReactionList reactions = re.getReactions();
100 hddm_r::ReactionList::iterator iter;
101 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
102 int ireaction = std::distance(reactions.begin(), iter);
104 hddm_r::VertexList vertices = iter->getVertices();
105 hddm_r::VertexList::iterator iter_vertex;
107 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
108 int ivertex = std::distance(vertices.begin(), iter_vertex);
110 hddm_r::ProductList products = iter_vertex->getProducts();
111 hddm_r::ProductList::iterator iter_product;
113 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
114 int iparticle = std::distance(products.begin(), iter_product);
116 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
117 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
118 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
120 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
126 idLambda = iter_product->getId();
132 if(iter_product->getParentId() == idLambda){
136 idProton = iter_product->getId();
142 if(iter_product->getPdgtype() ==
PDGtype(
Eta)){
144 idEta = iter_product->getId();
150 if(iter_product->getParentId() == idEta){
153 idPhoton[nPhoton] = iter_product->getId();
161 if(foundProton && nPhoton==2){
172 else if(select_type==3){
173 const double VERTEXDIFFMAX = 10.0;
175 bool foundLambda =
false;
178 bool foundProton =
false;
179 double vertexLambda[3];
180 double vertexProton[3];
181 for(
int v=0;v<3;v++){
182 vertexLambda[v] = -999;
183 vertexProton[v] = 999;
185 double vertexdiff = 999;
187 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
188 hddm_r::ReactionList reactions = re.getReactions();
190 hddm_r::ReactionList::iterator iter;
191 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
192 int ireaction = std::distance(reactions.begin(), iter);
194 hddm_r::VertexList vertices = iter->getVertices();
195 hddm_r::VertexList::iterator iter_vertex;
197 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
198 int ivertex = std::distance(vertices.begin(), iter_vertex);
200 hddm_r::ProductList products = iter_vertex->getProducts();
201 hddm_r::ProductList::iterator iter_product;
203 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
204 int iparticle = std::distance(products.begin(), iter_product);
206 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
207 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
208 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
210 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
216 idLambda = iter_product->getId();
217 vertexLambda[0] = iter_vertex->getOrigin().getVx();
218 vertexLambda[1] = iter_vertex->getOrigin().getVy();
219 vertexLambda[2] = iter_vertex->getOrigin().getVz();
225 if(iter_product->getParentId() == idLambda){
229 idProton = iter_product->getId();
230 vertexProton[0] = iter_vertex->getOrigin().getVx();
231 vertexProton[1] = iter_vertex->getOrigin().getVy();
232 vertexProton[2] = iter_vertex->getOrigin().getVz();
233 vertexdiff =
sqrt(pow(vertexLambda[0] - vertexProton[0],2.) + pow(vertexLambda[1] - vertexProton[1],2.) + pow(vertexLambda[2] - vertexProton[2],2.));
241 if(foundProton && vertexdiff < VERTEXDIFFMAX){
252 else if(select_type==5){
253 bool foundProton =
false;
254 bool foundPip =
false;
255 bool foundPim =
false;
259 TLorentzVector p4photon_init;
261 TLorentzVector p4proton;
262 TLorentzVector p4pip;
263 TLorentzVector p4pim;
264 TLorentzVector p4diff;
266 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
267 hddm_r::ReactionList reactions = re.getReactions();
269 hddm_r::ReactionList::iterator iter;
270 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
271 int ireaction = std::distance(reactions.begin(), iter);
273 hddm_r::VertexList vertices = iter->getVertices();
274 hddm_r::VertexList::iterator iter_vertex;
276 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
277 int ivertex = std::distance(vertices.begin(), iter_vertex);
279 hddm_r::ProductList products = iter_vertex->getProducts();
280 hddm_r::ProductList::iterator iter_product;
282 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
283 int iparticle = std::distance(products.begin(), iter_product);
285 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
286 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
287 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
289 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
290 <<
" Ebeam: " << iter_vertex->getOrigin().getEbeam()
293 p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
298 idProton = iter_product->getId();
299 p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
305 idPip = iter_product->getId();
306 p4pip.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
312 idPim = iter_product->getId();
313 p4pim.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
321 if(foundProton && foundPip && foundPim){
325 p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim;
328 if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
337 else if(select_type==6){
342 Int_t nPimFromLambda = 0;
343 Int_t nPimFromXiMinus = 0;
345 int idXiMinus = -999;
347 TLorentzVector p4photon_init;
349 TLorentzVector p4kp1;
350 TLorentzVector p4kp2;
351 TLorentzVector p4XiMinus;
352 TLorentzVector p4Lambda;
353 TLorentzVector p4proton;
354 TLorentzVector p4pimFromXiMinus;
355 TLorentzVector p4pimFromLambda;
356 TLorentzVector p4diff;
357 TLorentzVector p4total;
359 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
360 hddm_r::ReactionList reactions = re.getReactions();
362 hddm_r::ReactionList::iterator iter;
363 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
364 int ireaction = std::distance(reactions.begin(), iter);
366 hddm_r::VertexList vertices = iter->getVertices();
367 hddm_r::VertexList::iterator iter_vertex;
369 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
370 int ivertex = std::distance(vertices.begin(), iter_vertex);
372 hddm_r::ProductList products = iter_vertex->getProducts();
373 hddm_r::ProductList::iterator iter_product;
375 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
376 int iparticle = std::distance(products.begin(), iter_product);
378 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
379 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
380 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
382 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
383 <<
" Ebeam: " << iter_vertex->getOrigin().getEbeam()
386 p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
389 if(nKp==0 && iter_product->getPdgtype() ==
PDGtype(
KPlus)){
391 p4kp1.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
393 cout <<
"K+ 1 : (" << p4kp1.X() <<
", " << p4kp1.Y() <<
", " << p4kp1.Z() <<
", " << p4kp1.M() <<
")" << endl;
401 if(nKp==1 && iter_product->getPdgtype() ==
PDGtype(
KPlus)){
403 p4kp1.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
405 cout <<
"K+ 2 : (" << p4kp2.X() <<
", " << p4kp2.Y() <<
", " << p4kp2.Z() <<
", " << p4kp2.M() <<
")" << endl;
411 idXiMinus = iter_product->getId();
412 p4XiMinus.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
416 if(iter_product->getPdgtype() ==
PDGtype(
PiMinus) && iter_product->getParentId() == idXiMinus){
418 p4pimFromXiMinus.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
422 if(iter_product->getPdgtype() ==
PDGtype(
Lambda) && iter_product->getParentId() == idXiMinus){
424 idLambda = iter_product->getId();
425 p4Lambda.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
429 if(iter_product->getPdgtype() ==
PDGtype(
Proton) && iter_product->getParentId() == idLambda){
431 p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
435 if(iter_product->getPdgtype() ==
PDGtype(
PiMinus) && iter_product->getParentId() == idLambda){
437 p4pimFromLambda.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
445 if(nKp==2 && nXiMinus==1 && nLambda==1 && nProton==1 && nPimFromXiMinus==1 && nPimFromLambda==1){
449 p4diff = p4photon_init + p4proton_init - p4kp1 - p4kp2 - p4proton - p4pimFromXiMinus - p4pimFromLambda;
450 p4total = p4kp1 + p4kp2 + p4proton + p4pimFromXiMinus + p4pimFromLambda;
453 cout <<
"photon: (" << p4photon_init.X() <<
", " << p4photon_init.Y() <<
", " << p4photon_init.Z() <<
", " << p4photon_init.M() <<
")" << endl;
454 cout <<
"proton: (" << p4proton_init.X() <<
", " << p4proton_init.Y() <<
", " << p4proton_init.Z() <<
", " << p4proton_init.M() <<
")" << endl;
455 cout <<
"Xi- : (" << p4XiMinus.X() <<
", " << p4XiMinus.Y() <<
", " << p4XiMinus.Z() <<
", " << p4XiMinus.M() <<
")" << endl;
456 cout <<
"Lambda: (" << p4Lambda.X() <<
", " << p4Lambda.Y() <<
", " << p4Lambda.Z() <<
", " << p4Lambda.M() <<
")" << endl;
457 cout <<
"K+ 1 : (" << p4kp1.X() <<
", " << p4kp1.Y() <<
", " << p4kp1.Z() <<
", " << p4kp1.M() <<
")" << endl;
458 cout <<
"K+ 2 : (" << p4kp2.X() <<
", " << p4kp2.Y() <<
", " << p4kp2.Z() <<
", " << p4kp2.M() <<
")" << endl;
459 cout <<
"proton: (" << p4proton.X() <<
", " << p4proton.Y() <<
", " << p4proton.Z() <<
", " << p4proton.M() <<
")" << endl;
460 cout <<
"pi- 1 : (" << p4pimFromXiMinus.X() <<
", " << p4pimFromXiMinus.Y() <<
", " << p4pimFromXiMinus.Z() <<
", " << p4pimFromXiMinus.M() <<
")" << endl;
461 cout <<
"pi- 2 : (" << p4pimFromLambda.X() <<
", " << p4pimFromLambda.Y() <<
", " << p4pimFromLambda.Z() <<
", " << p4pimFromLambda.M() <<
")" << endl;
465 if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
468 cout <<
"Added 4-mom of final state particles does not match generated by more than 0.1 MeV!!!!" << endl;
469 cout <<
"diff in |p| = " << p4diff.P()* 1000. <<
" MeV/c" << endl;
470 cout <<
"diff = (" << p4diff.X() <<
", " << p4diff.Y() <<
", " << p4diff.Z() <<
", " << p4diff.E() <<
")" << endl;
471 cout <<
"total = (" << p4total.X() <<
", " << p4total.Y() <<
", " << p4total.Z() <<
", " << p4total.E() <<
")" << endl;
480 else if(select_type==7){
481 bool foundProton =
false;
482 bool foundPip =
false;
483 bool foundPim =
false;
484 bool foundPi0 =
false;
489 TLorentzVector p4photon_init;
491 TLorentzVector p4proton;
492 TLorentzVector p4pip;
493 TLorentzVector p4pim;
494 TLorentzVector p4pi0;
495 TLorentzVector p4diff;
497 hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
498 hddm_r::ReactionList reactions = re.getReactions();
500 hddm_r::ReactionList::iterator iter;
501 for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
502 int ireaction = std::distance(reactions.begin(), iter);
504 hddm_r::VertexList vertices = iter->getVertices();
505 hddm_r::VertexList::iterator iter_vertex;
507 for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
508 int ivertex = std::distance(vertices.begin(), iter_vertex);
510 hddm_r::ProductList products = iter_vertex->getProducts();
511 hddm_r::ProductList::iterator iter_product;
513 for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
514 int iparticle = std::distance(products.begin(), iter_product);
516 if(debug) cout <<
"reaction: " << ireaction <<
" vertex: " << ivertex <<
" particle: " << iparticle
517 <<
" id: " << setw(3) << iter_product->getId() <<
" parent: " << iter_product->getParentId()
518 <<
" PDG type: " << setw(7) << iter_product->getPdgtype()
520 <<
" vertex: (" << iter_vertex->getOrigin().getVx() <<
", " << iter_vertex->getOrigin().getVy() <<
", " << iter_vertex->getOrigin().getVz() <<
")"
521 <<
" Ebeam: " << iter_vertex->getOrigin().getEbeam()
524 p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
529 idProton = iter_product->getId();
530 p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
536 idPip = iter_product->getId();
537 p4pip.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
543 idPim = iter_product->getId();
544 p4pim.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
548 if(iter_product->getPdgtype() ==
PDGtype(
Pi0)){
550 idPi0 = iter_product->getId();
551 p4pi0.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
560 if(foundProton && foundPip && foundPim && foundPi0){
564 p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim - p4pi0;
567 if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
static Particle_t PDGtoPType(int locPDG_PID)
bool selectEvent_r(int select_type, hddm_r::HDDM &record, int nevents, bool debug)
static char * ParticleType(Particle_t p)
static int PDGtype(Particle_t p)
static double ParticleMass(Particle_t p)