GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
PhotonRESGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Alfonso Garcia <aagarciasoto \at km3net.de>
7  IFIC & Harvard University
8 */
9 //____________________________________________________________________________
10 
11 #include <cstring>
12 
13 #include <RVersion.h>
14 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6)
15 #include <TMCParticle.h>
16 #else
17 #include <TMCParticle6.h>
18 #endif
19 #include <TClonesArray.h>
20 #include <TMath.h>
21 
36 
37 using namespace genie;
38 using namespace genie::constants;
39 using namespace genie::utils::math;
40 
41 //___________________________________________________________________________
43 EventRecordVisitorI("genie::PhotonRESGenerator")
44 {
45  this->Initialize();
46  born = new Born();
47 }
48 //___________________________________________________________________________
50 EventRecordVisitorI("genie::PhotonRESGenerator", config)
51 {
52  this->Initialize();
53 }
54 //___________________________________________________________________________
56 {
57 
58 }
59 //____________________________________________________________________________
61 {
62  fPythia = TPythia6::Instance();
63 
64  // sync GENIE/PYTHIA6 seed number
66 }
67 //___________________________________________________________________________
69 {
70 
71  Interaction * interaction = evrec->Summary();
72  const InitialState & init_state = interaction->InitState();
73 
74  //incoming v & struck particle & remnant nucleus
75  GHepParticle * nu = evrec->Probe();
76  GHepParticle * el = evrec->HitNucleon();
77 
78  GHepParticle * target = evrec -> TargetNucleus();
79  if(target) evrec->AddParticle(target->Pdg(), kIStFinalStateNuclearRemnant, 1,-1,-1,-1, *(target->P4()), *(target->X4()) );
80 
81  TVector3 unit_nu = nu->P4()->Vect().Unit();
82 
83  int probepdg = init_state.ProbePdg();
84 
85  long double Mtarget = init_state.Tgt().HitNucMass();
86  long double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton
87 
88  long double Enuin = init_state.ProbeE(kRfLab);
89  long double s = born->GetS(Mtarget,Enuin);
90 
91  long double n1 = interaction->Kine().GetKV(kKVn1);
92  long double n2 = interaction->Kine().GetKV(kKVn2);
93 
94  long double costhCM = n1;
95  long double sinthCM = sqrtl(1-costhCM*costhCM);
96 
97  long double xmin = fQ2PDFmin/2./Enuin/Mtarget;
98  long double x = expl( logl(xmin) + (logl(1.0)-logl(xmin))*n2 );
99  long double s_r = s*x;
100 
101  // Boost velocity CM -> LAB
102  long double EnuinCM = sqrtl(s_r)/2.;
103  long double beta = (powl(Enuin,2)-powl(EnuinCM,2))/(powl(Enuin,2)+powl(EnuinCM,2));
104 
105  // Final state primary lepton PDG code
106  int pdgl = interaction->FSPrimLeptonPdg();
107  assert(pdgl!=0);
108 
109  if ( pdg::IsElectron(TMath::Abs(pdgl)) || pdg::IsMuon(TMath::Abs(pdgl)) || pdg::IsTau(TMath::Abs(pdgl)) ) {
110 
111  long double ElpoutCM = (s_r+mlout*mlout)/sqrtl(s_r)/2.;
112  long double EnuoutCM = (s_r-mlout*mlout)/sqrtl(s_r)/2.;
113  LongLorentzVector p4_lpout( 0., EnuoutCM*sinthCM, EnuoutCM*costhCM, ElpoutCM );
114  LongLorentzVector p4_nuout( 0., -EnuoutCM*sinthCM, -EnuoutCM*costhCM, EnuoutCM );
115 
116  p4_lpout.BoostZ(beta);
117  p4_nuout.BoostZ(beta);
118 
119  TLorentzVector p4lp_o( (double)p4_lpout.Px(), (double)p4_lpout.Py(), (double)p4_lpout.Pz(), (double)p4_lpout.E() );
120  TLorentzVector p4nu_o( (double)p4_nuout.Px(), (double)p4_nuout.Py(), (double)p4_nuout.Pz(), (double)p4_nuout.E() );
121 
122  // Randomize transverse components
123  RandomGen * rnd = RandomGen::Instance();
124  double phi = 2* kPi * rnd->RndLep().Rndm();
125  p4lp_o.RotateZ(phi);
126  p4nu_o.RotateZ(phi);
127 
128  //rotate from LAB=[0,0,Ev,Ev]->[px,py,pz,E]
129  p4lp_o.RotateUz(unit_nu);
130  p4nu_o.RotateUz(unit_nu);
131 
132  int pdgvout = 0;
133  if ( pdg::IsElectron(pdgl) ) pdgvout = kPdgAntiNuE;
134  else if ( pdg::IsPositron(pdgl) ) pdgvout = kPdgNuE;
135  else if ( pdg::IsMuon(pdgl) ) pdgvout = kPdgAntiNuMu;
136  else if ( pdg::IsAntiMuon(pdgl) ) pdgvout = kPdgNuMu;
137  else if ( pdg::IsTau(pdgl) ) pdgvout = kPdgAntiNuTau;
138  else if ( pdg::IsAntiTau(pdgl) ) pdgvout = kPdgNuTau;
139 
140  int pdgboson = pdg::IsNeutrino(probepdg) ? kPdgWP : kPdgWM;
141 
142  // Create a GHepParticle and add it to the event record
143  evrec->AddParticle( pdgboson, kIStDecayedState, 0, -1, 5, 6, *nu->P4()+*el->P4(), *(nu->X4()) ); //W [mothers: nuebar_in,e_in][daugthers: nulbar_out,lp_out]
144  evrec->AddParticle( pdgl, kIStStableFinalState, 4, -1, -1, -1, p4lp_o, *(nu->X4()) );
145  evrec->AddParticle( pdgvout, kIStStableFinalState, 4, -1, -1, -1, p4nu_o, *(nu->X4()) );
146  evrec->Summary()->KinePtr()->SetFSLeptonP4(p4lp_o);
147 
148  }
149  else {
150 
151  char p6frame[10];
152  strcpy(p6frame, "CMS" );
153 
154  char p6nu[10], p6tgt[10];
155  if (pdg::IsAntiNeutrino(nu->Pdg())) { strcpy(p6nu, "nu_ebar"); strcpy(p6tgt, "e-"); }
156  else if (pdg::IsNeutrino (nu->Pdg())) { strcpy(p6nu, "nu_e"); strcpy(p6tgt, "e+"); }
157 
158  int def61 = fPythia->GetMSTP(61);
159  int def71 = fPythia->GetMSTP(71);
160  int def206 = fPythia->GetMDME(206,1);
161  int def207 = fPythia->GetMDME(207,1);
162  int def208 = fPythia->GetMDME(208,1);
163  fPythia->SetMSTP(61,0); // (Default=2) master switch for initial-state QCD and QED radiation.
164  fPythia->SetMSTP(71,0); // (Default=2) master switch for initial-state QCD and QED radiation.
165  fPythia->SetMDME(206,1,0); //swicht off W decay leptonic modes
166  fPythia->SetMDME(207,1,0);
167  fPythia->SetMDME(208,1,0);
168 
169  fPythia->Pyinit(p6frame, p6nu, p6tgt, sqrtl(s_r));
170  fPythia->Pyevnt();
171 
172  fPythia->SetMSTP(61,def61);
173  fPythia->SetMSTP(71,def71);
174  fPythia->SetMDME(206,1,def206);
175  fPythia->SetMDME(207,1,def207);
176  fPythia->SetMDME(208,1,def208);
177 
178  // get LUJETS record
179  fPythia->GetPrimaries();
180  TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All");
181  int np = pythia_particles->GetEntries();
182  assert(np>0);
183 
184  TMCParticle * particle = 0;
185  TIter piter(pythia_particles);
186  while( (particle = (TMCParticle *) piter.Next()) ) {
187 
188  int pdgc = particle->GetKF();
189  int ks = particle->GetKS();
190 
191  if ( ks==21 ) { continue; } //we dont want to save first particles from pythia (init states)
192 
193  LongLorentzVector p4longo(particle->GetPx(), particle->GetPy(), particle->GetPz(), particle->GetEnergy());
194  p4longo.BoostZ(beta);
195 
196  TLorentzVector p4o( (double)p4longo.Px(), (double)p4longo.Py(), (double)p4longo.Pz(), (double)p4longo.E() );
197  p4o.RotateUz(unit_nu);
198 
199  TParticlePDG * part = PDGLibrary::Instance()->Find(pdgc);
200  if ( (ks==1 || ks==4) && p4o.E() < part->Mass() ) {
201  LOG("PhotonRESGenerator", pWARN) << "Putting at rest one stable particle generated by PYTHIA because E < m";
202  LOG("PhotonRESGenerator", pWARN) << "PDG = " << pdgc << " // State = " << ks;
203  LOG("PhotonRESGenerator", pWARN) << "E = " << p4o.E() << " // |p| = " << TMath::Sqrt(p4o.P());
204  LOG("PhotonRESGenerator", pWARN) << "p = [ " << p4o.Px() << " , " << p4o.Py() << " , " << p4o.Pz() << " ]";
205  LOG("PhotonRESGenerator", pWARN) << "m = " << p4o.M() << " // mpdg = " << part->Mass();
206  p4o.SetXYZT(0,0,0,part->Mass());
207  }
208 
209  // copy final state particles to the event record
211 
212  // fix numbering scheme used for mother/daughter assignments
213  int firstmother = -1;
214  int lastmother = -1;
215  int firstchild = -1;
216  int lastchild = -1;
217 
218  if ( particle->GetParent() < 10 ) {
219  if ( TMath::Abs(pdgc)<7 ) { //outgoing quarks: mother will be the boson (saved in position 4)
220  firstmother = 4;
221  firstchild = particle->GetFirstChild() - 6;
222  lastchild = particle->GetLastChild() - 6;
223  }
224  else if ( TMath::Abs(pdgc)==24 ) { //produced W boson: mother will be the incoming neutrino
225  firstmother = 0;
226  firstchild = particle->GetFirstChild() - 6;
227  lastchild = particle->GetLastChild() - 6;
228  }
229  else if ( pdgc==22 ) { //radiative photons: mother will be the incoming electron
230  firstmother = 2;
231  }
232  }
233  else { //rest
234  firstmother = particle->GetParent() - 6; //shift to match boson position
235  firstchild = (particle->GetFirstChild()==0) ? particle->GetFirstChild() - 1 : particle->GetFirstChild() - 6;
236  lastchild = (particle->GetLastChild()==0) ? particle->GetLastChild() - 1 : particle->GetLastChild() - 6;
237  }
238 
239  double lightspeed = 299792458e3; //c in mm/s. Used for time in PYTHIA t[s]=t_pythia[mm]/c[mm/s]
240  double vx = nu->X4()->X() + particle->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm]
241  double vy = nu->X4()->Y() + particle->GetVy()*1e12;
242  double vz = nu->X4()->Z() + particle->GetVz()*1e12;
243  double vt = nu->X4()->T() + particle->GetTime()/lightspeed;
244  TLorentzVector pos( vx, vy, vz, vt );
245 
246  evrec->AddParticle(pdgc, ist, firstmother, lastmother, firstchild, lastchild, p4o, pos );
247 
248  }
249 
250  delete particle;
251  pythia_particles->Clear("C");
252 
253  }
254 
255 }
256 //___________________________________________________________________________
258 {
259  Algorithm::Configure(config);
260  this->LoadConfig();
261 }
262 //____________________________________________________________________________
263 void PhotonRESGenerator::Configure(string config)
264 {
265  Algorithm::Configure(config);
266  this->LoadConfig();
267 }
268 //____________________________________________________________________________
270 {
271 
272  // PYTHIA parameters only valid for HEDIS
273  double wmin; GetParam( "Xsec-Wmin", wmin ) ;
274  int warnings; GetParam( "Warnings", warnings ) ;
275  int errors; GetParam( "Errors", errors ) ;
276  int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ;
277  fPythia->SetPARP(2, wmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes)
278  fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed
279  fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed
280  fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term.
281 
282  fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385)
283  fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation
284  fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation
285 
286  GetParam( "Q2Grid-Min", fQ2PDFmin );
287 
288 }
289 //____________________________________________________________________________
const int kPdgNuE
Definition: PDGCodes.h:28
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:91
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:78
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Defines the EventRecordVisitorI interface. Concrete implementations of this interface use the &#39;Visito...
const int kPdgAntiNuE
Definition: PDGCodes.h:29
double HitNucMass(void) const
Definition: Target.cxx:233
const int kPdgWM
Definition: PDGCodes.h:192
const int kPdgNuMu
Definition: PDGCodes.h:30
static constexpr double s
Definition: Units.h:95
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int FSPrimLeptonPdg(void) const
final state primary lepton pdg
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
void SetFSLeptonP4(const TLorentzVector &p4)
Definition: Kinematics.cxx:297
TPythia6 * fPythia
PYTHIA6 wrapper class.
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
const Kinematics & Kine(void) const
Definition: Interaction.h:71
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
int ProbePdg(void) const
Definition: InitialState.h:64
bool IsAntiMuon(int pdgc)
Definition: PDGUtils.cxx:203
bool IsTau(int pdgc)
Definition: PDGUtils.cxx:208
void Configure(const Registry &config)
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
bool IsPositron(int pdgc)
Definition: PDGUtils.cxx:193
bool IsMuon(int pdgc)
Definition: PDGUtils.cxx:198
double GetS(double mlin, double Enuin)
Definition: Born.cxx:195
const int kPdgAntiNuTau
Definition: PDGCodes.h:33
const int kPdgAntiNuMu
Definition: PDGCodes.h:31
#define pWARN
Definition: Messenger.h:60
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const int kPdgNuTau
Definition: PDGCodes.h:32
Born level nu-electron cross section.
Definition: Born.h:26
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const InitialState & InitState(void) const
Definition: Interaction.h:69
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
double ProbeE(RefFrame_t rf) const
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
const int kPdgWP
Definition: PDGCodes.h:191
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
bool IsAntiTau(int pdgc)
Definition: PDGUtils.cxx:213
bool IsElectron(int pdgc)
Definition: PDGUtils.cxx:188
enum genie::EGHepStatus GHepStatus_t
void ProcessEventRecord(GHepRecord *evrec) const
Initial State information.
Definition: InitialState.h:48