1+ #if !defined(__CLING__ ) || defined(__ROOTCLING__ )
2+ #include "FairGenerator.h"
3+ #include "FairPrimaryGenerator.h"
4+ #include "Generators/GeneratorPythia8.h"
5+ #include "Pythia8/Pythia.h"
6+ #include "TDatabasePDG.h"
7+ #include "TMath.h"
8+ #include "TParticlePDG.h"
9+ #include "TRandom3.h"
10+ #include "TSystem.h"
11+ #include "TVector2.h"
12+ #include "fairlogger/Logger.h"
13+ #include <cmath>
14+ #include <fstream>
15+ #include <string>
16+ #include "TLorentzVector.h"
17+ #include <vector>
18+ using namespace Pythia8 ;
19+ #endif
20+
21+ /// Event generator using Pythia ropes
22+ /// Triggers events containing at least two generated phi(1020) mesons.
23+
24+ class GeneratorPythia8DoublePhi : public o2 ::eventgen ::GeneratorPythia8
25+ {
26+ public :
27+ /// Constructor
28+ GeneratorPythia8DoublePhi (int gapSize = 0 , double minPt = 0.0 , double maxPt = 100.0 , double maxEta = 0.8 )
29+ : o2 ::eventgen ::GeneratorPythia8 (),
30+ mGapSize (gapSize ),
31+ mMinPt (minPt ),
32+ mMaxPt (maxPt ),
33+ mMaxEta (maxEta )
34+ {
35+ fmt ::printf (">> Pythia8 generator: two phi(1020) mesons, gap = %d, minPtPhi = %f, maxPtPhi = %f, |etaPhi| < %f\n" , gapSize , minPt , maxPt , maxEta );
36+ }
37+ /// Destructor
38+ ~GeneratorPythia8DoublePhi () = default ;
39+
40+ bool Init () override
41+ {
42+ addSubGenerator (0 , "Pythia8 events with two phi(1020) mesons" );
43+ return o2 ::eventgen ::GeneratorPythia8 ::Init ();
44+ }
45+
46+ protected :
47+ bool isPhiFromHFDecay (const Pythia8 ::Particle & p , const Pythia8 ::Event & event )
48+ {
49+
50+ // Walk up ancestry
51+ int motherId = p .mother1 ();
52+
53+ while (motherId > 0 )
54+ {
55+ // Get mother
56+ const auto & mother = event [motherId ];
57+ const int absMotherPdg = std ::abs (mother .id ());
58+
59+ // Check if particle is from HF decay
60+ if (((absMotherPdg / 100 ) % 10 == 4 ) ||
61+ ((absMotherPdg / 100 ) % 10 == 5 ) ||
62+ ((absMotherPdg / 1000 ) % 10 == 4 ) ||
63+ ((absMotherPdg / 1000 ) % 10 == 5 ))
64+ {
65+ return true;
66+ }
67+
68+ motherId = mother .mother1 ();
69+ }
70+ return false;
71+ }
72+
73+ bool generateEvent () override
74+ {
75+ // fmt::printf(">> Generating event %d\n", mGeneratedEvents);
76+
77+ bool genOk = false;
78+ int localCounter {0 };
79+ constexpr int kMaxTries {300000 };
80+
81+ // If mGapSize <= 0, filter ALL events to contain two phis.
82+ // Otherwise, generate mGapSize gap events before 1 triggered event.
83+ if (mGapSize > 0 && (mGeneratedEvents % (mGapSize + 1 ) < mGapSize ))
84+ {
85+ genOk = GeneratorPythia8 ::generateEvent ();
86+ // fmt::printf(">> Gap-event (no phi check)\n");
87+ }
88+ else
89+ {
90+ while (!genOk && localCounter < kMaxTries )
91+ {
92+ if (GeneratorPythia8 ::generateEvent ())
93+ {
94+ genOk = selectEvent (mPythia .event );
95+ }
96+ localCounter ++ ;
97+ }
98+ if (!genOk )
99+ {
100+ fmt ::printf ("Failed to generate triggered event after %d tries\n" , kMaxTries );
101+ return false;
102+ }
103+ fmt ::printf (">> Triggered event: event accepted after %d iterations (double phi(1020))\n" , localCounter );
104+ }
105+
106+ notifySubGenerator (0 );
107+ mGeneratedEvents ++ ;
108+ return true;
109+ }
110+
111+ bool selectEvent (Pythia8 ::Event & event )
112+ {
113+ std ::vector < TLorentzVector > phiCandidates ;
114+
115+ for (int i = 0 ; i < event .size (); i ++ )
116+ {
117+ const auto & p = event [i ];
118+
119+ if (std ::abs (p .id ()) != 333 )
120+ continue ;
121+
122+ if (p .pT () < mMinPt || p .pT () > mMaxPt )
123+ continue ;
124+
125+ if (std ::abs (p .eta ()) > mMaxEta )
126+ continue ;
127+
128+ if (isPhiFromHFDecay (p , event ))
129+ continue ;
130+
131+ // Avoid double-counting copy/history entries:
132+ // Ensure this is the physical produced phi (e.g. check if its daughter is a copy of itself)
133+ int d1 = p .daughter1 ();
134+ int d2 = p .daughter2 ();
135+ if (d1 > 0 && d1 == d2 && std ::abs (event [d1 ].id ()) == 333 )
136+ {
137+ // p decayed into another copy of phi, so skip this intermediate entry
138+ continue ;
139+ }
140+
141+ TLorentzVector phi ;
142+ phi .SetPtEtaPhiM (p .pT (), p .eta (), p .phi (), p .m ());
143+
144+ phiCandidates .push_back (phi );
145+ }
146+ if (phiCandidates .size () < 2 )
147+ return false;
148+
149+ // Check all possible phi-phi pairs
150+ for (size_t i = 0 ; i < phiCandidates .size (); i ++ )
151+ {
152+ for (size_t j = i + 1 ; j < phiCandidates .size (); j ++ )
153+ {
154+ TLorentzVector phiPhi = phiCandidates [i ] + phiCandidates [j ];
155+
156+ double mass = phiPhi .M ();
157+ double pt = phiPhi .Pt ();
158+
159+ if (mass > 2.4 && pt > 4.0 )
160+ return true;
161+ }
162+ }
163+
164+ return false;
165+ }
166+
167+ private :
168+ int mGapSize {0 };
169+ double mMinPt {0.0 };
170+ double mMaxPt {100.0 };
171+ double mMaxEta {0.8 };
172+ uint64_t mGeneratedEvents {0 };
173+ };
174+
175+ ///___________________________________________________________
176+ FairGenerator * generateDoublePhi (int gap = 0 , double minPt = 0.0 , double maxPt = 100.0 , double maxEta = 0.8 )
177+ {
178+ auto myGenerator = new GeneratorPythia8DoublePhi (gap , minPt , maxPt , maxEta );
179+
180+ myGenerator -> readString ("333:onMode = off" );
181+ myGenerator -> readString ("333:onIfMatch = 321 -321" );
182+
183+ auto seed = (gRandom -> TRandom ::GetSeed () % 900000000 );
184+ myGenerator -> readString ("Random:setSeed on" );
185+ myGenerator -> readString ("Random:seed " + std ::to_string (seed ));
186+
187+ return myGenerator ;
188+ }
0 commit comments