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 <vector>
17+ using namespace Pythia8 ;
18+ #endif
19+
20+ /// Event of interest:
21+ /// an event that contains at least one Sigma- or Sigma+
22+ /// paired with at least one hadron of a specific PDG code,
23+ /// and the pair has k* < kStarMax
24+
25+ class GeneratorPythia8SigmaHadron : public o2 ::eventgen ::GeneratorPythia8
26+ {
27+ public :
28+ /// Constructor
29+ GeneratorPythia8SigmaHadron (int hadronPdg , int gapSize = 4 , double minPt = 0.2 ,
30+ double maxPt = 10 , double maxEta = 0.8 , double kStarMax = 1.0 )
31+ : o2 ::eventgen ::GeneratorPythia8 (),
32+ mHadronPdg (hadronPdg ),
33+ mGapSize (gapSize ),
34+ mMinPt (minPt ),
35+ mMaxPt (maxPt ),
36+ mMaxEta (maxEta ),
37+ mKStarMax (kStarMax )
38+ {
39+ fmt ::printf (
40+ ">> Pythia8 generator: Sigma± + hadron(PDG=%d), gap = %d, minPt = %f, maxPt = %f, |eta| < %f, k* < %f\n" ,
41+ hadronPdg , gapSize , minPt , maxPt , maxEta , kStarMax );
42+ }
43+
44+ ~GeneratorPythia8SigmaHadron () = default ;
45+
46+ bool Init () override
47+ {
48+ addSubGenerator (0 , "Pythia8 events with Sigma± and a specific hadron" );
49+ return o2 ::eventgen ::GeneratorPythia8 ::Init ();
50+ }
51+
52+ protected :
53+ /// Check whether particle descends from a Sigma+ or Sigma-
54+ bool isFromSigmaDecay (const Pythia8 ::Particle & p , const Pythia8 ::Event & event )
55+ {
56+ int motherId = p .mother1 ();
57+
58+ while (motherId > 0 ) {
59+ const auto& mother = event [motherId ];
60+ const int absMotherPdg = std ::abs (mother .id ());
61+
62+ if (absMotherPdg == 3112 || absMotherPdg == 3222 ) {
63+ return true;
64+ }
65+
66+ motherId = mother .mother1 ();
67+ }
68+
69+ return false;
70+ }
71+
72+ /// k* of a pair from invariant masses and 4-momenta
73+ /// k* = momentum of either particle in the pair rest frame
74+ double computeKStar (const Pythia8 ::Particle & p1 , const Pythia8 ::Particle & p2 ) const
75+ {
76+ const double e = p1 .e () + p2 .e ();
77+ const double px = p1 .px () + p2 .px ();
78+ const double py = p1 .py () + p2 .py ();
79+ const double pz = p1 .pz () + p2 .pz ();
80+ const double s = e * e - px * px - py * py - pz * pz ;
81+ if (s <= 0. ) {
82+ return 1.e9 ;
83+ }
84+ const double m1 = p1 .m ();
85+ const double m2 = p2 .m ();
86+ const double term1 = s - std ::pow (m1 + m2 , 2 );
87+ const double term2 = s - std ::pow (m1 - m2 , 2 );
88+ const double lambda = term1 * term2 ;
89+
90+ if (lambda <= 0. ) {
91+ return 0. ;
92+ }
93+ return 0.5 * std ::sqrt (lambda / s );
94+ }
95+
96+ bool generateEvent () override
97+ {
98+ fmt ::printf (">> Generating event %llu\n" , mGeneratedEvents );
99+
100+ bool genOk = false;
101+ int localCounter {0 };
102+ constexpr int kMaxTries {100000 };
103+
104+ // Accept mGapSize events unconditionally, then one triggered event
105+ if (mGeneratedEvents % (mGapSize + 1 ) < mGapSize ) {
106+ genOk = GeneratorPythia8 ::generateEvent ();
107+ fmt ::printf (">> Gap-event (no trigger check)\n" );
108+ } else {
109+ while (!genOk && localCounter < kMaxTries ) {
110+ if (GeneratorPythia8 ::generateEvent ()) {
111+ genOk = selectEvent (mPythia .event );
112+ }
113+ localCounter ++ ;
114+ }
115+
116+ if (!genOk ) {
117+ fmt ::printf ("Failed to generate triggered event after %d tries\n" , kMaxTries );
118+ return false;
119+ }
120+
121+ fmt ::printf (">> Triggered event: accepted after %d iterations (Sigma± + hadron PDG=%d, k* < %f)\n" ,
122+ localCounter , mHadronPdg , mKStarMax );
123+ }
124+
125+ notifySubGenerator (0 );
126+ mGeneratedEvents ++ ;
127+ return true;
128+ }
129+
130+ bool selectEvent (Pythia8 ::Event & event )
131+ {
132+ std ::vector < int > sigmaIndices ;
133+ std ::vector < int > hadronIndices ;
134+
135+ for (int i = 0 ; i < event .size (); i ++ ) {
136+ const auto& p = event [i ];
137+
138+ if (std ::abs (p .eta ()) > mMaxEta || p .pT () < mMinPt || p .pT () > mMaxPt ) {
139+ continue ;
140+ }
141+
142+ const int pdg = p .id ();
143+ const int absPdg = std ::abs (pdg );
144+
145+ // Sigma- or Sigma+
146+ if (absPdg == 3112 || absPdg == 3222 ) {
147+ sigmaIndices .push_back (i );
148+ }
149+
150+ if (std ::abs (pdg ) == mHadronPdg && !isFromSigmaDecay (p , event )) {
151+ hadronIndices .push_back (i );
152+ }
153+ }
154+ if (sigmaIndices .empty () || hadronIndices .empty ()) {
155+ return false;
156+ }
157+
158+ for (const auto iSigma : sigmaIndices ) {
159+ for (const auto iHadron : hadronIndices ) {
160+
161+ if (iSigma == iHadron ) {
162+ continue ;
163+ }
164+
165+ const auto& sigma = event [iSigma ];
166+ const auto & hadron = event [iHadron ];
167+
168+ const double kStar = computeKStar (sigma , hadron );
169+ if (kStar < mKStarMax ) {
170+ return true;
171+ }
172+ }
173+ }
174+
175+ return false;
176+ }
177+
178+ private :
179+ int mHadronPdg {211 };
180+ int mGapSize {4 };
181+ double mMinPt {0.2 };
182+ double mMaxPt {10.0 };
183+ double mMaxEta {0.8 };
184+ double mKStarMax {1.0 };
185+ uint64_t mGeneratedEvents {0 };
186+ };
187+
188+ ///___________________________________________________________
189+ FairGenerator * generateSigmaHadron (int hadronPdg , int gap = 4 , double minPt = 0.2 ,
190+ double maxPt = 10 , double maxEta = 0.8 , double kStarMax = 1.0 )
191+ {
192+ auto myGenerator = new GeneratorPythia8SigmaHadron (hadronPdg , gap , minPt , maxPt , maxEta , kStarMax );
193+ auto seed = (gRandom -> TRandom ::GetSeed () % 900000000 );
194+ myGenerator -> readString ("Random:setSeed on" );
195+ myGenerator -> readString ("Random:seed " + std ::to_string (seed ));
196+ return myGenerator ;
197+ }
0 commit comments