Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
selectEvents_r.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 <HDDM/hddm_r.hpp>
7 #include "particleType.h"
8 
9 using namespace std;
10 using std::distance;
11 
12 //-----------------------------------------------------------------
13 // selectEvent_s
14 //-----------------------------------------------------------------
15 bool selectEvent_r(int select_type, hddm_r::HDDM &record, int nevents, bool debug){
16 
17  // Select Lambda -> p pi- events
18  if(select_type==1){
19 
20  bool foundLambda = false;
21  int idLambda = -999;
22  int idProton = -999;
23  bool foundProton = false;
24  if(debug) cout << "event: " << nevents << "----------------------------------------------------------------" << endl;
25 
26  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
27  hddm_r::ReactionList reactions = re.getReactions();
28  // Loop over ReactionList
29  hddm_r::ReactionList::iterator iter;
30  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
31  int ireaction = std::distance(reactions.begin(), iter);
32 
33  hddm_r::VertexList vertices = iter->getVertices();
34  hddm_r::VertexList::iterator iter_vertex;
35  // Loop over VertexList
36  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
37  int ivertex = std::distance(vertices.begin(), iter_vertex);
38 
39  hddm_r::ProductList products = iter_vertex->getProducts();
40  hddm_r::ProductList::iterator iter_product;
41  // Loop over ProductList
42  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
43  int iparticle = std::distance(products.begin(), iter_product);
44 
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()
48  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
49  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
50  << endl;
51 
52  // Find Lambda
53  if(iter_product->getPdgtype() == PDGtype(Lambda)){
54  foundLambda = true;
55  idLambda = iter_product->getId();
56  }
57 
58  // Find decay products of Lambda
59  if(foundLambda){
60  // if parent was Lambda
61  if(iter_product->getParentId() == idLambda){
62  // find what type it is
63  if(iter_product->getPdgtype() == PDGtype(Proton)){
64  foundProton = true;
65  idProton = iter_product->getId();
66  }
67  }
68  }
69 
70  // If we found a proton decaying from Lambda,
71  // save this event
72  if(foundProton){
73  return true;
74  }
75 
76  } // end of loop over ProductList
77  } // end of loop over VertexList
78  } // end of loop over ReactiionList
79  }
80  //__________________________________________________________________________________________________
81 
82  // select Lambda -> p pi-, eta -> gamma gamma events
83  else if(select_type==2){
84 
85  bool foundLambda = false;
86  int idLambda = -999;
87  int idProton = -999;
88  bool foundProton = false;
89 
90  bool foundEta = false;
91  int idEta = -999;
92  int idPhoton[2];
93  int nPhoton = 0;
94 
95  if(debug) cout << "event: " << nevents << "----------------------------------------------------------------" << endl;
96 
97  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
98  hddm_r::ReactionList reactions = re.getReactions();
99  // Loop over ReactionList
100  hddm_r::ReactionList::iterator iter;
101  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
102  int ireaction = std::distance(reactions.begin(), iter);
103 
104  hddm_r::VertexList vertices = iter->getVertices();
105  hddm_r::VertexList::iterator iter_vertex;
106  // Loop over VertexList
107  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
108  int ivertex = std::distance(vertices.begin(), iter_vertex);
109 
110  hddm_r::ProductList products = iter_vertex->getProducts();
111  hddm_r::ProductList::iterator iter_product;
112  // Loop over ProductList
113  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
114  int iparticle = std::distance(products.begin(), iter_product);
115 
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()
119  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
120  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
121  << endl;
122 
123  // Find Lambda
124  if(iter_product->getPdgtype() == PDGtype(Lambda)){
125  foundLambda = true;
126  idLambda = iter_product->getId();
127  }
128 
129  // Find decay products of Lambda
130  if(foundLambda){
131  // if parent was Lambda
132  if(iter_product->getParentId() == idLambda){
133  // find what type it is
134  if(iter_product->getPdgtype() == PDGtype(Proton)){
135  foundProton = true;
136  idProton = iter_product->getId();
137  }
138  }
139  }
140 
141  // Find Eta
142  if(iter_product->getPdgtype() == PDGtype(Eta)){
143  foundEta = true;
144  idEta = iter_product->getId();
145  }
146 
147  // Find decay products of Eta
148  if(foundEta){
149  // if parent was Eta
150  if(iter_product->getParentId() == idEta){
151  // find what type it is
152  if(iter_product->getPdgtype() == PDGtype(Gamma)){
153  idPhoton[nPhoton] = iter_product->getId();
154  nPhoton++;
155  }
156  }
157  }
158 
159  // If we found a proton decaying from Lambda,
160  // save this event
161  if(foundProton && nPhoton==2){
162  return true;
163  }
164 
165  } // end of loop over ProductList
166  } // end of loop over VertexList
167  } // end of loop over ReactiionList
168  }
169  //__________________________________________________________________________________________________
170 
171  // Select Lambda -> p pi- events, decay length is less than VERTEXDIFFMAX
172  else if(select_type==3){
173  const double VERTEXDIFFMAX = 10.0;
174 
175  bool foundLambda = false;
176  int idLambda = -999;
177  int idProton = -999;
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;
184  }
185  double vertexdiff = 999;
186 
187  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
188  hddm_r::ReactionList reactions = re.getReactions();
189  // Loop over ReactionList
190  hddm_r::ReactionList::iterator iter;
191  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
192  int ireaction = std::distance(reactions.begin(), iter);
193 
194  hddm_r::VertexList vertices = iter->getVertices();
195  hddm_r::VertexList::iterator iter_vertex;
196  // Loop over VertexList
197  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
198  int ivertex = std::distance(vertices.begin(), iter_vertex);
199 
200  hddm_r::ProductList products = iter_vertex->getProducts();
201  hddm_r::ProductList::iterator iter_product;
202  // Loop over ProductList
203  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
204  int iparticle = std::distance(products.begin(), iter_product);
205 
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()
209  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
210  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
211  << endl;
212 
213  // Find Lambda
214  if(iter_product->getPdgtype() == PDGtype(Lambda)){
215  foundLambda = true;
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();
220  }
221 
222  // Find decay products of Lambda
223  if(foundLambda){
224  // if parent was Lambda
225  if(iter_product->getParentId() == idLambda){
226  // find what type it is
227  if(iter_product->getPdgtype() == PDGtype(Proton)){
228  foundProton = true;
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.));
234  }
235  }
236  }
237 
238  // If we found a proton decaying from Lambda,
239  // and the vertex difference is less than VERTEXDIFFMAX,
240  // save this event
241  if(foundProton && vertexdiff < VERTEXDIFFMAX){
242  return true;
243  }
244 
245  } // end of loop over ProductList
246  } // end of loop over VertexList
247  } // end of loop over ReactiionList
248  } // end of select_type 3
249  //__________________________________________________________________________________________________
250 
251  // Select p pi+ pi- events
252  else if(select_type==5){
253  bool foundProton = false;
254  bool foundPip = false;
255  bool foundPim = false;
256  int idProton = -999;
257  int idPip = -999;
258  int idPim = -999;
259  TLorentzVector p4photon_init;
260  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
261  TLorentzVector p4proton;
262  TLorentzVector p4pip;
263  TLorentzVector p4pim;
264  TLorentzVector p4diff;
265 
266  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
267  hddm_r::ReactionList reactions = re.getReactions();
268  // Loop over ReactionList
269  hddm_r::ReactionList::iterator iter;
270  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
271  int ireaction = std::distance(reactions.begin(), iter);
272 
273  hddm_r::VertexList vertices = iter->getVertices();
274  hddm_r::VertexList::iterator iter_vertex;
275  // Loop over VertexList
276  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
277  int ivertex = std::distance(vertices.begin(), iter_vertex);
278 
279  hddm_r::ProductList products = iter_vertex->getProducts();
280  hddm_r::ProductList::iterator iter_product;
281  // Loop over ProductList
282  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
283  int iparticle = std::distance(products.begin(), iter_product);
284 
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()
288  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
289  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
290  << " Ebeam: " << iter_vertex->getOrigin().getEbeam()
291  << endl;
292 
293  p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
294 
295  // Find proton
296  if(iter_product->getPdgtype() == PDGtype(Proton)){
297  foundProton = true;
298  idProton = iter_product->getId();
299  p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
300  }
301 
302  // Find pi+
303  if(iter_product->getPdgtype() == PDGtype(PiPlus)){
304  foundPip = true;
305  idPip = iter_product->getId();
306  p4pip.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
307  }
308 
309  // Find pi-
310  if(iter_product->getPdgtype() == PDGtype(PiMinus)){
311  foundPim = true;
312  idPim = iter_product->getId();
313  p4pim.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
314  }
315  } // end of loop over ProductList
316  } // end of loop over VertexList
317  } // end of loop over ReactiionList
318 
319  // If we found a proton decaying from Lambda,
320  // save this event
321  if(foundProton && foundPip && foundPim){
322 
323  // make sure that total 4-mom of p, pi+, pi- are close to
324  // initial state
325  p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim;
326 
327  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
328  if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
329  return true;
330  }
331  } // found proton, pi+, pi-
332 
333  } // end of select_type 5
334  //__________________________________________________________________________________________________
335 
336  // Select K+ K+ Xi- -> K+ K+ p pi- pi- events
337  else if(select_type==6){
338  Int_t nKp = 0;
339  Int_t nXiMinus = 0;
340  Int_t nLambda = 0;
341  Int_t nProton = 0;
342  Int_t nPimFromLambda = 0;
343  Int_t nPimFromXiMinus = 0;
344  // need intermediate IDs
345  int idXiMinus = -999;
346  int idLambda = -999;
347  TLorentzVector p4photon_init;
348  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
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;
358 
359  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
360  hddm_r::ReactionList reactions = re.getReactions();
361  // Loop over ReactionList
362  hddm_r::ReactionList::iterator iter;
363  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
364  int ireaction = std::distance(reactions.begin(), iter);
365 
366  hddm_r::VertexList vertices = iter->getVertices();
367  hddm_r::VertexList::iterator iter_vertex;
368  // Loop over VertexList
369  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
370  int ivertex = std::distance(vertices.begin(), iter_vertex);
371 
372  hddm_r::ProductList products = iter_vertex->getProducts();
373  hddm_r::ProductList::iterator iter_product;
374  // Loop over ProductList
375  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
376  int iparticle = std::distance(products.begin(), iter_product);
377 
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()
381  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
382  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
383  << " Ebeam: " << iter_vertex->getOrigin().getEbeam()
384  << endl;
385 
386  p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
387 
388  // Find first K+
389  if(nKp==0 && iter_product->getPdgtype() == PDGtype(KPlus)){
390  nKp++;
391  p4kp1.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
392  if(debug)
393  cout << "K+ 1 : (" << p4kp1.X() << ", " << p4kp1.Y() << ", " << p4kp1.Z() << ", " << p4kp1.M() << ")" << endl;
394 
395  // Need this break since without it we would count the
396  // the first K+ as the second K+ too
397  continue;
398  }
399 
400  // Find second K+
401  if(nKp==1 && iter_product->getPdgtype() == PDGtype(KPlus)){
402  nKp++;
403  p4kp1.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
404  if(debug)
405  cout << "K+ 2 : (" << p4kp2.X() << ", " << p4kp2.Y() << ", " << p4kp2.Z() << ", " << p4kp2.M() << ")" << endl;
406  }
407 
408  // Find Xi-
409  if(iter_product->getPdgtype() == PDGtype(XiMinus)){
410  nXiMinus++;
411  idXiMinus = iter_product->getId();
412  p4XiMinus.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
413  }
414 
415  // Find pi- from Xi-
416  if(iter_product->getPdgtype() == PDGtype(PiMinus) && iter_product->getParentId() == idXiMinus){
417  nPimFromXiMinus++;
418  p4pimFromXiMinus.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
419  }
420 
421  // Find Lambda
422  if(iter_product->getPdgtype() == PDGtype(Lambda) && iter_product->getParentId() == idXiMinus){
423  nLambda++;
424  idLambda = iter_product->getId();
425  p4Lambda.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
426  }
427 
428  // Find proton
429  if(iter_product->getPdgtype() == PDGtype(Proton) && iter_product->getParentId() == idLambda){
430  nProton++;
431  p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
432  }
433 
434  // Find pi- from Lambda
435  if(iter_product->getPdgtype() == PDGtype(PiMinus) && iter_product->getParentId() == idLambda){
436  nPimFromLambda++;
437  p4pimFromLambda.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
438  }
439  } // end of loop over ProductList
440  } // end of loop over VertexList
441  } // end of loop over ReactiionList
442 
443  // If we n a proton decaying from Lambda,
444  // save this event
445  if(nKp==2 && nXiMinus==1 && nLambda==1 && nProton==1 && nPimFromXiMinus==1 && nPimFromLambda==1){
446 
447  // make sure that total 4-mom of p, pi+, pi- are close to
448  // initial state
449  p4diff = p4photon_init + p4proton_init - p4kp1 - p4kp2 - p4proton - p4pimFromXiMinus - p4pimFromLambda;
450  p4total = p4kp1 + p4kp2 + p4proton + p4pimFromXiMinus + p4pimFromLambda;
451 
452  if(debug){
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;
462  }
463 
464  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
465  if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
466  return true;
467  }else{
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;
472  abort();
473  }
474  } // found proton, pi+, pi-
475 
476  } // end of select_type 6
477  //__________________________________________________________________________________________________
478 
479  // Select p pi+ pi- pi0 events
480  else if(select_type==7){
481  bool foundProton = false;
482  bool foundPip = false;
483  bool foundPim = false;
484  bool foundPi0 = false;
485  int idProton = -999;
486  int idPip = -999;
487  int idPim = -999;
488  int idPi0 = -999;
489  TLorentzVector p4photon_init;
490  TLorentzVector p4proton_init(0,0,0,ParticleMass(Proton));
491  TLorentzVector p4proton;
492  TLorentzVector p4pip;
493  TLorentzVector p4pim;
494  TLorentzVector p4pi0;
495  TLorentzVector p4diff;
496 
497  hddm_r::ReconstructedPhysicsEvent &re = record.getReconstructedPhysicsEvent();
498  hddm_r::ReactionList reactions = re.getReactions();
499  // Loop over ReactionList
500  hddm_r::ReactionList::iterator iter;
501  for (iter = reactions.begin(); iter != reactions.end(); ++iter) {
502  int ireaction = std::distance(reactions.begin(), iter);
503 
504  hddm_r::VertexList vertices = iter->getVertices();
505  hddm_r::VertexList::iterator iter_vertex;
506  // Loop over VertexList
507  for (iter_vertex = vertices.begin(); iter_vertex != vertices.end(); ++iter_vertex){
508  int ivertex = std::distance(vertices.begin(), iter_vertex);
509 
510  hddm_r::ProductList products = iter_vertex->getProducts();
511  hddm_r::ProductList::iterator iter_product;
512  // Loop over ProductList
513  for (iter_product = products.begin(); iter_product != products.end(); ++iter_product){
514  int iparticle = std::distance(products.begin(), iter_product);
515 
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()
519  << " type: " << ParticleType((Particle_t)PDGtoPType((Particle_t)iter_product->getPdgtype()))
520  << " vertex: (" << iter_vertex->getOrigin().getVx() << ", " << iter_vertex->getOrigin().getVy() << ", " << iter_vertex->getOrigin().getVz() << ")"
521  << " Ebeam: " << iter_vertex->getOrigin().getEbeam()
522  << endl;
523 
524  p4photon_init.SetXYZT(0,0,iter_vertex->getOrigin().getEbeam(),iter_vertex->getOrigin().getEbeam());
525 
526  // Find proton
527  if(iter_product->getPdgtype() == PDGtype(Proton)){
528  foundProton = true;
529  idProton = iter_product->getId();
530  p4proton.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
531  }
532 
533  // Find pi+
534  if(iter_product->getPdgtype() == PDGtype(PiPlus)){
535  foundPip = true;
536  idPip = iter_product->getId();
537  p4pip.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
538  }
539 
540  // Find pi-
541  if(iter_product->getPdgtype() == PDGtype(PiMinus)){
542  foundPim = true;
543  idPim = iter_product->getId();
544  p4pim.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
545  }
546 
547  // Find pi-
548  if(iter_product->getPdgtype() == PDGtype(Pi0)){
549  foundPi0 = true;
550  idPi0 = iter_product->getId();
551  p4pi0.SetXYZT(iter_product->getMomentum().getPx(),iter_product->getMomentum().getPy(),iter_product->getMomentum().getPz(),iter_product->getMomentum().getE());
552  }
553 
554  } // end of loop over ProductList
555  } // end of loop over VertexList
556  } // end of loop over ReactiionList
557 
558  // If we found a proton decaying from Lambda,
559  // save this event
560  if(foundProton && foundPip && foundPim && foundPi0){
561 
562  // make sure that total 4-mom of p, pi+, pi- are close to
563  // initial state
564  p4diff = p4photon_init + p4proton_init - p4proton - p4pip - p4pim - p4pi0;
565 
566  // Make sure that remaining |p| < 1 MeV, M < 0.1 MeV
567  if(p4diff.P() < 0.001 && p4diff.M() < 0.0001){
568  return true;
569  }
570  } // found proton, pi+, pi-
571 
572  } // end of select_type 7
573  //__________________________________________________________________________________________________
574 
575  // default is to return false
576  return false;
577 }
hddm_select_events.h
PDGtoPType
static Particle_t PDGtoPType(int locPDG_PID)
Definition: particleType.h:1246
selectEvent_r
bool selectEvent_r(int select_type, hddm_r::HDDM &record, int nevents, bool debug)
Definition: selectEvents_r.cc:15
ParticleType
static char * ParticleType(Particle_t p)
Definition: particleType.h:142
PDGtype
static int PDGtype(Particle_t p)
Definition: particleType.h:1167
ParticleMass
static double ParticleMass(Particle_t p)
Definition: particleType.h:1002
sqrt
double sqrt(double)
debug
bool debug
Definition: hddm_select_events.cc:17
Particle_t
Particle_t
Definition: particleType.h:12