GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
COHHadronicSystemGenerator.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  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <cstdlib>
12 
13 #include <TVector3.h>
14 
30 
31 using namespace genie;
32 using namespace genie::constants;
33 
34 //___________________________________________________________________________
36  HadronicSystemGenerator("genie::COHHadronicSystemGenerator")
37 {
38 
39 }
40 //___________________________________________________________________________
42  HadronicSystemGenerator("genie::COHHadronicSystemGenerator", config)
43 {
44 
45 }
46 //___________________________________________________________________________
48 {
49 
50 }
51 //___________________________________________________________________________
53 {
54  // Access cross section algorithm for running thread
56  const EventGeneratorI * evg = rtinfo->RunningThread();
57  const XSecAlgorithmI *fXSecModel = evg->CrossSectionAlg();
58  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
60  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")) {
62  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")) {
64  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
66  }
67  else {
68  LOG("COHHadronicSystemGenerator",pFATAL) <<
69  "ProcessEventRecord >> Cannot calculate hadronic system for " <<
70  fXSecModel->Id().Name();
71  }
72 }
73 //___________________________________________________________________________
75 {
76  int pion_pdgc = 0;
77  if (xcls_tag.NPi0() == 1) pion_pdgc = kPdgPi0;
78  else if (xcls_tag.NPiPlus() == 1) pion_pdgc = kPdgPiP;
79  else if (xcls_tag.NPiMinus() == 1) pion_pdgc = kPdgPiM;
80  else {
81  LOG("COHHadronicVtx", pFATAL)
82  << "No final state pion information in XclsTag!";
83  exit(1);
84  }
85  return pion_pdgc;
86 }
87 //___________________________________________________________________________
89 {
90  // Treatment of the hadronic side is identical to Rein-Sehgal if we assume an infinite
91  // mass for the nucleus.
93 }
94 //___________________________________________________________________________
96 {
97  //
98  // This method generates the final state hadronic system (pion + nucleus) in
99  // COH interactions
100  //
101  RandomGen * rnd = RandomGen::Instance();
102 
103  Interaction * interaction = evrec->Summary();
104  const XclsTag & xcls_tag = interaction->ExclTag();
105  const InitialState & init_state = interaction -> InitState();
106 
107  //-- Access neutrino, initial nucleus and final state prim. lepton entries
108  GHepParticle * nu = evrec->Probe();
109  GHepParticle * Ni = evrec->TargetNucleus();
110  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
111  assert(nu);
112  assert(Ni);
113  assert(fsl);
114 
115  const TLorentzVector & vtx = *(nu->X4());
116  const TLorentzVector & p4nu = *(nu ->P4());
117  const TLorentzVector & p4fsl = *(fsl->P4());
118 
119  //-- Determine the pdg code of the final state pion & nucleus
120  int nucl_pdgc = Ni->Pdg(); // same as the initial nucleus
121  int pion_pdgc = getPionPDGCodeFromXclTag(xcls_tag);
122 
123  //-- basic kinematic inputs
124  double E = nu->E();
125  double Q2 = interaction->Kine().Q2(true);
126  double y = interaction->Kine().y(true);
127  double t = interaction->Kine().t(true);
128  double MA = init_state.Tgt().Mass();
129  // double MA2 = TMath::Power(MA, 2.); // Unused
130  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
131  double mpi2 = TMath::Power(mpi,2);
132 
133  SLOG("COHHadronicVtx", pINFO)
134  << "Ev = "<< E << ", Q^2 = " << Q2
135  << ", y = " << y << ", t = " << t;
136 
137  double Epi = y * E - t / (2 * MA);
138  double ppi2 = Epi * Epi - mpi2;
139  double ppi = ppi2 > 0.0 ? TMath::Sqrt(ppi2) : 0.0;
140 
141  double costheta = (t - Q2 - mpi2) / (2 * ( (y *E - Epi) * Epi -
142  ppi * sqrt(TMath::Power(y * E - Epi, 2.) + t) ) );
143 
144  if ((costheta > 1.0) || (costheta < -1.0)) {
145  SLOG("COHHadronicVtx", pERROR)
146  << "Unphysical pion angle!";
147  }
148 
149  double sintheta = TMath::Sqrt(1 - costheta * costheta);
150 
151  //-- first work in the c.m.s. frame
152  // double S = 2 * MA * nuh - Q2 + MA2;
153  // double S_2 = S >= 0 ? TMath::Sqrt(S) : 0.0; // TODO - Error here?
154  // double Pcm = MA * TMath::Sqrt( (nuh*nuh + Q2)/S );
155  // double Epi = (S + mpi2 - MA2)/(2 * S_2);
156  // double EAprime = (S - mpi2 + MA2)/(2 * S_2);
157  // double EA = (S + MA2 + Q2)/(2 * S_2);
158  // double PAprime2 = TMath::Power(EAprime,2.0) - MA2;
159  // double PAprime = TMath::Sqrt(PAprime2);
160  // double tA = TMath::Power((EAprime - EA),2.0) - TMath::Power(PAprime,2.0) -
161  // TMath::Power(Pcm, 2.0);
162  // double tB = 2 * Pcm * PAprime;
163  // double cosT = (t - tA)/tB;
164  // double sinT = TMath::Sqrt(1 - cosT*cosT);
165  // double PAz = PAprime * cosT;
166  // double PAperp = PAprime * sinT;
167  // double PPiz = -PAz;
168 
169  // Randomize transverse components
170  double phi = 2 * kPi * rnd->RndHadro().Rndm();
171  double ppix = ppi * sintheta * TMath::Cos(phi);
172  double ppiy = ppi * sintheta * TMath::Sin(phi);
173  double ppiz = ppi * costheta;
174 
175  // boost back to the lab frame
176  // double beta = TMath::Sqrt( nuh*nuh + Q2 )/(nuh + MA);
177  // double gamma = (nuh + MA)/TMath::Sqrt(S);
178  // double betagamma = beta * gamma;
179 
180  // double epi = gamma*Epi + betagamma*PPiz;
181  // double ppiz = betagamma*Epi + gamma*PPiz;
182 
183  // double ea = gamma*EAprime + betagamma*PAz;
184  // double paz = betagamma*EAprime + gamma*PAz;
185 
186  // Now rotate so our axes are aligned with the lab instead of q
187  TLorentzVector q = p4nu - p4fsl;
188  TVector3 ppi3(ppix, ppiy, ppiz);
189  ppi3.RotateUz(q.Vect().Unit());
190 
191  // Nucleus...
192  // TVector3 pa(PAx,PAy,paz);
193  // pa.RotateUz(q.Vect().Unit());
194 
195  // now figure out the f/s nucleus 4-p
196 
197  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
198  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
199  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
200  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
201 
202  //-- Save the particles at the GHEP record
203 
204  int mom = evrec->TargetNucleusPosition();
205 
206  // Nucleus - need to balance overall 4-momentum
207  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
208  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
209 
210  // evrec->AddParticle(
211  // nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
212  // pa.Px(), pa.Py(), pa.Pz(), ea, 0, 0, 0, 0);
213 
214  evrec->AddParticle(
215  pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
216  ppi3.Px(), ppi3.Py(), ppi3.Pz(), Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
217 }
218 //___________________________________________________________________________
220 {
221  //
222  // This method generates the final state hadronic system (pion + nucleus) in
223  // COH interactions
224  //
225  RandomGen * rnd = RandomGen::Instance();
226 
227  Interaction * interaction = evrec->Summary();
228  const XclsTag & xcls_tag = interaction->ExclTag();
229 
230  //-- Access neutrino, initial nucleus and final state prim. lepton entries
231  GHepParticle * nu = evrec->Probe();
232  GHepParticle * Ni = evrec->TargetNucleus();
233  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
234  assert(nu);
235  assert(Ni);
236  assert(fsl);
237 
238  const TLorentzVector & vtx = *(nu->X4());
239  const TLorentzVector & p4nu = *(nu ->P4());
240  const TLorentzVector & p4fsl = *(fsl->P4());
241 
242  //-- Determine the pdg code of the final state pion & nucleus
243  int nucl_pdgc = Ni->Pdg(); // same as the initial nucleus
244  int pion_pdgc = getPionPDGCodeFromXclTag(xcls_tag);
245 
246  //-- basic kinematic inputs
247  double E = nu->E();
248  double M = kNucleonMass;
249  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
250  double mpi2 = TMath::Power(mpi,2);
251  double xo = interaction->Kine().x(true);
252  double yo = interaction->Kine().y(true);
253  double to = interaction->Kine().t(true);
254 
255  SLOG("COHHadronicVtx", pINFO)
256  << "Ev = "<< E << ", xo = " << xo
257  << ", yo = " << yo << ", to = " << to;
258 
259  //-- compute pion energy and |momentum|
260  double Epi = yo * E;
261  double Epi2 = TMath::Power(Epi,2);
262  double ppi2 = Epi2-mpi2;
263  double ppi = TMath::Sqrt(TMath::Max(0.,ppi2));
264 
265  SLOG("COHHadronicVtx", pINFO)
266  << "f/s pion E = " << Epi << ", |p| = " << ppi;
267  assert(Epi>mpi);
268 
269  //-- 4-momentum transfer q=p(neutrino) - p(f/s lepton)
270  // Note: m^2 = q^2 < 0
271  // Also, since the nucleus is heavy all energy loss is comminicated to
272  // the outgoing pion Rein & Sehgal, Nucl.Phys.B223.29-44(1983), p.35
273 
274  TLorentzVector q = p4nu - p4fsl;
275 
276  SLOG("COHHadronicVtx", pINFO)
277  << "\n 4-p transfer q @ LAB: " << utils::print::P4AsString(&q);
278 
279  //-- find angle theta between q and ppi (xi=costheta)
280  // note: t=|(ppi-q)^2|, Rein & Sehgal, Nucl.Phys.B223.29-44(1983), p.36
281 
282  double xi = 1. + M*xo/Epi - 0.5*mpi2/Epi2 - 0.5*to/Epi2;
283  xi /= TMath::Sqrt((1.+2.*M*xo/Epi)*(1.-mpi2/Epi2));
284 
285  double costheta = xi;
286  double sintheta = TMath::Sqrt(TMath::Max(0.,1.-xi*xi));
287 
288  SLOG("COHHadronicVtx", pINFO) << "cos(pion, q) = " << costheta;
289 
290  // compute transverse and longitudinal ppi components along q
291  double ppiL = ppi*costheta;
292  double ppiT = ppi*sintheta;
293 
294  double phi = 2*kPi* rnd->RndHadro().Rndm();
295 
296  TVector3 ppi3(0,ppiT,ppiL);
297 
298  ppi3.RotateZ(phi); // randomize transverse components
299  ppi3.RotateUz(q.Vect().Unit()); // align longit. component with q in LAB
300 
301  SLOG("COHHadronicVtx", pINFO)
302  << "Pion 3-p @ LAB: " << utils::print::Vec3AsString(&ppi3);
303 
304  // now figure out the f/s nucleus 4-p
305 
306  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
307  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
308  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
309  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
310 
311  //-- Save the particles at the GHEP record
312 
313  int mom = evrec->TargetNucleusPosition();
314 
315  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
316  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
317 
318  evrec->AddParticle(pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
319  ppi3.Px(), ppi3.Py(),ppi3.Pz(),Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
320 }
321 //___________________________________________________________________________
323 {
324  Interaction * interaction = evrec->Summary();
325  const Kinematics & kinematics = interaction -> Kine();
326  GHepParticle * nu = evrec->Probe();
327  GHepParticle * Ni = evrec->TargetNucleus();
328  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
329 
330  // Pion
331  const TLorentzVector ppi = kinematics.HadSystP4();
332  const TVector3 ppi3 = ppi.Vect();
333  const double Epi = ppi.E();
334  int pion_pdgc=0;
335  if ( interaction->ProcInfo().IsWeakCC() ) {
336  if( nu->Pdg() > 0 ){ // neutrino
337  pion_pdgc = kPdgPiP;
338  }
339  else{ // anti-neutrino
340  pion_pdgc = kPdgPiM;
341  }
342  }
343  else if ( interaction->ProcInfo().IsWeakNC() ) {
344  pion_pdgc = kPdgPi0;
345  }
346  else{
347  LOG("COHHadronicSystemGeneratorAR", pFATAL)
348  << "Could not determine pion involved in interaction";
349  exit(1);
350  }
351 
352  //
353  // Nucleus
354  int nucl_pdgc = Ni->Pdg(); // pdg of final nucleus same as the initial nucleus
355  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
356  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
357  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
358  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
359  //
360  // Both
361  const TLorentzVector & vtx = *(nu->X4());
362  int mom = evrec->TargetNucleusPosition();
363 
364  //
365  // Fill the records
366  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
367  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
368 
369  evrec->AddParticle(pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
370  ppi3.Px(), ppi3.Py(),ppi3.Pz(),Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
371 }
Cross Section Calculation Interface.
bool IsWeakCC(void) const
void CalculateHadronicSystem_AlvarezRuso(GHepRecord *event_rec) const
#define pERROR
Definition: Messenger.h:59
double E(void) const
Get energy.
Definition: GHepParticle.h:91
static const double kNucleonMass
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int getPionPDGCodeFromXclTag(const XclsTag &xcls_tag) const
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
#define pFATAL
Definition: Messenger.h:56
Defines the EventGeneratorI interface.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
int NPiMinus(void) const
Definition: XclsTag.h:60
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:27
void CalculateHadronicSystem_ReinSehgal(GHepRecord *event_rec) const
double x(bool selected=false) const
Definition: Kinematics.cxx:99
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:66
double Mass(void) const
Definition: Target.cxx:224
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:284
double y(bool selected=false) const
Definition: Kinematics.cxx:112
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
bool IsWeakNC(void) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string Name(void) const
Definition: AlgId.h:44
void CalculateHadronicSystem_BergerSehgalFM(GHepRecord *event_rec) const
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:333
const Kinematics & Kine(void) const
Definition: Interaction.h:71
int NPiPlus(void) const
Definition: XclsTag.h:59
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:293
#define pINFO
Definition: Messenger.h:62
int NPi0(void) const
Definition: XclsTag.h:58
void CalculateHadronicSystem_BergerSehgal(GHepRecord *event_rec) const
static RunningThreadInfo * Instance(void)
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
const int kPdgPiM
Definition: PDGCodes.h:159
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double t(bool selected=false) const
Definition: Kinematics.cxx:170
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
void ProcessEventRecord(GHepRecord *event_rec) const
const EventGeneratorI * RunningThread(void)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:45
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
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:80
Keep info on the event generation thread currently on charge. This is used so that event generation m...
virtual int TargetNucleusPosition(void) const
Definition: GHepRecord.cxx:370
Initial State information.
Definition: InitialState.h:48
double Py(void) const
Get Py.
Definition: GHepParticle.h:89