@@ -41,6 +41,8 @@ public:
4141
4242 void addHadronPDGs (int pdg ) { mHadronsPDGs .push_back (pdg ); mRejFactorPrompt .push_back (1.0 ); mRejFactorNonPrompt .push_back (1.0 );}
4343
44+ void decrementGeneratedEvents () {mGeneratedEvents -- ;}
45+
4446 void setRejFactorPrompt (int pdg , float rejFactor ) {
4547 for (size_t i = 0 ; i < mHadronsPDGs .size (); i ++ ) {
4648 if (pdg == mHadronsPDGs [i ]) {
@@ -68,6 +70,10 @@ public:
6870 void setConfigMBdecays (TString val ){mConfigMBdecays = val ;}
6971
7072 void setVerbose (bool val ) { mVerbose = val ; };
73+
74+ int getNGeneratedEvents () { return mGeneratedEvents ;}
75+
76+ int getTriggerGap () { return mInverseTriggerRatio ;}
7177
7278protected :
7379
@@ -123,6 +129,13 @@ bool isOpenBhadron(int pdg) {
123129 return ((abs (pdg ) >= 500 && abs (pdg ) < 599 ) || (abs (pdg ) >= 5000 && abs (pdg ) < 5999 )) && pdg != 553 ;
124130}
125131
132+ bool isCharmonium (int pdg ) {
133+ int firstQuark = ( abs (pdg ) % 100 - abs (pdg ) % 10 ) / 10 ;
134+ int secondQuark = ( abs (pdg ) % 1000 - abs (pdg ) % 100 ) / 100 ;
135+ int thirdQuark = ( abs (pdg ) % 10000 - abs (pdg ) % 1000 ) / 1000 ;
136+ return ( firstQuark == 4 && secondQuark == 4 && thirdQuark == 0 );
137+ }
138+
126139// search for the presence of at least one of the required hadrons in a selected rapidity window
127140bool findHadrons (Pythia8 ::Event & event ) {
128141 int ihad = 0 ;
@@ -133,21 +146,42 @@ bool findHadrons(Pythia8::Event& event) {
133146 for (auto ida : daughterList ) {
134147 ihad = 0 ;
135148 for (int pdg : mHadronsPDGs ) { // check that at least one of the pdg code is found in the event
136- if (event [ida ].id () == pdg ) {
149+ if (abs ( event [ida ].id () ) == pdg ) {
137150 if ((event [ida ].y () > mRapidityMin ) && (event [ida ].y () < mRapidityMax )) {
138- cout << "============= Found jpsi y,pt,pdg " << event [ida ].y () << ", " << event [ida ].pT () << ", " << event [ida ].pdg () << endl ;
151+ cout << "============= Found jpsi y,pt,pdg " << event [ida ].y () << ", " << event [ida ].pT () << ", " << event [ida ].id () << endl ;
139152 std ::vector < int > daughters = event [ida ].daughterList ();
140153 for (int d : daughters ) {
141154 cout << "###### daughter " << d << ": code " << event [d ].id () << ", pt " << event [d ].pT () << endl ;
142155 }
156+ if (event [ida ].daughter1 () == event [ida ].daughter2 () && event [ida ].daughter1 () > 0 ) {
157+ continue ; // particle has a carbon-copy as daughter, its daughter will already be considered for triggering
158+ }
143159
144160 // check whether particle is prompt or non-prompt, since rejection factor can depend on it
145161 bool isNonPrompt = false;
146162 if (isOpenBhadron (pdg )) {
147163 isNonPrompt = true;
148- cout << "particle is non-prompt" << endl
164+ LOGF ( info , "particle is non-prompt" );
149165 } else {
150- // we check the mother
166+ // check history
167+ int currentIdx = ida ;
168+ int currentPdg = pdg ;
169+ cout << "particle history: " ;
170+ while (isCharmonium (currentPdg ) && currentIdx >= 0 ) {
171+ currentIdx = event [currentIdx ].mother1 ();
172+ if (currentIdx < 0 ) {
173+ break ;
174+ }
175+ currentPdg = abs (event [currentIdx ].id ());
176+ cout << currentPdg << " " ;
177+ if (isOpenBhadron (currentPdg )) {
178+ isNonPrompt = true;
179+ LOGF (info , "particle is non-prompt" );
180+ break ;
181+ }
182+ }
183+ cout << endl ;
184+ /*// we check the mother
151185 int mother = event[ida].mother1();
152186 if (mother >= 0 && isOpenBhadron(event[mother].id())) {
153187 isNonPrompt = true;
@@ -164,12 +198,14 @@ bool findHadrons(Pythia8::Event& event) {
164198 isNonPrompt = false;
165199 cout << "particle is prompt, mother pdg: " << event[mother].id() << ", grand-mother pdg: "<< event[grandmother].id() << endl;
166200 }
167- }
201+ }*/
168202 }
169203
170204 // rejection factor given in the ini file
171205 float randomNumber = gRandom -> Rndm ();
172- if ((isPrompt && (randomNumber <= mRejFactorPrompt [ihad ])) || (isNonPrompt && (randomNumber <= mRejFactorNonPrompt [ihad ]))) {
206+ cout << randomNumber << " rej factor: " << (isNonPrompt ? mRejFactorNonPrompt [ihad ] : mRejFactorPrompt [ihad ]) << endl ;
207+ if ((!isNonPrompt && (randomNumber <= mRejFactorPrompt [ihad ])) || (isNonPrompt && (randomNumber <= mRejFactorNonPrompt [ihad ]))) {
208+ cout << "event triggered " << endl ;
173209 return true;
174210 }
175211 }
@@ -204,6 +240,25 @@ private:
204240
205241}
206242
243+
244+ o2 ::eventgen ::Trigger triggerPDGRap (double rapMin , double rapMax , int pdg , GeneratorPythia8HadronTriggeredWithGap * gen ) {
245+ auto trigger = [rapMin , rapMax , pdg , gen ](const std ::vector < TParticle > & particles ) -> bool {
246+ if (gen -> getTriggerGap () != 1 && gen -> getNGeneratedEvents () % gen -> getTriggerGap () != 1 ) {
247+ // this is a MB event
248+ return true;
249+ }
250+ for (const auto& p : particles ) {
251+ if (p .Y () > rapMin && p .Y () < rapMax && std ::abs (p .GetPdgCode ()) == pdg ) {
252+ return true;
253+ }
254+ }
255+ return false;
256+ };
257+ return trigger ;
258+ }
259+
260+
261+
207262// Predefined generators:
208263FairGenerator *
209264 GeneratorInclusiveJpsi_EvtGenMidY (int triggerGap , double rapidityMin = -1.5 , double rapidityMax = 1.5 , bool verbose = false)
@@ -337,7 +392,9 @@ GeneratorInclusiveAllQuarkonia_EvtGenMidY(int triggerGap, double rapidityMin = -
337392
338393 auto gen = new o2 ::eventgen ::GeneratorEvtGen < o2 ::eventgen ::GeneratorPythia8HadronTriggeredWithGap > (
339394 gen -> setTriggerGap (triggerGap );
340- gen -> setRapidityRange (rapidityMin , rapidityMax );
395+ // this is a trigger before EvtGen decays, after which the rapidities of the particles are modified.
396+ // The rapidity cut is then only applied in the triggerEvent after the decays
397+ gen -> setRapidityRange (rapidityMin - 1. , rapidityMax + 1. );
341398 // specify particles to be triggered
342399 for (int i = 0 ; i < 16 ; i ++ ) {
343400 gen -> addHadronPDGs (particleList [i ]);
@@ -350,42 +407,45 @@ GeneratorInclusiveAllQuarkonia_EvtGenMidY(int triggerGap, double rapidityMin = -
350407 // or the keyword can be used for only one particle (e.g. "pdg1:rejFactor1:prompt,pdg1:rejFactor2:non-prompt")
351408 TObjArray * objArray = rejFactors .Tokenize ("," );
352409 for (int i = 0 ; i < objArray -> GetEntries (); i ++ ) {
353- TString rejStr = objArray -> At (i );
410+ TString rejStr = TString ( objArray -> At (i ) -> GetName () );
354411 TObjArray * objArrayCurrent = rejStr .Tokenize (":" );
355412 if (objArrayCurrent -> GetEntries () != 2 && objArrayCurrent -> GetEntries () != 3 ) {
356413 LOGF (fatal , "Problem when configuring string for particle rejection factors: %s, incorrect length" , rejStr .Data ());
357414 }
358- if (!objArrayCurrent [1 ].IsFloat ()) {
415+ TString str0 = TString (objArrayCurrent -> At (0 )-> GetName ());
416+ TString str1 = TString (objArrayCurrent -> At (1 )-> GetName ());
417+ if (!str1 .IsFloat ()) {
359418 LOGF (fatal , "Problem when configuring string for particle rejection factors: %s, is not float" , rejStr .Data ());
360419 }
361- if (objArrayCurrent [ 0 ] .CompareTo ("prompt" ) == 0 ) {
420+ if (str0 .CompareTo ("prompt" ) == 0 ) {
362421 // Common switch for all prompt particles
363- for (int ihad = 0 ; ihad < 16 ; i ++ ) {
364- gen -> setRejFactorPrompt (particleList [i ], objArrayCurrent [ 1 ] .Atof ());
422+ for (int ihad = 0 ; ihad < 16 ; ihad ++ ) {
423+ gen -> setRejFactorPrompt (particleList [ihad ], str1 .Atof ());
365424 }
366425 continue ;
367426 }
368- if (objArrayCurrent [ 0 ] .CompareTo ("non-prompt" ) == 0 ) {
427+ if (str0 .CompareTo ("non-prompt" ) == 0 ) {
369428 // Common switch for all non-prompt particles
370- for (int ihad = 0 ; ihad < 16 ; i ++ ) {
371- gen -> setRejFactorNonPrompt (particleList [i ], objArrayCurrent [ 1 ] .Atof ());
429+ for (int ihad = 0 ; ihad < 16 ; ihad ++ ) {
430+ gen -> setRejFactorNonPrompt (particleList [ihad ], str1 .Atof ());
372431 }
373432 continue ;
374433 }
375- if (objArrayCurrent [ 0 ] .IsDigit ()) {
434+ if (str0 .IsDigit ()) {
376435 // Setting the rejection factor for a specific particle
377436 if (objArrayCurrent -> GetEntries () == 2 ) {
378- gen -> setRejFactorPrompt (objArrayCurrent [ 0 ] .Atoi (), objArrayCurrent [ 1 ] .Atof ());
379- gen -> setRejFactorNonPrompt (objArrayCurrent [ 0 ] .Atoi (), objArrayCurrent [ 1 ] .Atof ());
437+ gen -> setRejFactorPrompt (str0 .Atoi (), str1 .Atof ());
438+ gen -> setRejFactorNonPrompt (str0 .Atoi (), str1 .Atof ());
380439 continue ;
381440 }
382441 else {
383- if (objArrayCurrent [2 ].CompareTo ("prompt" ) == 0 ) {
384- gen -> setRejFactorPrompt (objArrayCurrent [0 ].Atoi (), objArrayCurrent [1 ].Atof ());
442+ TString str2 = TString (objArrayCurrent -> At (2 )-> GetName ());
443+ if (str2 .CompareTo ("prompt" ) == 0 ) {
444+ gen -> setRejFactorPrompt (str0 .Atoi (), str1 .Atof ());
385445 continue ;
386446 }
387- if (objArrayCurrent [ 2 ] .CompareTo ("non-prompt" ) == 0 ) {
388- gen -> setRejFactorNonPrompt (objArrayCurrent [ 0 ] .Atoi (), objArrayCurrent [ 1 ] .Atof ());
447+ if (str2 .CompareTo ("non-prompt" ) == 0 ) {
448+ gen -> setRejFactorNonPrompt (str0 .Atoi (), str1 .Atof ());
389449 continue ;
390450 }
391451 }
@@ -417,6 +477,19 @@ GeneratorInclusiveAllQuarkonia_EvtGenMidY(int triggerGap, double rapidityMin = -
417477
418478 // print debug
419479 // gen->PrintDebug();
480+
481+ // add trigger on the correct rapidity range after EvtGen decays
482+ gen -> setTriggerMode (Generator ::kTriggerOR );
483+ for (int i = 0 ; i < 16 ; i ++ ) {
484+ gen -> addTrigger (triggerPDGRap (rapidityMin , rapidityMax , particleList [i ], gen ));
485+ }
486+
487+ // what to do if the trigger was rejected
488+ gen -> setTriggerFalseHook (
489+ [gen ](std ::vector < TParticle > const & p , int eventCount ) {
490+ gen -> decrementGeneratedEvents ();
491+ }
492+ );
420493
421494 return gen ;
422495}
0 commit comments