Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHSartre.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHSartre.cc
1 #include "PHSartre.h"
2 #include "PHSartreGenTrigger.h"
3 
4 #include <sartre/Enumerations.h> // for incoherent
5 #include <sartre/Event.h> // for Particle, Event
6 #include <sartre/EventGeneratorSettings.h> // for EventGeneratorSettings
7 
8 
9 #include <phhepmc/PHHepMCGenHelper.h> // for PHHepMCGenHelper
10 
12 #include <fun4all/SubsysReco.h> // for SubsysReco
13 
14 #include <sartre/Sartre.h>
15 
16 #include <TGenPhaseSpace.h>
17 #include <TLorentzVector.h> // for TLorentzVector
18 #include <TParticlePDG.h> // for TParticlePDG
19 
20 #include <CLHEP/Vector/LorentzVector.h>
21 
22 #include <HepMC/GenEvent.h>
23 #include <HepMC/GenParticle.h> // for GenParticle
24 #include <HepMC/GenVertex.h> // for GenVertex
25 #include <HepMC/PdfInfo.h> // for PdfInfo
26 #include <HepMC/SimpleVector.h> // for FourVector
27 #include <HepMC/Units.h> // for GEV, MM
28 
29 #include <gsl/gsl_rng.h> // for gsl_rng_uniform
30 
31 #include <cfloat> // for FLT_EPSILON
32 #include <cmath> // for M_PI
33 #include <cstdlib> // for getenv
34 #include <iostream> // for operator<<, endl, basic_o...
35 #include <memory> // for allocator_traits<>::value...
36 
37 class PHHepMCGenEvent;
38 
39 using namespace std;
40 
42  : SubsysReco(name)
43 {
44  char *charPath = getenv("SARTRE_DIR");
45  if (!charPath)
46  {
47  cout << "PHSartre::Could not find $SARTRE_DIR path!" << endl;
48  return;
49  }
50 
51  //
52  // Create the generator and initialize it.
53  // Once initialized you cannot (should not) change
54  // the settings w/o re-initialing sartre.
55  //
56  _sartre = new Sartre();
57 
58  PHHepMCGenHelper::set_embedding_id(1); // default embedding ID to 1
59 }
60 
62 {
63  delete _sartre;
64 }
65 
67 {
68  if (!_configFile.empty())
69  {
70  bool ok = _sartre->init(_configFile);
71  if (!ok)
72  {
73  cerr << "Initialization of sartre failed." << endl;
75  }
76  }
77  else
78  {
79  cerr << "Sartre configuration file must be specified" << endl;
81  }
82 
83  settings = _sartre->runSettings();
84  settings->list();
85 
86  create_node_tree(topNode);
87 
88  // event numbering will start from 1
89  _eventcount = 0;
90 
91  decay = new TGenPhaseSpace(); // for VM decays
92  daughterID = settings->userInt();
93  if (daughterID && (settings->vectorMesonId() != 22))
94  {
95  doPerformDecay = true;
96  daughterMasses[0] = settings->lookupPDG(daughterID)->Mass();
97  daughterMasses[1] = settings->lookupPDG(-daughterID)->Mass();
98  cout << "PHSartre: "
99  << "Will decay vector meson: ";
100  cout << "PHSartre: " << settings->lookupPDG(settings->vectorMesonId())->GetName();
101  cout << "PHSartre: "
102  << " -> ";
103  cout << "PHSartre: " << settings->lookupPDG(daughterID)->GetName();
104  cout << " ";
105  cout << "PHSartre: " << settings->lookupPDG(-daughterID)->GetName();
106  cout << endl;
107  }
108 
110 }
111 
113 {
114  if (Verbosity() > 1) cout << "PHSartre::End - I'm here!" << endl;
115 
116  cout << "PHSartre: "
117  << " Total cross-section: " << _sartre->totalCrossSection() << " nb" << endl;
118  _sartre->listStatus();
119 
120  cout << " *------- Begin PHSARTRE Trigger Statistics ----------------------"
121  << "-------------------------------------------------* " << endl;
122  cout << " | "
123  << " | " << endl;
124  cout << " PHSartre::End - " << _eventcount
125  << " events passed trigger" << endl;
126  cout << " Fraction passed: " << _eventcount
127  << "/" << _gencount
128  << " = " << _eventcount / float(_gencount) << endl;
129  cout << " *------- End PHSARTRE Trigger Statistics ------------------------"
130  << "-------------------------------------------------* " << endl;
131 
133 }
134 
135 //-* print pythia config info
137 {
138  settings->list();
139 }
140 
142 {
143  if (Verbosity() > 1) cout << "PHSartre::process_event - event: " << _eventcount << endl;
144 
145  bool passedTrigger = false;
146  Event *event = nullptr;
147 
148  TLorentzVector *eIn = nullptr;
149  TLorentzVector *pIn = nullptr;
150  TLorentzVector *eOut = nullptr;
151  TLorentzVector *gamma = nullptr;
152  TLorentzVector *vm = nullptr;
153  TLorentzVector *PomOut = nullptr;
154  TLorentzVector *pOut = nullptr;
155  TLorentzVector *vmDecay1 = nullptr;
156  TLorentzVector *vmDecay2 = nullptr;
157  unsigned int preVMDecaySize = 0;
158 
159  while (!passedTrigger)
160  {
161  ++_gencount;
162 
163  // Generate a Sartre event
164  event = _sartre->generateEvent();
165 
166  //
167  // If Sartre is run in UPC mode, half of the events needs to be
168  // rotated around and axis perpendicular to z:
169  // (only for symmetric events)
170  //
171  if (settings->UPC() and settings->A() == settings->UPCA())
172  {
174  }
175 
176  // for sPHENIX/RHIC p+Au
177  // (see comments in ReverseBeams)
178  // reverse when the proton emits the virtual photon
179 
180  if (settings->UPC() and settings->A() == 197)
181  {
183  }
184 
185  // Set pointers to the parts of the event we will need:
186 
187  eIn = &event->particles[0].p;
188  pIn = &event->particles[1].p;
189  eOut = &event->particles[2].p;
190  gamma = &event->particles[3].p;
191  vm = &event->particles[4].p;
192  PomOut = &event->particles[5].p;
193  pOut = &event->particles[6].p;
194 
195  // To allow the triggering to work properly, we need to decay the vector meson here
196 
197  preVMDecaySize = event->particles.size();
198 
199  if (doPerformDecay)
200  {
201  if (decay->SetDecay(*vm, 2, daughterMasses))
202  {
203  double weight = decay->Generate(); // weight is always 1 here
204  if ((weight - 1) > FLT_EPSILON)
205  {
206  cout << "PHSartre: Warning decay weight != 1, weight = " << weight << endl;
207  }
208  TLorentzVector *vmDaughter1 = decay->GetDecay(0);
209  TLorentzVector *vmDaughter2 = decay->GetDecay(1);
210 
211  event->particles[4].status = 2; // set VM status
212 
213  Particle vmDC1;
214  vmDC1.index = event->particles.size();
215  vmDC1.pdgId = daughterID;
216  vmDC1.status = 1; // final state
217  vmDC1.p = *vmDaughter1;
218  vmDC1.parents.push_back(4);
219  event->particles.push_back(vmDC1);
220  vmDecay1 = &event->particles[event->particles.size() - 1].p;
221 
222  Particle vmDC2;
223  vmDC2.index = event->particles.size();
224  vmDC2.pdgId = -daughterID;
225  vmDC2.status = 1; // final state
226  vmDC2.p = *vmDaughter2;
227  vmDC2.parents.push_back(4);
228  event->particles.push_back(vmDC2);
229  vmDecay2 = &event->particles[event->particles.size() - 1].p;
230  }
231  else
232  {
233  cout << "PHSartre: WARNING: Kinematics of Vector Meson does not allow decay!" << endl;
234  }
235  }
236 
237  // test trigger logic
238 
239  bool andScoreKeeper = true;
240  if (Verbosity() > 2)
241  {
242  cout << "PHSartre::process_event - triggersize: " << _registeredTriggers.size() << endl;
243  }
244 
245  for (unsigned int tr = 0; tr < _registeredTriggers.size(); tr++)
246  {
247  bool trigResult = _registeredTriggers[tr]->Apply(event);
248 
249  if (Verbosity() > 2)
250  {
251  cout << "PHSartre::process_event trigger: "
252  << _registeredTriggers[tr]->GetName() << " " << trigResult << endl;
253  }
254 
255  if (_triggersOR && trigResult)
256  {
257  passedTrigger = true;
258  break;
259  }
260  else if (_triggersAND)
261  {
262  andScoreKeeper &= trigResult;
263  }
264 
265  if (Verbosity() > 2 && !passedTrigger)
266  {
267  cout << "PHSartre::process_event - failed trigger: "
268  << _registeredTriggers[tr]->GetName() << endl;
269  }
270  }
271 
272  if ((andScoreKeeper && _triggersAND) || (_registeredTriggers.size() == 0))
273  {
274  passedTrigger = true;
275  }
276  }
277 
278  // fill HepMC object with event
279 
280  HepMC::GenEvent *genevent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
281 
282  // add some information to the event
283  genevent->set_event_number(_eventcount);
284 
285  // Set the PDF information
286  HepMC::PdfInfo pdfinfo;
287  pdfinfo.set_scalePDF(event->Q2);
288  genevent->set_pdf_info(pdfinfo);
289 
290  // We would also like to save:
291  //
292  // event->t;
293  // event->x;
294  // event->y;
295  // event->s;
296  // event->W;
297  // event->xpom;
298  // (event->polarization == transverse ? 0 : 1);
299  // (event->diffractiveMode == coherent ? 0 : 1);
300  //
301  // but there doesn't seem to be a good place to do so
302  // within the HepMC event information?
303  //
304  // t, W and Q^2 form a minial set of good variables for diffractive events
305  // Maybe what I do is record the input particles to the event at the HepMC
306  // vertices and reconstruct the kinematics from there?
307 
308  // Create HepMC vertices and add final state particles to them
309 
310  // First, the emitter(electron)-virtual photon vertex:
311 
312  HepMC::GenVertex *egammavtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
313  genevent->add_vertex(egammavtx);
314 
315  egammavtx->add_particle_in(
316  new HepMC::GenParticle(CLHEP::HepLorentzVector(eIn->Px(),
317  eIn->Py(),
318  eIn->Pz(),
319  eIn->E()),
320  event->particles[0].pdgId,
321  3));
322 
323  HepMC::GenParticle *hgamma = new HepMC::GenParticle(CLHEP::HepLorentzVector(gamma->Px(),
324  gamma->Py(),
325  gamma->Pz(),
326  gamma->E()),
327  event->particles[3].pdgId,
328  3);
329 
330  egammavtx->add_particle_out(hgamma);
331 
332  egammavtx->add_particle_out(
333  new HepMC::GenParticle(CLHEP::HepLorentzVector(eOut->Px(),
334  eOut->Py(),
335  eOut->Pz(),
336  eOut->E()),
337  event->particles[2].pdgId,
338  1));
339 
340  // Next, the hadron-pomeron vertex:
341 
342  HepMC::GenVertex *ppomvtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
343  genevent->add_vertex(ppomvtx);
344 
345  ppomvtx->add_particle_in(
346  new HepMC::GenParticle(CLHEP::HepLorentzVector(pIn->Px(),
347  pIn->Py(),
348  pIn->Pz(),
349  pIn->E()),
350  event->particles[1].pdgId,
351  3));
352 
353  HepMC::GenParticle *hPomOut = new HepMC::GenParticle(CLHEP::HepLorentzVector(PomOut->Px(),
354  PomOut->Py(),
355  PomOut->Pz(),
356  PomOut->E()),
357  event->particles[5].pdgId,
358  3);
359 
360  ppomvtx->add_particle_out(hPomOut);
361 
362  // If this is a nuclear breakup, add in the nuclear fragments
363  // Otherwise, add in the outgoing hadron
364 
365  //If the event is incoherent, and nuclear breakup is enabled, fill the remnants to the tree
366  if (settings->enableNuclearBreakup() and event->diffractiveMode == incoherent)
367  {
368  for (unsigned int iParticle = 7; iParticle < preVMDecaySize; iParticle++)
369  {
370  if (event->particles[iParticle].status == 1)
371  { // Final-state particle
372  const Particle &particle = event->particles[iParticle];
373  ppomvtx->add_particle_out(
374  new HepMC::GenParticle(CLHEP::HepLorentzVector(particle.p.Px(),
375  particle.p.Py(),
376  particle.p.Pz(),
377  particle.p.E()),
378  particle.pdgId,
379  1));
380  }
381  }
382  }
383  else
384  {
385  ppomvtx->add_particle_out(
386  new HepMC::GenParticle(CLHEP::HepLorentzVector(pOut->Px(),
387  pOut->Py(),
388  pOut->Pz(),
389  pOut->E()),
390  event->particles[6].pdgId,
391  1));
392  }
393 
394  // The Pomeron-Photon vertex
395 
396  HepMC::GenVertex *gammapomvtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
397  genevent->add_vertex(gammapomvtx);
398 
399  gammapomvtx->add_particle_in(hgamma);
400  gammapomvtx->add_particle_in(hPomOut);
401 
402  int isVMFinal = 1;
403  if (doPerformDecay) isVMFinal = 2;
404 
405  HepMC::GenParticle *hvm = new HepMC::GenParticle(CLHEP::HepLorentzVector(vm->Px(),
406  vm->Py(),
407  vm->Pz(),
408  vm->E()),
409  event->particles[4].pdgId,
410  isVMFinal);
411 
412  gammapomvtx->add_particle_out(hvm);
413 
414  // Add the VM decay to the event
415 
416  if (doPerformDecay)
417  {
418  if (vmDecay1 && vmDecay2)
419  {
420  HepMC::GenVertex *fvtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
421  genevent->add_vertex(fvtx);
422 
423  fvtx->add_particle_in(hvm);
424 
425  fvtx->add_particle_out(
426  new HepMC::GenParticle(CLHEP::HepLorentzVector(vmDecay1->Px(),
427  vmDecay1->Py(),
428  vmDecay1->Pz(),
429  vmDecay1->E()),
430  daughterID,
431  1));
432  fvtx->add_particle_out(
433  new HepMC::GenParticle(CLHEP::HepLorentzVector(vmDecay2->Px(),
434  vmDecay2->Py(),
435  vmDecay2->Pz(),
436  vmDecay2->E()),
437  -daughterID,
438  1));
439  }
440  else
441  {
442  cout << "PHSartre: WARNING: Kinematics of Vector Meson does not allow decay!" << endl;
443  }
444  }
445 
446  // pass HepMC to PHNode
447 
448  PHHepMCGenEvent *success = PHHepMCGenHelper::insert_event(genevent);
449  if (!success)
450  {
451  cout << "PHSartre::process_event - Failed to add event to HepMC record!" << endl;
453  }
454 
455  // print outs
456 
457  if (Verbosity() > 2) cout << "PHSartre::process_event - FINISHED WHOLE EVENT" << endl;
458 
459  ++_eventcount;
461 }
462 
463 
464 
466 {
468 }
469 
471 {
472  if (Verbosity() > 1) cout << "PHSartre::registerTrigger - trigger " << theTrigger->GetName() << " registered" << endl;
473  _registeredTriggers.push_back(theTrigger);
474 }
475 
476 // UPC only
478 {
479  if (gsl_rng_uniform(PHHepMCGenHelper::get_random_generator()) > 0.5)
480  {
481  for (unsigned int i = 0; i < myEvent->particles.size(); i++)
482  myEvent->particles.at(i).p.RotateX(M_PI);
483  }
484 }
485 
486 // Used to rotate into the sPHENIX/RHIC frame
487 // The photon emitting beam is always pz<0, and in sPHENIX this will be
488 // the ion direction. So what you want to run are two files with:
489 // A=1 UPCA=197 (no reversal, Au ion emits photon)
490 // A=197 UPCA=1 (reversal required, proton emits photon)
492 {
493  for (unsigned int i = 0; i < myEvent->particles.size(); i++)
494  myEvent->particles.at(i).p.RotateX(M_PI);
495 }