4 #include <sartre/Enumerations.h>
5 #include <sartre/Event.h>
6 #include <sartre/EventGeneratorSettings.h>
14 #include <sartre/Sartre.h>
16 #include <TGenPhaseSpace.h>
17 #include <TLorentzVector.h>
18 #include <TParticlePDG.h>
20 #include <CLHEP/Vector/LorentzVector.h>
22 #include <HepMC/GenEvent.h>
23 #include <HepMC/GenParticle.h>
24 #include <HepMC/GenVertex.h>
25 #include <HepMC/PdfInfo.h>
26 #include <HepMC/SimpleVector.h>
27 #include <HepMC/Units.h>
29 #include <gsl/gsl_rng.h>
44 char *charPath = getenv(
"SARTRE_DIR");
47 cout <<
"PHSartre::Could not find $SARTRE_DIR path!" << endl;
73 cerr <<
"Initialization of sartre failed." << endl;
79 cerr <<
"Sartre configuration file must be specified" << endl;
91 decay =
new TGenPhaseSpace();
99 <<
"Will decay vector meson: ";
100 cout <<
"PHSartre: " <<
settings->lookupPDG(
settings->vectorMesonId())->GetName();
114 if (
Verbosity() > 1) cout <<
"PHSartre::End - I'm here!" << endl;
117 <<
" Total cross-section: " <<
_sartre->totalCrossSection() <<
" nb" << endl;
120 cout <<
" *------- Begin PHSARTRE Trigger Statistics ----------------------"
121 <<
"-------------------------------------------------* " << endl;
125 <<
" events passed trigger" << endl;
129 cout <<
" *------- End PHSARTRE Trigger Statistics ------------------------"
130 <<
"-------------------------------------------------* " << endl;
145 bool passedTrigger =
false;
146 Event *
event =
nullptr;
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;
159 while (!passedTrigger)
164 event =
_sartre->generateEvent();
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;
197 preVMDecaySize =
event->particles.size();
203 double weight =
decay->Generate();
204 if ((weight - 1) > FLT_EPSILON)
206 cout <<
"PHSartre: Warning decay weight != 1, weight = " << weight << endl;
208 TLorentzVector *vmDaughter1 =
decay->GetDecay(0);
209 TLorentzVector *vmDaughter2 =
decay->GetDecay(1);
211 event->particles[4].status = 2;
214 vmDC1.index =
event->particles.size();
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;
223 vmDC2.index =
event->particles.size();
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;
233 cout <<
"PHSartre: WARNING: Kinematics of Vector Meson does not allow decay!" << endl;
239 bool andScoreKeeper =
true;
251 cout <<
"PHSartre::process_event trigger: "
257 passedTrigger =
true;
262 andScoreKeeper &= trigResult;
267 cout <<
"PHSartre::process_event - failed trigger: "
274 passedTrigger =
true;
280 HepMC::GenEvent *
genevent =
new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
286 HepMC::PdfInfo pdfinfo;
287 pdfinfo.set_scalePDF(
event->Q2);
288 genevent->set_pdf_info(pdfinfo);
312 HepMC::GenVertex *egammavtx =
new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
313 genevent->add_vertex(egammavtx);
315 egammavtx->add_particle_in(
316 new HepMC::GenParticle(CLHEP::HepLorentzVector(eIn->Px(),
320 event->particles[0].pdgId,
323 HepMC::GenParticle *hgamma =
new HepMC::GenParticle(CLHEP::HepLorentzVector(gamma->Px(),
327 event->particles[3].pdgId,
330 egammavtx->add_particle_out(hgamma);
332 egammavtx->add_particle_out(
333 new HepMC::GenParticle(CLHEP::HepLorentzVector(eOut->Px(),
337 event->particles[2].pdgId,
342 HepMC::GenVertex *ppomvtx =
new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
343 genevent->add_vertex(ppomvtx);
345 ppomvtx->add_particle_in(
346 new HepMC::GenParticle(CLHEP::HepLorentzVector(pIn->Px(),
350 event->particles[1].pdgId,
353 HepMC::GenParticle *hPomOut =
new HepMC::GenParticle(CLHEP::HepLorentzVector(PomOut->Px(),
357 event->particles[5].pdgId,
360 ppomvtx->add_particle_out(hPomOut);
366 if (
settings->enableNuclearBreakup() and
event->diffractiveMode == incoherent)
368 for (
unsigned int iParticle = 7; iParticle < preVMDecaySize; iParticle++)
370 if (
event->particles[iParticle].status == 1)
372 const Particle &
particle =
event->particles[iParticle];
373 ppomvtx->add_particle_out(
374 new HepMC::GenParticle(CLHEP::HepLorentzVector(particle.p.Px(),
385 ppomvtx->add_particle_out(
386 new HepMC::GenParticle(CLHEP::HepLorentzVector(pOut->Px(),
390 event->particles[6].pdgId,
396 HepMC::GenVertex *gammapomvtx =
new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
397 genevent->add_vertex(gammapomvtx);
399 gammapomvtx->add_particle_in(hgamma);
400 gammapomvtx->add_particle_in(hPomOut);
405 HepMC::GenParticle *hvm =
new HepMC::GenParticle(CLHEP::HepLorentzVector(vm->Px(),
409 event->particles[4].pdgId,
412 gammapomvtx->add_particle_out(hvm);
418 if (vmDecay1 && vmDecay2)
420 HepMC::GenVertex *fvtx =
new HepMC::GenVertex(CLHEP::HepLorentzVector(0.0, 0.0, 0.0, 0.0));
421 genevent->add_vertex(fvtx);
423 fvtx->add_particle_in(hvm);
425 fvtx->add_particle_out(
426 new HepMC::GenParticle(CLHEP::HepLorentzVector(vmDecay1->Px(),
432 fvtx->add_particle_out(
433 new HepMC::GenParticle(CLHEP::HepLorentzVector(vmDecay2->Px(),
442 cout <<
"PHSartre: WARNING: Kinematics of Vector Meson does not allow decay!" << endl;
451 cout <<
"PHSartre::process_event - Failed to add event to HepMC record!" << endl;
457 if (
Verbosity() > 2) cout <<
"PHSartre::process_event - FINISHED WHOLE EVENT" << endl;
472 if (
Verbosity() > 1) cout <<
"PHSartre::registerTrigger - trigger " << theTrigger->
GetName() <<
" registered" << endl;
481 for (
unsigned int i = 0;
i < myEvent->particles.size();
i++)
482 myEvent->particles.at(
i).p.RotateX(M_PI);
493 for (
unsigned int i = 0;
i < myEvent->particles.size();
i++)
494 myEvent->particles.at(
i).p.RotateX(M_PI);