10 extern TRandom2 *
rndm;
12 #define SQR(X) ((X)*(X))
18 int nevents,
bool debug)
20 if (record.getPhysicsEvents().size() == 0) {
21 std::cout <<
"PhysicsEvent not given" << std::endl;
26 if (select_type == 1) {
27 bool foundLambda =
false;
30 bool foundProton =
false;
32 std::cout <<
"---------------------------------"
33 <<
"---------------------------------"
36 hddm_s::VertexList vs = record.getVertices();
37 hddm_s::VertexList::iterator viter;
38 for (viter = vs.begin(); viter != vs.end(); ++viter) {
39 hddm_s::OriginList orig = viter->getOrigins();
40 hddm_s::ProductList ps = viter->getProducts();
41 if (ps.size() > 0 && orig.size() > 0) {
42 hddm_s::ProductList::iterator piter;
43 for (piter = ps.begin(); piter != ps.end(); ++piter) {
44 int type = (int)piter->getType();
45 int id = piter->getId();
46 int parentid = piter->getParentid();
47 double E = piter->getMomentum().getE();
48 double px = piter->getMomentum().getPx();
49 double py = piter->getMomentum().getPy();
50 double pz = piter->getMomentum().getPz();
51 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
55 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
56 std::cout <<
"\t\t\ttype = " << type <<
" "
58 <<
" parentid = " << parentid << std::endl;
59 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
60 << E <<
"), mass = " << mass << std::endl;
61 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
62 << orig().getVy() <<
"," << orig().getVz() <<
","
63 << orig().getT() <<
")" << std::endl;
75 if (parentid == idLambda) {
96 else if (select_type == 2) {
97 bool foundLambda =
false;
100 bool foundProton =
false;
101 bool foundEta =
false;
106 std::cout <<
"------------------------------"
107 "------------------------------------"
110 hddm_s::VertexList vs = record.getVertices();
111 hddm_s::VertexList::iterator viter;
112 for (viter = vs.begin(); viter != vs.end(); ++viter) {
113 hddm_s::OriginList orig = viter->getOrigins();
114 hddm_s::ProductList ps = viter->getProducts();
115 if (ps.size() > 0 && orig.size() > 0) {
116 hddm_s::ProductList::iterator piter;
117 for (piter = ps.begin(); piter != ps.end(); ++piter) {
118 int type = (int)piter->getType();
119 int id = piter->getId();
120 int parentid = piter->getParentid();
121 double E = piter->getMomentum().getE();
122 double px = piter->getMomentum().getPx();
123 double py = piter->getMomentum().getPy();
124 double pz = piter->getMomentum().getPz();
125 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
129 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
130 std::cout <<
"\t\t\ttype = " << type <<
" "
132 <<
" parentid = " << parentid << std::endl;
133 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
134 << E <<
"), mass = " << mass << std::endl;
135 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
136 << orig().getVy() <<
"," << orig().getVz() <<
","
137 << orig().getT() <<
")" << std::endl;
149 if (parentid == idLambda) {
167 if (parentid == idEta) {
170 idPhoton[nPhoton] = id;
181 if (foundProton && nPhoton == 2) {
189 else if (select_type == 3) {
190 const double VERTEXDIFFMAX = 2.0;
191 bool foundLambda =
false;
194 bool foundProton =
false;
195 double vertexLambda[3];
196 double vertexProton[3];
197 for (
int v=0; v<3; v++) {
198 vertexLambda[v] = -999;
199 vertexProton[v] = 999;
201 double vertexdiff = 999;
203 std::cout <<
"----------------------------------"
204 "--------------------------------" << std::endl;
206 hddm_s::VertexList vs = record.getVertices();
207 hddm_s::VertexList::iterator viter;
208 for (viter = vs.begin(); viter != vs.end(); ++viter) {
209 hddm_s::OriginList orig = viter->getOrigins();
210 hddm_s::ProductList ps = viter->getProducts();
211 if (ps.size() > 0 && orig.size() > 0) {
212 hddm_s::ProductList::iterator piter;
213 for (piter = ps.begin(); piter != ps.end(); ++piter) {
214 int type = (int)piter->getType();
215 int id = piter->getId();
216 int parentid = piter->getParentid();
217 double E = piter->getMomentum().getE();
218 double px = piter->getMomentum().getPx();
219 double py = piter->getMomentum().getPy();
220 double pz = piter->getMomentum().getPz();
221 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
225 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
226 std::cout <<
"\t\t\ttype = " << type <<
" "
228 <<
" parentid = " << parentid << std::endl;
229 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
230 << E <<
"), mass = " << mass << std::endl;
231 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
232 << orig().getVy() <<
"," << orig().getVz() <<
","
233 << orig().getT() <<
")" << std::endl;
240 vertexLambda[0] = orig().getVx();
241 vertexLambda[1] = orig().getVy();
242 vertexLambda[2] = orig().getVz();
248 if (parentid == idLambda) {
253 vertexProton[0] = orig().getVx();
254 vertexProton[1] = orig().getVy();
255 vertexProton[2] = orig().getVz();
256 vertexdiff =
sqrt(
SQR(vertexLambda[0] - vertexProton[0]) +
257 SQR(vertexLambda[1] - vertexProton[1]) +
258 SQR(vertexLambda[2] - vertexProton[2]));
266 if (foundProton && vertexdiff < VERTEXDIFFMAX) {
277 else if (select_type == 4) {
278 bool foundKPlus =
false;
279 bool foundLambda =
false;
282 bool foundProton =
false;
283 TLorentzVector p4photon_init;
286 TLorentzVector p4Lambda;
287 TLorentzVector p4proton;
288 TVector3 beamdir(0,0,1);
292 std::cout <<
"---------------------------------"
293 "---------------------------------" << std::endl;
295 hddm_s::ReactionList res = record.getReactions();
296 hddm_s::ReactionList::iterator riter;
297 for (riter = res.begin(); riter != res.end(); ++riter) {
299 p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
300 riter->getBeam().getMomentum().getPy(),
301 riter->getBeam().getMomentum().getPz(),
302 riter->getBeam().getMomentum().getE());
304 std::cout <<
"photon. : (" << setw(5) << p4photon_init.X()
305 <<
", " << setw(5) << p4photon_init.Y() <<
", "
306 << p4photon_init.Z() <<
", "
307 << setw(5) << p4photon_init.E() <<
")" << std::endl;
309 hddm_s::VertexList vs = res().getVertices();
310 hddm_s::VertexList::iterator viter;
311 for (viter = vs.begin(); viter != vs.end(); ++viter) {
312 hddm_s::OriginList orig = viter->getOrigins();
313 hddm_s::ProductList ps = viter->getProducts();
314 if (ps.size() > 0 && orig.size() > 0) {
315 hddm_s::ProductList::iterator piter;
316 for (piter = ps.begin(); piter != ps.end(); ++piter) {
317 int type = (int)piter->getType();
318 int id = piter->getId();
319 int parentid = piter->getParentid();
320 double E = piter->getMomentum().getE();
321 double px = piter->getMomentum().getPx();
322 double py = piter->getMomentum().getPy();
323 double pz = piter->getMomentum().getPz();
324 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
328 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
329 std::cout <<
"\t\t\ttype = " << type <<
" "
331 <<
" parentid = " << parentid << std::endl;
332 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
333 << E <<
"), mass = " << mass << std::endl;
334 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
335 << orig().getVy() <<
"," << orig().getVz() <<
","
336 << orig().getT() <<
")" << std::endl;
342 p4kp.SetXYZT(px,py,pz,E);
343 normaldir = p4photon_init.Vect().Cross(p4kp.Vect()).Unit();
344 phi_n = normaldir.Phi();
347 std::cout <<
"K+ dir. : (" << setw(5) << p4kp.X()
348 <<
", " << setw(5) << p4kp.Y() <<
", "
349 << p4kp.Z() <<
")" << std::endl;
350 std::cout <<
"normal dir.: (" << setw(5)
351 << normaldir.X() <<
", " << setw(5)
352 << normaldir.Y() <<
", "
353 << normaldir.Z() <<
")" << std::endl;
361 p4Lambda.SetXYZT(px,py,pz,E);
367 if (parentid == idLambda) {
372 p4proton.SetXYZT(px,py,pz,E);
374 std::cout <<
"proton : ("
375 << setw(5) << p4proton.X() <<
", "
376 << setw(5) << p4proton.Y() <<
", "
377 << p4proton.Z() <<
")" << std::endl;
388 if (foundKPlus && foundProton) {
391 p4Lambda.Boost(-(p4photon_init + p4proton_init).BoostVector());
392 p4proton.Boost(-(p4photon_init + p4proton_init).BoostVector());
396 p4Lambda.RotateZ(TMath::Pi()/2. - phi_n);
397 p4proton.RotateZ(TMath::Pi()/2. - phi_n);
401 p4Lambda.RotateX(TMath::Pi()/2);
402 p4proton.RotateX(TMath::Pi()/2);
405 p4proton.Boost(-p4Lambda.BoostVector());
407 double costheta_proton = p4proton.CosTheta();
410 double rand =
rndm->Rndm();
412 std::cout <<
"random number: " << rand << std::endl;
413 std::cout <<
"compare to : "
414 << (1. +
alpha * costheta_proton) / (1. +
alpha)
417 if (rand < (1. +
alpha * costheta_proton) / (1. +
alpha)) {
418 std::cout <<
" --- " << costheta_proton << std::endl;
426 else if (select_type == 5) {
427 bool foundProton =
false;
428 bool foundPip =
false;
429 bool foundPim =
false;
433 TLorentzVector p4photon_init;
435 TLorentzVector p4proton;
436 TLorentzVector p4pip;
437 TLorentzVector p4pim;
438 TLorentzVector p4diff;
440 std::cout <<
"----------------------------------"
441 "--------------------------------" << std::endl;
443 hddm_s::ReactionList res = record.getReactions();
444 hddm_s::ReactionList::iterator riter;
445 for (riter = res.begin(); riter != res.end(); ++riter) {
447 p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
448 riter->getBeam().getMomentum().getPy(),
449 riter->getBeam().getMomentum().getPz(),
450 riter->getBeam().getMomentum().getE());
452 std::cout <<
"photon. : (" << setw(5)
453 << p4photon_init.X() <<
", " << setw(5)
454 << p4photon_init.Y() <<
", " << p4photon_init.Z()
455 <<
", " << setw(5) << p4photon_init.E() <<
")"
458 hddm_s::VertexList vs = res().getVertices();
459 hddm_s::VertexList::iterator viter;
460 for (viter = vs.begin(); viter != vs.end(); ++viter) {
461 hddm_s::OriginList orig = viter->getOrigins();
462 hddm_s::ProductList ps = viter->getProducts();
463 if (ps.size() > 0 && orig.size() > 0) {
464 hddm_s::ProductList::iterator piter;
465 for (piter = ps.begin(); piter != ps.end(); ++piter) {
466 int type = (int)piter->getType();
467 int id = piter->getId();
468 int parentid = piter->getParentid();
469 double E = piter->getMomentum().getE();
470 double px = piter->getMomentum().getPx();
471 double py = piter->getMomentum().getPy();
472 double pz = piter->getMomentum().getPz();
473 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
477 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
478 std::cout <<
"\t\t\ttype = " << type <<
" "
480 <<
" parentid = " << parentid << std::endl;
481 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
482 << E <<
"), mass = " << mass << std::endl;
483 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
484 << orig().getVy() <<
"," << orig().getVz() <<
","
485 << orig().getT() <<
")" << std::endl;
492 p4proton.SetXYZT(px,py,pz,E);
499 p4pip.SetXYZT(px,py,pz,E);
506 p4pim.SetXYZT(px,py,pz,E);
515 if (foundProton && foundPip && foundPim) {
519 p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim;
522 if (p4diff.P() < 0.001 && p4diff.M() < 0.0001)
530 else if (select_type == 6) {
535 Int_t nPimFromLambda = 0;
536 Int_t nPimFromXiMinus = 0;
538 int idXiMinus = -999;
540 TLorentzVector p4photon_init;
542 TLorentzVector p4kp1;
543 TLorentzVector p4kp2;
544 TLorentzVector p4XiMinus;
545 TLorentzVector p4Lambda;
546 TLorentzVector p4proton;
547 TLorentzVector p4pimFromXiMinus;
548 TLorentzVector p4pimFromLambda;
549 TLorentzVector p4diff;
550 TLorentzVector p4total;
552 std::cout <<
"-------------------------------"
553 "-----------------------------------" << std::endl;
555 hddm_s::ReactionList res = record.getReactions();
556 hddm_s::ReactionList::iterator riter;
557 for (riter = res.begin(); riter != res.end(); ++riter) {
559 p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
560 riter->getBeam().getMomentum().getPy(),
561 riter->getBeam().getMomentum().getPz(),
562 riter->getBeam().getMomentum().getE());
564 std::cout <<
"photon. : (" << setw(5)
565 << p4photon_init.X() <<
", " << setw(5)
566 << p4photon_init.Y() <<
", " << p4photon_init.Z()
567 <<
", " << setw(5) << p4photon_init.E() <<
")"
570 hddm_s::VertexList vs = res().getVertices();
571 hddm_s::VertexList::iterator viter;
572 for (viter = vs.begin(); viter != vs.end(); ++viter) {
573 hddm_s::OriginList orig = viter->getOrigins();
574 hddm_s::ProductList ps = viter->getProducts();
575 if (ps.size() > 0 && orig.size() > 0) {
576 hddm_s::ProductList::iterator piter;
577 for (piter = ps.begin(); piter != ps.end(); ++piter) {
578 int type = (int)piter->getType();
579 int id = piter->getId();
580 int parentid = piter->getParentid();
581 double E = piter->getMomentum().getE();
582 double px = piter->getMomentum().getPx();
583 double py = piter->getMomentum().getPy();
584 double pz = piter->getMomentum().getPz();
585 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
589 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
590 std::cout <<
"\t\t\ttype = " << type <<
" "
592 <<
" parentid = " << parentid << std::endl;
593 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
594 << E <<
"), mass = " << mass << std::endl;
595 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
596 << orig().getVy() <<
"," << orig().getVz() <<
","
597 << orig().getT() <<
")" << std::endl;
601 if (nKp == 0 && type ==
KPlus) {
603 p4kp1.SetXYZT(px,py,pz,E);
605 std::cout <<
"K+ 1 : (" << p4kp1.X() <<
", "
606 << p4kp1.Y() <<
", " << p4kp1.Z() <<
", "
607 << p4kp1.M() <<
")" << std::endl;
615 if (nKp == 1 && type ==
KPlus) {
617 p4kp2.SetXYZT(px,py,pz,E);
619 std::cout <<
"K+ 2 : (" << p4kp2.X() <<
", "
620 << p4kp2.Y() <<
", " << p4kp2.Z() <<
", "
621 << p4kp2.M() <<
")" << std::endl;
628 p4XiMinus.SetXYZT(px,py,pz,E);
632 if (type ==
PiMinus && parentid == idXiMinus) {
634 p4pimFromXiMinus.SetXYZT(px,py,pz,E);
638 if (type ==
Lambda && parentid == idXiMinus) {
641 p4Lambda.SetXYZT(px,py,pz,E);
645 if (type ==
Proton && parentid == idLambda) {
647 p4proton.SetXYZT(px,py,pz,E);
651 if (type ==
PiMinus && parentid == idLambda) {
653 p4pimFromLambda.SetXYZT(px,py,pz,E);
662 if (nKp == 2 && nXiMinus == 1 && nLambda == 1 && nProton == 1 &&
663 nPimFromXiMinus == 1 && nPimFromLambda == 1)
665 std::cout <<
"found all particles" << std::endl;
669 p4diff = p4photon_init + p4proton_init - p4kp1 - p4kp2 -
670 p4proton - p4pimFromXiMinus - p4pimFromLambda;
671 p4total = p4kp1 + p4kp2 + p4proton + p4pimFromXiMinus + p4pimFromLambda;
674 std::cout <<
"photon: (" << p4photon_init.X() <<
", "
675 << p4photon_init.Y() <<
", " << p4photon_init.Z()
676 <<
", " << p4photon_init.M() <<
")" << std::endl;
677 std::cout <<
"proton: (" << p4proton_init.X() <<
", "
678 << p4proton_init.Y() <<
", " << p4proton_init.Z()
679 <<
", " << p4proton_init.M() <<
")" << std::endl;
680 std::cout <<
"Xi- : (" << p4XiMinus.X() <<
", "
681 << p4XiMinus.Y() <<
", " << p4XiMinus.Z()
682 <<
", " << p4XiMinus.M() <<
")" << std::endl;
683 std::cout <<
"Lambda: (" << p4Lambda.X() <<
", "
684 << p4Lambda.Y() <<
", " << p4Lambda.Z() <<
", "
685 << p4Lambda.M() <<
")" << std::endl;
686 std::cout <<
"K+ 1 : (" << p4kp1.X() <<
", " << p4kp1.Y()
687 <<
", " << p4kp1.Z() <<
", " << p4kp1.M() <<
")"
689 std::cout <<
"K+ 2 : (" << p4kp2.X() <<
", " << p4kp2.Y()
690 <<
", " << p4kp2.Z() <<
", " << p4kp2.M() <<
")"
692 std::cout <<
"proton: (" << p4proton.X() <<
", " << p4proton.Y()
693 <<
", " << p4proton.Z() <<
", " << p4proton.M() <<
")"
695 std::cout <<
"pi- 1 : (" << p4pimFromXiMinus.X() <<
", "
696 << p4pimFromXiMinus.Y() <<
", " << p4pimFromXiMinus.Z()
697 <<
", " << p4pimFromXiMinus.M() <<
")" << std::endl;
698 std::cout <<
"pi- 2 : (" << p4pimFromLambda.X() <<
", "
699 << p4pimFromLambda.Y() <<
", " << p4pimFromLambda.Z()
700 <<
", " << p4pimFromLambda.M() <<
")" << std::endl;
704 if (p4diff.P() < 0.001 && p4diff.M() < 0.0001) {
705 std::cout <<
"returning true" << std::endl;
709 std::cout <<
"Added 4-mom of final state particles does not match"
710 " generated by more than 0.1 MeV!!!!" << std::endl;
711 std::cout <<
"diff in |p| = " << p4diff.P()* 1000. <<
" MeV/c"
713 std::cout <<
"diff = (" << p4diff.X() <<
", " << p4diff.Y()
714 <<
", " << p4diff.Z() <<
", " << p4diff.E() <<
")"
716 std::cout <<
"total = (" << p4total.X() <<
", " << p4total.Y()
717 <<
", " << p4total.Z() <<
", " << p4total.E() <<
")"
726 else if (select_type == 7) {
727 bool foundProton =
false;
728 bool foundPip =
false;
729 bool foundPim =
false;
730 bool foundPi0 =
false;
735 TLorentzVector p4photon_init;
737 TLorentzVector p4proton;
738 TLorentzVector p4pip;
739 TLorentzVector p4pim;
740 TLorentzVector p4pi0;
741 TLorentzVector p4diff;
743 std::cout <<
"---------------------------------"
744 "---------------------------------" << std::endl;
746 hddm_s::ReactionList res = record.getReactions();
747 hddm_s::ReactionList::iterator riter;
748 for (riter = res.begin(); riter != res.end(); ++riter) {
750 p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
751 riter->getBeam().getMomentum().getPy(),
752 riter->getBeam().getMomentum().getPz(),
753 riter->getBeam().getMomentum().getE());
755 std::cout <<
"photon. : (" << setw(5)
756 << p4photon_init.X() <<
", " << setw(5)
757 << p4photon_init.Y() <<
", " << p4photon_init.Z()
758 <<
", " << setw(5) << p4photon_init.E() <<
")"
761 hddm_s::VertexList vs = res().getVertices();
762 hddm_s::VertexList::iterator viter;
763 for (viter = vs.begin(); viter != vs.end(); ++viter) {
764 hddm_s::OriginList orig = viter->getOrigins();
765 hddm_s::ProductList ps = viter->getProducts();
766 if (ps.size() > 0 && orig.size() > 0) {
767 hddm_s::ProductList::iterator piter;
768 for (piter = ps.begin(); piter != ps.end(); ++piter) {
769 int type = (int)piter->getType();
770 int id = piter->getId();
771 int parentid = piter->getParentid();
772 double E = piter->getMomentum().getE();
773 double px = piter->getMomentum().getPx();
774 double py = piter->getMomentum().getPy();
775 double pz = piter->getMomentum().getPz();
776 double mass =
sqrt(E*E - (px*px + py*py + pz*pz));
780 std::cout <<
"\t\t\t--- new particle ---" << std::endl;
781 std::cout <<
"\t\t\ttype = " << type <<
" "
783 <<
" parentid = " << parentid << std::endl;
784 std::cout <<
"\t\t\t(" << px <<
"," << py <<
"," << pz <<
","
785 << E <<
"), mass = " << mass << std::endl;
786 std::cout <<
"\t\t\t(" << orig().getVx() <<
","
787 << orig().getVy() <<
"," << orig().getVz() <<
","
788 << orig().getT() <<
")" << std::endl;
795 p4proton.SetXYZT(px,py,pz,E);
802 p4pip.SetXYZT(px,py,pz,E);
809 p4pim.SetXYZT(px,py,pz,E);
816 p4pi0.SetXYZT(px,py,pz,E);
826 if (foundProton && foundPip && foundPim && foundPi0) {
830 p4diff = p4photon_init + p4proton_init - p4proton -
831 p4pip - p4pim - p4pi0;
834 if (p4diff.P() < 0.001 && p4diff.M() < 0.0001)
static char * ParticleType(Particle_t p)
static double ParticleMass(Particle_t p)
bool selectEvent_s(int select_type, hddm_s::HDDM &record, int nevents, bool debug)