1+ #include "FairGenerator.h"
2+ #include "Generators/GeneratorPythia8.h"
3+ #include "Pythia8/Pythia.h"
4+ #include "TRandom.h"
5+
6+ R__ADD_INCLUDE_PATH ($O2DPG_MC_CONFIG_ROOT /MC /config /PWGDQ /EvtGen )
7+ #include "GeneratorEvtGen.C"
8+
9+ #include <string>
10+
11+ using namespace o2 ::eventgen ;
12+
13+ namespace o2
14+ {
15+ namespace eventgen
16+ {
17+
18+ class GeneratorPythia8HadronTriggeredWithGap : public o2 ::eventgen ::GeneratorPythia8 {
19+ public :
20+
21+ /// constructor
22+ GeneratorPythia8HadronTriggeredWithGap (int inputTriggerRatio = 5 ) {
23+
24+ mGeneratedEvents = 0 ;
25+ mInverseTriggerRatio = inputTriggerRatio ;
26+ // define minimum bias event generator
27+ auto seed = (gRandom -> TRandom ::GetSeed () % 900000000 );
28+ // main physics option for the min bias pythia events: SoftQCD:Inelastic
29+ TString pathconfigMB = gSystem -> ExpandPathName ("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/generator/pythia8_inel_triggerGap.cfg" );
30+ pythiaMBgen .readFile (pathconfigMB .Data ());
31+ pythiaMBgen .readString ("Random:setSeed on" );
32+ pythiaMBgen .readString ("Random:seed " + std ::to_string (seed ));
33+ mConfigMBdecays = "" ;
34+ mRapidityMin = -1. ;
35+ mRapidityMax = 1. ;
36+ mVerbose = false;
37+ }
38+
39+ /// Destructor
40+ ~GeneratorPythia8HadronTriggeredWithGap () = default ;
41+
42+ void addHadronPDGs (int pdg ) { mHadronsPDGs .push_back (pdg ); };
43+
44+ void setRapidityRange (double valMin , double valMax )
45+ {
46+ mRapidityMin = valMin ;
47+ mRapidityMax = valMax ;
48+ };
49+
50+ void setTriggerGap (int triggerGap ) {mInverseTriggerRatio = triggerGap ;}
51+
52+ void setConfigMBdecays (TString val ){mConfigMBdecays = val ;}
53+
54+ void setVerbose (bool val ) { mVerbose = val ; };
55+
56+ protected :
57+
58+ bool generateEvent () override {
59+ // reset event
60+ bool genOk = false;
61+ if (mGeneratedEvents % mInverseTriggerRatio == 0 ) {
62+ bool found = false;
63+ while (! (genOk && found )) {
64+ /// reset event
65+ mPythia .event .reset ();
66+ genOk = GeneratorPythia8 ::generateEvent ();
67+ // find the q-qbar or single hadron ancestor
68+ found = findHadrons (mPythia .event );
69+ }
70+ notifySubGenerator (1 );
71+ } else {
72+ /// reset event
73+ pythiaMBgen .event .reset ();
74+ while (!genOk ) {
75+ genOk = pythiaMBgen .next ();
76+ }
77+ mPythia .event = pythiaMBgen .event ;
78+ notifySubGenerator (0 );
79+ }
80+ mGeneratedEvents ++ ;
81+ if (mVerbose ) {
82+ mOutputEvent .list ();
83+ }
84+ return true;
85+ }
86+
87+ bool Init () override {
88+
89+ if (mConfigMBdecays .Contains ("cfg" )) {
90+ pythiaMBgen .readFile (mConfigMBdecays .Data ());
91+ }
92+ addSubGenerator (0 , "Minimum bias" );
93+ addSubGenerator (1 , "Hadron triggered" );
94+ GeneratorPythia8 ::Init ();
95+ pythiaMBgen .init ();
96+ return true;
97+ }
98+
99+ // search for the presence of at least one of the required hadrons in a selected rapidity window
100+ bool findHadrons (Pythia8 ::Event & event ) {
101+
102+ for (int ipa = 0 ; ipa < event .size (); ++ ipa ) {
103+
104+ auto daughterList = event [ipa ].daughterList ();
105+
106+ for (auto ida : daughterList ) {
107+ for (int pdg : mHadronsPDGs ) { // check that at least one of the pdg code is found in the event
108+ if (event [ida ].id () == pdg ) {
109+ if ((event [ida ].y () > mRapidityMin ) && (event [ida ].y () < mRapidityMax )) {
110+ cout << "============= Found jpsi y,pt " << event [ida ].y () << ", " << event [ida ].pT () << endl ;
111+ std ::vector < int > daughters = event [ida ].daughterList ();
112+ for (int d : daughters ) {
113+ cout << "###### daughter " << d << ": code " << event [d ].id () << ", pt " << event [d ].pT () << endl ;
114+ }
115+ return true;
116+ }
117+ }
118+ }
119+ }
120+ }
121+
122+ return false;
123+ };
124+
125+
126+ private :
127+ // Interface to override import particles
128+ Pythia8 ::Event mOutputEvent ;
129+
130+ // Control gap-triggering
131+ unsigned long long mGeneratedEvents ;
132+ int mInverseTriggerRatio ;
133+ Pythia8 ::Pythia pythiaMBgen ; // minimum bias event
134+ TString mConfigMBdecays ;
135+ std ::vector < int > mHadronsPDGs ;
136+ double mRapidityMin ;
137+ double mRapidityMax ;
138+ bool mVerbose ;
139+ };
140+
141+ }
142+
143+ }
144+
145+ // Predefined generators:
146+ FairGenerator *
147+ GeneratorInclusiveJpsi_EvtGenMidY (int triggerGap , double rapidityMin = -1.5 , double rapidityMax = 1.5 , bool verbose = false)
148+ {
149+ auto gen = new o2 ::eventgen ::GeneratorEvtGen < o2 ::eventgen ::GeneratorPythia8HadronTriggeredWithGap > ();
150+ gen -> setTriggerGap (triggerGap );
151+ gen -> setRapidityRange (rapidityMin , rapidityMax );
152+ gen -> addHadronPDGs (443 );
153+ gen -> setVerbose (verbose );
154+
155+ TString pathO2table = gSystem -> ExpandPathName ("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGDQ/pythia8/decayer/switchOffJpsi.cfg" );
156+ gen -> readFile (pathO2table .Data ());
157+ gen -> setConfigMBdecays (pathO2table );
158+ gen -> PrintDebug (true);
159+
160+ gen -> SetSizePdg (1 );
161+ gen -> AddPdg (443 , 0 );
162+
163+ gen -> SetForceDecay (kEvtDiElectron );
164+
165+ // set random seed
166+ gen -> readString ("Random:setSeed on" );
167+ uint random_seed ;
168+ unsigned long long int random_value = 0 ;
169+ ifstream urandom ("/dev/urandom" , ios ::in |ios ::binary );
170+ urandom .read (reinterpret_cast < char * > (& random_value ), sizeof (random_seed ));
171+ gen -> readString (Form ("Random:seed = %llu" , random_value % 900000001 ));
172+
173+ // print debug
174+ // gen->PrintDebug();
175+
176+ return gen ;
177+ }
0 commit comments