COHHadronicSystemGenerator.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5  or see $GENIE/LICENSE
6 
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 09, 2009 - CA
14  Moved into the new Coherent package from its previous location (EVGModules
15  package)
16  @ Mar 03, 2009 - CA
17  Renamed COHPiHadronicSystemGenerator -> COHHadronicSystemGenerator in
18  anticipation of reusing the code for simulating coherent production of
19  vector mesons.
20  @ Apr 02, 2009 - CA,HG,PK
21  Bug fix: Reverse the order of the pion momentum rotations: Randomize the
22  transverse component direction in the x'y' plane before aligning z' with
23  the direction of the momentum transfer q in the LAB.
24 */
25 //____________________________________________________________________________
26 
27 #include <cstdlib>
28 
29 #include <TVector3.h>
30 
46 
47 using namespace genie;
48 using namespace genie::constants;
49 
50 //___________________________________________________________________________
52  HadronicSystemGenerator("genie::COHHadronicSystemGenerator")
53 {
54 
55 }
56 //___________________________________________________________________________
58  HadronicSystemGenerator("genie::COHHadronicSystemGenerator", config)
59 {
60 
61 }
62 //___________________________________________________________________________
64 {
65 
66 }
67 //___________________________________________________________________________
69 {
70  // Access cross section algorithm for running thread
72  const EventGeneratorI * evg = rtinfo->RunningThread();
73  const XSecAlgorithmI *fXSecModel = evg->CrossSectionAlg();
74  if (fXSecModel->Id().Name() == "genie::ReinSehgalCOHPiPXSec") {
76  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiPXSec2015")) {
78  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalFMCOHPiPXSec2015")) {
80  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHPiPXSec")) {
82  }
83  else {
84  LOG("COHHadronicSystemGenerator",pFATAL) <<
85  "ProcessEventRecord >> Cannot calculate hadronic system for " <<
86  fXSecModel->Id().Name();
87  }
88 }
89 //___________________________________________________________________________
91 {
92  int pion_pdgc = 0;
93  if (xcls_tag.NPi0() == 1) pion_pdgc = kPdgPi0;
94  else if (xcls_tag.NPiPlus() == 1) pion_pdgc = kPdgPiP;
95  else if (xcls_tag.NPiMinus() == 1) pion_pdgc = kPdgPiM;
96  else {
97  LOG("COHHadronicVtx", pFATAL)
98  << "No final state pion information in XclsTag!";
99  exit(1);
100  }
101  return pion_pdgc;
102 }
103 //___________________________________________________________________________
105 {
106  // Treatment of the hadronic side is identical to Rein-Sehgal if we assume an infinite
107  // mass for the nucleus.
109 }
110 //___________________________________________________________________________
112 {
113  //
114  // This method generates the final state hadronic system (pion + nucleus) in
115  // COH interactions
116  //
118 
119  Interaction * interaction = evrec->Summary();
120  const XclsTag & xcls_tag = interaction->ExclTag();
121  const InitialState & init_state = interaction -> InitState();
122 
123  //-- Access neutrino, initial nucleus and final state prim. lepton entries
124  GHepParticle * nu = evrec->Probe();
125  GHepParticle * Ni = evrec->TargetNucleus();
126  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
127  assert(nu);
128  assert(Ni);
129  assert(fsl);
130 
131  const TLorentzVector & vtx = *(nu->X4());
132  const TLorentzVector & p4nu = *(nu ->P4());
133  const TLorentzVector & p4fsl = *(fsl->P4());
134 
135  //-- Determine the pdg code of the final state pion & nucleus
136  int nucl_pdgc = Ni->Pdg(); // same as the initial nucleus
137  int pion_pdgc = getPionPDGCodeFromXclTag(xcls_tag);
138 
139  //-- basic kinematic inputs
140  double E = nu->E();
141  double Q2 = interaction->Kine().Q2(true);
142  double y = interaction->Kine().y(true);
143  double t = interaction->Kine().t(true);
144  double MA = init_state.Tgt().Mass();
145  // double MA2 = TMath::Power(MA, 2.); // Unused
146  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
147  double mpi2 = TMath::Power(mpi,2);
148 
149  SLOG("COHHadronicVtx", pINFO)
150  << "Ev = "<< E << ", Q^2 = " << Q2
151  << ", y = " << y << ", t = " << t;
152 
153  double Epi = y * E - t / (2 * MA);
154  double ppi2 = Epi * Epi - mpi2;
155  double ppi = ppi2 > 0.0 ? TMath::Sqrt(ppi2) : 0.0;
156 
157  double costheta = (t - Q2 - mpi2) / (2 * ( (y *E - Epi) * Epi -
158  ppi * sqrt(TMath::Power(y * E - Epi, 2.) + t) ) );
159 
160  if ((costheta > 1.0) || (costheta < -1.0)) {
161  SLOG("COHHadronicVtx", pERROR)
162  << "Unphysical pion angle!";
163  }
164 
165  double sintheta = TMath::Sqrt(1 - costheta * costheta);
166 
167  //-- first work in the c.m.s. frame
168  // double S = 2 * MA * nuh - Q2 + MA2;
169  // double S_2 = S >= 0 ? TMath::Sqrt(S) : 0.0; // TODO - Error here?
170  // double Pcm = MA * TMath::Sqrt( (nuh*nuh + Q2)/S );
171  // double Epi = (S + mpi2 - MA2)/(2 * S_2);
172  // double EAprime = (S - mpi2 + MA2)/(2 * S_2);
173  // double EA = (S + MA2 + Q2)/(2 * S_2);
174  // double PAprime2 = TMath::Power(EAprime,2.0) - MA2;
175  // double PAprime = TMath::Sqrt(PAprime2);
176  // double tA = TMath::Power((EAprime - EA),2.0) - TMath::Power(PAprime,2.0) -
177  // TMath::Power(Pcm, 2.0);
178  // double tB = 2 * Pcm * PAprime;
179  // double cosT = (t - tA)/tB;
180  // double sinT = TMath::Sqrt(1 - cosT*cosT);
181  // double PAz = PAprime * cosT;
182  // double PAperp = PAprime * sinT;
183  // double PPiz = -PAz;
184 
185  // Randomize transverse components
186  double phi = 2 * kPi * rnd->RndHadro().Rndm();
187  double ppix = ppi * sintheta * TMath::Cos(phi);
188  double ppiy = ppi * sintheta * TMath::Sin(phi);
189  double ppiz = ppi * costheta;
190 
191  // boost back to the lab frame
192  // double beta = TMath::Sqrt( nuh*nuh + Q2 )/(nuh + MA);
193  // double gamma = (nuh + MA)/TMath::Sqrt(S);
194  // double betagamma = beta * gamma;
195 
196  // double epi = gamma*Epi + betagamma*PPiz;
197  // double ppiz = betagamma*Epi + gamma*PPiz;
198 
199  // double ea = gamma*EAprime + betagamma*PAz;
200  // double paz = betagamma*EAprime + gamma*PAz;
201 
202  // Now rotate so our axes are aligned with the lab instead of q
203  TLorentzVector q = p4nu - p4fsl;
204  TVector3 ppi3(ppix, ppiy, ppiz);
205  ppi3.RotateUz(q.Vect().Unit());
206 
207  // Nucleus...
208  // TVector3 pa(PAx,PAy,paz);
209  // pa.RotateUz(q.Vect().Unit());
210 
211  // now figure out the f/s nucleus 4-p
212 
213  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
214  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
215  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
216  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
217 
218  //-- Save the particles at the GHEP record
219 
220  int mom = evrec->TargetNucleusPosition();
221 
222  // Nucleus - need to balance overall 4-momentum
223  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
224  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
225 
226  // evrec->AddParticle(
227  // nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
228  // pa.Px(), pa.Py(), pa.Pz(), ea, 0, 0, 0, 0);
229 
230  evrec->AddParticle(
231  pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
232  ppi3.Px(), ppi3.Py(), ppi3.Pz(), Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
233 }
234 //___________________________________________________________________________
236 {
237  //
238  // This method generates the final state hadronic system (pion + nucleus) in
239  // COH interactions
240  //
242 
243  Interaction * interaction = evrec->Summary();
244  const XclsTag & xcls_tag = interaction->ExclTag();
245 
246  //-- Access neutrino, initial nucleus and final state prim. lepton entries
247  GHepParticle * nu = evrec->Probe();
248  GHepParticle * Ni = evrec->TargetNucleus();
249  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
250  assert(nu);
251  assert(Ni);
252  assert(fsl);
253 
254  const TLorentzVector & vtx = *(nu->X4());
255  const TLorentzVector & p4nu = *(nu ->P4());
256  const TLorentzVector & p4fsl = *(fsl->P4());
257 
258  //-- Determine the pdg code of the final state pion & nucleus
259  int nucl_pdgc = Ni->Pdg(); // same as the initial nucleus
260  int pion_pdgc = getPionPDGCodeFromXclTag(xcls_tag);
261 
262  //-- basic kinematic inputs
263  double E = nu->E();
264  double M = kNucleonMass;
265  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
266  double mpi2 = TMath::Power(mpi,2);
267  double xo = interaction->Kine().x(true);
268  double yo = interaction->Kine().y(true);
269  double to = interaction->Kine().t(true);
270 
271  SLOG("COHHadronicVtx", pINFO)
272  << "Ev = "<< E << ", xo = " << xo
273  << ", yo = " << yo << ", to = " << to;
274 
275  //-- compute pion energy and |momentum|
276  double Epi = yo * E;
277  double Epi2 = TMath::Power(Epi,2);
278  double ppi2 = Epi2-mpi2;
279  double ppi = TMath::Sqrt(TMath::Max(0.,ppi2));
280 
281  SLOG("COHHadronicVtx", pINFO)
282  << "f/s pion E = " << Epi << ", |p| = " << ppi;
283  assert(Epi>mpi);
284 
285  //-- 4-momentum transfer q=p(neutrino) - p(f/s lepton)
286  // Note: m^2 = q^2 < 0
287  // Also, since the nucleus is heavy all energy loss is comminicated to
288  // the outgoing pion Rein & Sehgal, Nucl.Phys.B223.29-44(1983), p.35
289 
290  TLorentzVector q = p4nu - p4fsl;
291 
292  SLOG("COHHadronicVtx", pINFO)
293  << "\n 4-p transfer q @ LAB: " << utils::print::P4AsString(&q);
294 
295  //-- find angle theta between q and ppi (xi=costheta)
296  // note: t=|(ppi-q)^2|, Rein & Sehgal, Nucl.Phys.B223.29-44(1983), p.36
297 
298  double xi = 1. + M*xo/Epi - 0.5*mpi2/Epi2 - 0.5*to/Epi2;
299  xi /= TMath::Sqrt((1.+2.*M*xo/Epi)*(1.-mpi2/Epi2));
300 
301  double costheta = xi;
302  double sintheta = TMath::Sqrt(TMath::Max(0.,1.-xi*xi));
303 
304  SLOG("COHHadronicVtx", pINFO) << "cos(pion, q) = " << costheta;
305 
306  // compute transverse and longitudinal ppi components along q
307  double ppiL = ppi*costheta;
308  double ppiT = ppi*sintheta;
309 
310  double phi = 2*kPi* rnd->RndHadro().Rndm();
311 
312  TVector3 ppi3(0,ppiT,ppiL);
313 
314  ppi3.RotateZ(phi); // randomize transverse components
315  ppi3.RotateUz(q.Vect().Unit()); // align longit. component with q in LAB
316 
317  SLOG("COHHadronicVtx", pINFO)
318  << "Pion 3-p @ LAB: " << utils::print::Vec3AsString(&ppi3);
319 
320  // now figure out the f/s nucleus 4-p
321 
322  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
323  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
324  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
325  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
326 
327  //-- Save the particles at the GHEP record
328 
329  int mom = evrec->TargetNucleusPosition();
330 
331  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
332  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
333 
334  evrec->AddParticle(pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
335  ppi3.Px(), ppi3.Py(),ppi3.Pz(),Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
336 }
337 //___________________________________________________________________________
339 {
340  Interaction * interaction = evrec->Summary();
341  const Kinematics & kinematics = interaction -> Kine();
342  GHepParticle * nu = evrec->Probe();
343  GHepParticle * Ni = evrec->TargetNucleus();
344  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
345 
346  // Pion
347  const TLorentzVector ppi = kinematics.HadSystP4();
348  const TVector3 ppi3 = ppi.Vect();
349  const double Epi = ppi.E();
350  int pion_pdgc=0;
351  if ( interaction->ProcInfo().IsWeakCC() ) {
352  if( nu->Pdg() > 0 ){ // neutrino
353  pion_pdgc = kPdgPiP;
354  }
355  else{ // anti-neutrino
356  pion_pdgc = kPdgPiM;
357  }
358  }
359  else if ( interaction->ProcInfo().IsWeakNC() ) {
360  pion_pdgc = kPdgPi0;
361  }
362  else{
363  LOG("COHHadronicSystemGeneratorAR", pFATAL)
364  << "Could not determine pion involved in interaction";
365  exit(1);
366  }
367 
368  //
369  // Nucleus
370  int nucl_pdgc = Ni->Pdg(); // pdg of final nucleus same as the initial nucleus
371  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - ppi3.Px();
372  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - ppi3.Py();
373  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - ppi3.Pz();
374  double ENf = nu->E() + Ni->E() - fsl->E() - Epi;
375  //
376  // Both
377  const TLorentzVector & vtx = *(nu->X4());
378  int mom = evrec->TargetNucleusPosition();
379 
380  //
381  // Fill the records
382  evrec->AddParticle(nucl_pdgc,kIStStableFinalState, mom,-1,-1,-1,
383  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
384 
385  evrec->AddParticle(pion_pdgc,kIStStableFinalState, mom,-1,-1,-1,
386  ppi3.Px(), ppi3.Py(),ppi3.Pz(),Epi, vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
387 }
388 
Cross Section Calculation Interface.
const double kPi
Basic constants.
bool IsWeakCC(void) const
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
void CalculateHadronicSystem_AlvarezRuso(GHepRecord *event_rec) const
#define pERROR
Definition: Messenger.h:60
double E(void) const
Get energy.
Definition: GHepParticle.h:92
static const double kNucleonMass
Definition: Constants.h:78
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int getPionPDGCodeFromXclTag(const XclsTag &xcls_tag) const
virtual Interaction * Summary(void) const
Definition: GHepRecord.cxx:135
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:79
const TLorentzVector * P4(void) const
Definition: GHepParticle.h:79
T sqrt(T number)
Definition: d0nt_math.hpp:156
#define pFATAL
Definition: Messenger.h:57
Defines the EventGeneratorI interface.
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
int NPiMinus(void) const
Definition: XclsTag.h:58
string P4AsString(const TLorentzVector *p)
Definition: PrintUtils.cxx:34
void CalculateHadronicSystem_ReinSehgal(GHepRecord *event_rec) const
double x(bool selected=false) const
Definition: Kinematics.cxx:109
const TLorentzVector & HadSystP4(void) const
Definition: Kinematics.h:67
Definition: config.py:1
double Mass(void) const
Definition: Target.cxx:241
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:30
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:91
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:37
double Px(void) const
Get Px.
Definition: GHepParticle.h:89
virtual GHepParticle * Probe(void) const
Definition: GHepRecord.cxx:321
double y(bool selected=false) const
Definition: Kinematics.cxx:122
int Pdg(void) const
Definition: GHepParticle.h:64
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:97
string Name(void) const
Definition: AlgId.h:45
Float_t E
Definition: plot.C:20
void CalculateHadronicSystem_BergerSehgalFM(GHepRecord *event_rec) const
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:370
const Kinematics & Kine(void) const
Definition: Interaction.h:71
int NPiPlus(void) const
Definition: XclsTag.h:57
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
virtual GHepParticle * TargetNucleus(void) const
Definition: GHepRecord.cxx:330
#define pINFO
Definition: Messenger.h:63
int NPi0(void) const
Definition: XclsTag.h:56
virtual const XSecAlgorithmI * CrossSectionAlg(void) const =0
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:54
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:80
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
TParticlePDG * Find(int pdgc)
Definition: PDGLibrary.cxx:61
exit(0)
assert(nhit_max >=nhit_nbins)
const int kPdgPiM
Definition: PDGCodes.h:136
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:535
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double t(bool selected=false) const
Definition: Kinematics.cxx:180
Abstract class. Is used to pass some commonly recurring methods to all concrete implementations of th...
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
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:85
void kinematics()
Definition: kinematics.C:10
GENIE&#39;s GHEP MC event record.
Definition: GHepRecord.h:46
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:40
string Vec3AsString(const TVector3 *vec)
Definition: PrintUtils.cxx:87
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:406
Initial State information.
Definition: InitialState.h:49
double Py(void) const
Get Py.
Definition: GHepParticle.h:90