Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
selectEvents_s.cc
Go to the documentation of this file.
1 // $Id$
2 //
3 // Created Oct 10, 2013 David Lawrence
4 
5 #include "hddm_select_events.h"
6 #include "particleType.h"
7 
8 using namespace std;
9 
10 extern TRandom2 *rndm;
11 
12 #define SQR(X) ((X)*(X))
13 
14 //-----------------------------------------------------------------
15 // selectEvent_s
16 //-----------------------------------------------------------------
17 bool selectEvent_s(int select_type, hddm_s::HDDM &record,
18  int nevents, bool debug)
19 {
20  if (record.getPhysicsEvents().size() == 0) {
21  std::cout << "PhysicsEvent not given" << std::endl;
22  return -1;
23  }
24 
25  // Select Lambda -> p pi- events
26  if (select_type == 1) {
27  bool foundLambda = false;
28  int idLambda = -999;
29  int idProton = -999;
30  bool foundProton = false;
31  if (debug)
32  std::cout << "---------------------------------"
33  << "---------------------------------"
34  << std::endl;
35 
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));
52  if (! finite(mass))
53  mass = 0.0;
54  if (debug) {
55  std::cout << "\t\t\t--- new particle ---" << std::endl;
56  std::cout << "\t\t\ttype = " << type << " "
57  << ParticleType((Particle_t)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;
64  }
65 
66  // Find Lambda
67  if (type == Lambda) {
68  foundLambda = true;
69  idLambda = id;
70  }
71 
72  // Find decay products of Lambda
73  if (foundLambda) {
74  // if parent was Lambda
75  if (parentid == idLambda) {
76  // find what type it is
77  if (type == Proton) {
78  foundProton = true;
79  idProton = id;
80  }
81  }
82  }
83 
84  // If we found a proton decaying from Lambda,
85  // save this event
86  if (foundProton) {
87  return true;
88  }
89  } // end of loop over products
90  } // end of having products and origin
91  } // end of loop over vertices
92  } // end of select_type 1
93  //___________________________________________________________________________
94 
95  // select Lambda -> p pi-, eta -> gamma gamma events
96  else if (select_type == 2) {
97  bool foundLambda = false;
98  int idLambda = -999;
99  int idProton = -999;
100  bool foundProton = false;
101  bool foundEta = false;
102  int idEta = -999;
103  int idPhoton[2];
104  int nPhoton = 0;
105  if (debug)
106  std::cout << "------------------------------"
107  "------------------------------------"
108  << std::endl;
109 
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));
126  if (! finite(mass))
127  mass = 0.0;
128  if (debug) {
129  std::cout << "\t\t\t--- new particle ---" << std::endl;
130  std::cout << "\t\t\ttype = " << type << " "
131  << ParticleType((Particle_t)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;
138  }
139 
140  // Find Lambda
141  if (type == Lambda) {
142  foundLambda = true;
143  idLambda = id;
144  }
145 
146  // Find decay products of Lambda
147  if (foundLambda) {
148  // if parent was Lambda
149  if (parentid == idLambda) {
150  // find what type it is
151  if (type == Proton) {
152  foundProton = true;
153  idProton = id;
154  }
155  }
156  }
157 
158  // Find eta
159  if (type == Eta) {
160  foundEta = true;
161  idEta = id;
162  }
163 
164  // Find decay products of Eta
165  if (foundEta) {
166  // if parent was Eta
167  if (parentid == idEta) {
168  // find what type it is
169  if (type == Gamma) {
170  idPhoton[nPhoton] = id;
171  nPhoton++;
172  }
173  }
174  }
175  }
176  }
177 
178  // If we found a proton decaying from Lambda,
179  // and we found 2 photons from eta decay,
180  // save this event
181  if (foundProton && nPhoton == 2) {
182  return true;
183  }
184  } // end of loop over vertices
185  }
186  //___________________________________________________________________________
187 
188  // Select Lambda -> p pi- events, decay length is less than VERTEXDIFFMAX
189  else if (select_type == 3) {
190  const double VERTEXDIFFMAX = 2.0;
191  bool foundLambda = false;
192  int idLambda = -999;
193  int idProton = -999;
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;
200  }
201  double vertexdiff = 999;
202  if (debug)
203  std::cout << "----------------------------------"
204  "--------------------------------" << std::endl;
205 
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));
222  if (! finite(mass))
223  mass = 0.0;
224  if (debug) {
225  std::cout << "\t\t\t--- new particle ---" << std::endl;
226  std::cout << "\t\t\ttype = " << type << " "
227  << ParticleType((Particle_t)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;
234  }
235 
236  // Find Lambda
237  if (type == Lambda) {
238  foundLambda = true;
239  idLambda = id;
240  vertexLambda[0] = orig().getVx();
241  vertexLambda[1] = orig().getVy();
242  vertexLambda[2] = orig().getVz();
243  }
244 
245  // Find decay products of Lambda
246  if (foundLambda) {
247  // if parent was Lambda
248  if (parentid == idLambda) {
249  // find what type it is
250  if (type == Proton) {
251  foundProton = true;
252  idProton = id;
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]));
259  }
260  }
261  }
262 
263  // If we found a proton decaying from Lambda,
264  // and the vertex difference is less than VERTEXDIFFMAX,
265  // save this event
266  if (foundProton && vertexdiff < VERTEXDIFFMAX) {
267  return true;
268  }
269 
270  } // end of loop over products
271  } // end of having products and origin
272  } // end of loop over vertices
273  } // end of select_type 3
274  //___________________________________________________________________________
275 
276  // Select Lambda -> p pi- events, carve out weak decay asymmetry
277  else if (select_type == 4) {
278  bool foundKPlus = false;
279  bool foundLambda = false;
280  int idLambda = -999;
281  int idProton = -999;
282  bool foundProton = false;
283  TLorentzVector p4photon_init;
284  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
285  TLorentzVector p4kp;
286  TLorentzVector p4Lambda;
287  TLorentzVector p4proton;
288  TVector3 beamdir(0,0,1);
289  TVector3 normaldir;
290  double phi_n = 0;
291  if (debug)
292  std::cout << "---------------------------------"
293  "---------------------------------" << std::endl;
294 
295  hddm_s::ReactionList res = record.getReactions();
296  hddm_s::ReactionList::iterator riter;
297  for (riter = res.begin(); riter != res.end(); ++riter) {
298  // Get initial photon
299  p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
300  riter->getBeam().getMomentum().getPy(),
301  riter->getBeam().getMomentum().getPz(),
302  riter->getBeam().getMomentum().getE());
303  if (debug)
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;
308 
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));
325  if (! finite(mass))
326  mass = 0.0;
327  if (debug) {
328  std::cout << "\t\t\t--- new particle ---" << std::endl;
329  std::cout << "\t\t\ttype = " << type << " "
330  << ParticleType((Particle_t)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;
337  }
338 
339  // Find K+, set up polarization axis
340  if (type == KPlus) {
341  foundKPlus = true;
342  p4kp.SetXYZT(px,py,pz,E);
343  normaldir = p4photon_init.Vect().Cross(p4kp.Vect()).Unit();
344  phi_n = normaldir.Phi();
345 
346  if (debug) {
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;
354  }
355  }
356 
357  // Find Lambda
358  if (type == Lambda) {
359  foundLambda = true;
360  idLambda = id;
361  p4Lambda.SetXYZT(px,py,pz,E);
362  }
363 
364  // Find decay products of Lambda
365  if (foundLambda) {
366  // if parent was Lambda
367  if (parentid == idLambda) {
368  // find what type it is
369  if (type == Proton) {
370  foundProton = true;
371  idProton = id;
372  p4proton.SetXYZT(px,py,pz,E);
373  if (debug)
374  std::cout << "proton : ("
375  << setw(5) << p4proton.X() << ", "
376  << setw(5) << p4proton.Y() << ", "
377  << p4proton.Z() << ")" << std::endl;
378  }
379  }
380  }
381  } // end of loop over products
382  } // end of having products and origin
383  } // end of loop over vertices
384  } // end of loop over reaction
385 
386  // If we found a proton decaying from Lambda,
387  // save this event
388  if (foundKPlus && foundProton) {
389  // Get proton angular dist
390  // First boost to c.m. frame
391  p4Lambda.Boost(-(p4photon_init + p4proton_init).BoostVector());
392  p4proton.Boost(-(p4photon_init + p4proton_init).BoostVector());
393 
394  // 2. Rotate around z by pi/2 - phi_n
395  // This takes K+ into xz plane, and normal into y-dir.
396  p4Lambda.RotateZ(TMath::Pi()/2. - phi_n);
397  p4proton.RotateZ(TMath::Pi()/2. - phi_n);
398 
399  // 3. Rotate around x by pi/2
400  // This takes K+ into xy plane, and normal into z-dir.
401  p4Lambda.RotateX(TMath::Pi()/2);
402  p4proton.RotateX(TMath::Pi()/2);
403 
404  // 4. Boost to Lambda rest frame.
405  p4proton.Boost(-p4Lambda.BoostVector());
406  // Now get decay angle
407  double costheta_proton = p4proton.CosTheta();
408 
409  // Throw random number and select events
410  double rand =rndm->Rndm();
411  if (debug) {
412  std::cout << "random number: " << rand << std::endl;
413  std::cout << "compare to : "
414  << (1. + alpha * costheta_proton) / (1. + alpha)
415  << std::endl;
416  }
417  if (rand < (1. + alpha * costheta_proton) / (1. + alpha)) {
418  std::cout << " --- " << costheta_proton << std::endl;
419  return true;
420  }
421  } // found K+ and proton
422  } // end of select_type 4
423  //___________________________________________________________________________
424 
425  // Select p pi+ pi- events
426  else if (select_type == 5) {
427  bool foundProton = false;
428  bool foundPip = false;
429  bool foundPim = false;
430  int idProton = -999;
431  int idPip = -999;
432  int idPim = -999;
433  TLorentzVector p4photon_init;
434  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
435  TLorentzVector p4proton;
436  TLorentzVector p4pip;
437  TLorentzVector p4pim;
438  TLorentzVector p4diff;
439  if (debug)
440  std::cout << "----------------------------------"
441  "--------------------------------" << std::endl;
442 
443  hddm_s::ReactionList res = record.getReactions();
444  hddm_s::ReactionList::iterator riter;
445  for (riter = res.begin(); riter != res.end(); ++riter) {
446  // Get initial photon
447  p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
448  riter->getBeam().getMomentum().getPy(),
449  riter->getBeam().getMomentum().getPz(),
450  riter->getBeam().getMomentum().getE());
451  if (debug)
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() << ")"
456  << std::endl;
457 
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));
474  if (! finite(mass))
475  mass = 0.0;
476  if (debug) {
477  std::cout << "\t\t\t--- new particle ---" << std::endl;
478  std::cout << "\t\t\ttype = " << type << " "
479  << ParticleType((Particle_t)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;
486  }
487 
488  // Find proton
489  if (type == Proton) {
490  foundProton = true;
491  idProton = id;
492  p4proton.SetXYZT(px,py,pz,E);
493  }
494 
495  // Find pi+
496  if (type == PiPlus) {
497  foundPip = true;
498  idPip = id;
499  p4pip.SetXYZT(px,py,pz,E);
500  }
501 
502  // Find pi-
503  if (type == PiMinus) {
504  foundPim = true;
505  idPim = id;
506  p4pim.SetXYZT(px,py,pz,E);
507  }
508  } // end of loop over products
509  } // end of having products and origin
510  } // end of loop over vertices
511  } // end of loop over reaction
512 
513  // If we found a proton decaying from Lambda,
514  // save this event
515  if (foundProton && foundPip && foundPim) {
516 
517  // make sure that total 4-mom of p, pi+, pi- are close to
518  // initial state
519  p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim;
520 
521  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
522  if (p4diff.P() < 0.001 && p4diff.M() < 0.0001)
523  return true;
524  } // found proton, pi+, pi-
525 
526  } // end of select_type 5
527  //___________________________________________________________________________
528 
529  // Select K+ K+ Xi- -> K+ K+ p pi- pi- events
530  else if (select_type == 6) {
531  Int_t nKp = 0;
532  Int_t nXiMinus = 0;
533  Int_t nLambda = 0;
534  Int_t nProton = 0;
535  Int_t nPimFromLambda = 0;
536  Int_t nPimFromXiMinus = 0;
537  // need intermediate IDs
538  int idXiMinus = -999;
539  int idLambda = -999;
540  TLorentzVector p4photon_init;
541  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
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;
551  if (debug)
552  std::cout << "-------------------------------"
553  "-----------------------------------" << std::endl;
554 
555  hddm_s::ReactionList res = record.getReactions();
556  hddm_s::ReactionList::iterator riter;
557  for (riter = res.begin(); riter != res.end(); ++riter) {
558  // Get initial photon
559  p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
560  riter->getBeam().getMomentum().getPy(),
561  riter->getBeam().getMomentum().getPz(),
562  riter->getBeam().getMomentum().getE());
563  if (debug)
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() << ")"
568  << std::endl;
569 
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));
586  if (! finite(mass))
587  mass = 0.0;
588  if (debug) {
589  std::cout << "\t\t\t--- new particle ---" << std::endl;
590  std::cout << "\t\t\ttype = " << type << " "
591  << ParticleType((Particle_t)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;
598  }
599 
600  // Find first K+
601  if (nKp == 0 && type == KPlus) {
602  nKp++;
603  p4kp1.SetXYZT(px,py,pz,E);
604  if (debug)
605  std::cout << "K+ 1 : (" << p4kp1.X() << ", "
606  << p4kp1.Y() << ", " << p4kp1.Z() << ", "
607  << p4kp1.M() << ")" << std::endl;
608 
609  // Need this break since without it we would count the
610  // the first K+ as the second K+ too
611  continue;
612  }
613 
614  // Find second K+
615  if (nKp == 1 && type == KPlus) {
616  nKp++;
617  p4kp2.SetXYZT(px,py,pz,E);
618  if (debug)
619  std::cout << "K+ 2 : (" << p4kp2.X() << ", "
620  << p4kp2.Y() << ", " << p4kp2.Z() << ", "
621  << p4kp2.M() << ")" << std::endl;
622  }
623 
624  // Find Xi-
625  if (type == XiMinus) {
626  nXiMinus++;
627  idXiMinus = id;
628  p4XiMinus.SetXYZT(px,py,pz,E);
629  }
630 
631  // Find pi- from Xi-
632  if (type == PiMinus && parentid == idXiMinus) {
633  nPimFromXiMinus++;
634  p4pimFromXiMinus.SetXYZT(px,py,pz,E);
635  }
636 
637  // Find Lambda from Xi-
638  if (type == Lambda && parentid == idXiMinus) {
639  nLambda++;
640  idLambda = id;
641  p4Lambda.SetXYZT(px,py,pz,E);
642  }
643 
644  // Find proton from Lambda
645  if (type == Proton && parentid == idLambda) {
646  nProton++;
647  p4proton.SetXYZT(px,py,pz,E);
648  }
649 
650  // Find pi- from Lambda
651  if (type == PiMinus && parentid == idLambda) {
652  nPimFromLambda++;
653  p4pimFromLambda.SetXYZT(px,py,pz,E);
654  }
655  } // end of loop over products
656  } // end of having products and origin
657  } // end of loop over vertices
658  } // end of loop over reaction
659 
660  // If we found a proton decaying from Lambda,
661  // save this event
662  if (nKp == 2 && nXiMinus == 1 && nLambda == 1 && nProton == 1 &&
663  nPimFromXiMinus == 1 && nPimFromLambda == 1)
664  {
665  std::cout << "found all particles" << std::endl;
666 
667  // make sure that total 4-mom of p, pi+, pi- are close to
668  // initial state
669  p4diff = p4photon_init + p4proton_init - p4kp1 - p4kp2 -
670  p4proton - p4pimFromXiMinus - p4pimFromLambda;
671  p4total = p4kp1 + p4kp2 + p4proton + p4pimFromXiMinus + p4pimFromLambda;
672 
673  if (debug) {
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() << ")"
688  << std::endl;
689  std::cout << "K+ 2 : (" << p4kp2.X() << ", " << p4kp2.Y()
690  << ", " << p4kp2.Z() << ", " << p4kp2.M() << ")"
691  << std::endl;
692  std::cout << "proton: (" << p4proton.X() << ", " << p4proton.Y()
693  << ", " << p4proton.Z() << ", " << p4proton.M() << ")"
694  << std::endl;
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;
701  }
702 
703  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
704  if (p4diff.P() < 0.001 && p4diff.M() < 0.0001) {
705  std::cout << "returning true" << std::endl;
706  return true;
707  }
708  else {
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"
712  << std::endl;
713  std::cout << "diff = (" << p4diff.X() << ", " << p4diff.Y()
714  << ", " << p4diff.Z() << ", " << p4diff.E() << ")"
715  << std::endl;
716  std::cout << "total = (" << p4total.X() << ", " << p4total.Y()
717  << ", " << p4total.Z() << ", " << p4total.E() << ")"
718  << std::endl;
719  abort();
720  }
721  } // found proton, pi+, pi-
722  } // end of select_type 6
723  //___________________________________________________________________________
724 
725  // Select p pi+ pi- pi0 events
726  else if (select_type == 7) {
727  bool foundProton = false;
728  bool foundPip = false;
729  bool foundPim = false;
730  bool foundPi0 = false;
731  int idProton = -999;
732  int idPip = -999;
733  int idPim = -999;
734  int idPi0 = -999;
735  TLorentzVector p4photon_init;
736  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
737  TLorentzVector p4proton;
738  TLorentzVector p4pip;
739  TLorentzVector p4pim;
740  TLorentzVector p4pi0;
741  TLorentzVector p4diff;
742  if (debug)
743  std::cout << "---------------------------------"
744  "---------------------------------" << std::endl;
745 
746  hddm_s::ReactionList res = record.getReactions();
747  hddm_s::ReactionList::iterator riter;
748  for (riter = res.begin(); riter != res.end(); ++riter) {
749  // Get initial photon
750  p4photon_init.SetXYZT(riter->getBeam().getMomentum().getPx(),
751  riter->getBeam().getMomentum().getPy(),
752  riter->getBeam().getMomentum().getPz(),
753  riter->getBeam().getMomentum().getE());
754  if (debug)
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() << ")"
759  << std::endl;
760 
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));
777  if (! finite(mass))
778  mass = 0.0;
779  if (debug) {
780  std::cout << "\t\t\t--- new particle ---" << std::endl;
781  std::cout << "\t\t\ttype = " << type << " "
782  << ParticleType((Particle_t)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;
789  }
790 
791  // Find proton
792  if (type == Proton) {
793  foundProton = true;
794  idProton = id;
795  p4proton.SetXYZT(px,py,pz,E);
796  }
797 
798  // Find pi+
799  if (type == PiPlus) {
800  foundPip = true;
801  idPip = id;
802  p4pip.SetXYZT(px,py,pz,E);
803  }
804 
805  // Find pi-
806  if (type == PiMinus) {
807  foundPim = true;
808  idPim = id;
809  p4pim.SetXYZT(px,py,pz,E);
810  }
811 
812  // Find pi0
813  if (type == Pi0) {
814  foundPi0 = true;
815  idPi0 = id;
816  p4pi0.SetXYZT(px,py,pz,E);
817  }
818 
819  } // end of loop over products
820  } // end of having products and origin
821  } // end of loop over vertices
822  } // end of loop over reaction
823 
824  // If we found a proton decaying from Lambda,
825  // save this event
826  if (foundProton && foundPip && foundPim && foundPi0) {
827 
828  // make sure that total 4-mom of p, pi+, pi- are close to
829  // initial state
830  p4diff = p4photon_init + p4proton_init - p4proton -
831  p4pip - p4pim - p4pi0;
832 
833  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
834  if (p4diff.P() < 0.001 && p4diff.M() < 0.0001)
835  return true;
836  } // found proton, pi+, pi-
837  } // end of select_type 7
838  //__________________________________________________________________________________________________
839 
840  // default is to return false
841  return false;
842 }
TRandom2 * rndm
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
const double alpha
TH1D * py[NCHANNELS]
static double ParticleMass(Particle_t p)
double sqrt(double)
bool selectEvent_s(int select_type, hddm_s::HDDM &record, int nevents, bool debug)
bool debug
#define SQR(X)
Particle_t
Definition: particleType.h:12