GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
DFRHadronicSystemGenerator.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 
28 
29 using namespace genie;
30 using namespace genie::constants;
31 
32 //___________________________________________________________________________
34 HadronicSystemGenerator("genie::DFRHadronicSystemGenerator")
35 {
36 
37 }
38 //___________________________________________________________________________
40 HadronicSystemGenerator("genie::DFRHadronicSystemGenerator", config)
41 {
42 
43 }
44 //___________________________________________________________________________
46 {
47 
48 }
49 //___________________________________________________________________________
51 {
52 // This method generates the final state hadronic system (meson + nucleus) in
53 // diffractive scattering
54 //
56 
57  //-- Access neutrino, initial nucleus and final state prim. lepton entries
58  GHepParticle * nu = evrec->Probe();
59  GHepParticle * Ni = evrec->HitNucleon();
60  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
61  assert(nu);
62  assert(Ni);
63  assert(fsl);
64 
65  Interaction * interaction = evrec->Summary();
66 
67  //-- Determine the pdg code of the final state pion
68  const XclsTag & xcls_tag = interaction->ExclTag();
69  int pion_pdgc = 0;
70  if (xcls_tag.NPi0() == 1) pion_pdgc = kPdgPi0;
71  else if (xcls_tag.NPiPlus() == 1) pion_pdgc = kPdgPiP;
72  else if (xcls_tag.NPiMinus() == 1) pion_pdgc = kPdgPiM;
73  else {
74  LOG("DFRHadronicVtx", pFATAL)
75  << "No final state pion information in XclsTag!";
76  exit(1);
77  }
78 
79  //-- Determine the pdg code of the recoil nucleon
80  int nucl_pdgc = Ni->Pdg(); // same as the hit nucleon
81 
82 
83  const TLorentzVector & p4nu = *(nu ->P4());
84  const TLorentzVector & p4Ni = *(Ni ->P4());
85  const TLorentzVector & p4fsl = *(fsl->P4());
86 
87  // q at LAB
88  TLorentzVector q = p4nu - p4fsl;
89 
90  // q at NRF
91  TLorentzVector qnrf(q);
92  TVector3 beta = p4Ni.BoostVector();
93  qnrf.Boost(-beta);
94 
95  const InitialState & init_state = interaction -> InitState();
96  const Target & target = init_state.Tgt();
97 
98  double E = init_state.ProbeE(kRfHitNucRest); // neutrino energy
99  double M = target.HitNucMass();
100  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
101  double mpi2 = TMath::Power(mpi,2);
102  double xo = interaction->Kine().x(true);
103  double yo = interaction->Kine().y(true);
104  double to = interaction->Kine().t(true);
105  double Epi = qnrf.E() - 0.5*to/M;
106  double Epi2 = TMath::Power(Epi,2);
107  double ppi2 = Epi2-mpi2;
108  double ppi = TMath::Sqrt(TMath::Max(0.,ppi2));
109 
110  SLOG("DFRHadronicVtx", pINFO)
111  << "Ev = "<< E
112  << ", xo = " << xo << ", yo = " << yo << ", to = " << to;
113  SLOG("DFRHadronicVtx", pINFO)
114  << "f/s pion E = " << Epi << ", |p| = " << ppi;
115 
116  // find angle theta between q and ppi (xi=costheta)
117 
118  double xi = -to - mpi2 - qnrf.M2() + 2*qnrf.E()*Epi;
119  xi /= (2 * TMath::Sqrt(Epi2-mpi2) * qnrf.Vect().Mag());
120 
121  if(xi < 0. || xi > 1. || Epi<= mpi) {
122  LOG("DFRKinematics", pWARN) << "Invalid selected kinematics; Attempt regenerating";
123  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
125  exception.SetReason("Invalid selected kinematics");
126  exception.SwitchOnStepBack();
127  exception.SetReturnStep(0);
128  throw exception;
129  }
130 
131  LOG("DFRHadronicVtx", pINFO) << "|to| = " << -to;
132  LOG("DFRHadronicVtx", pINFO) << "q2 = " << qnrf.M2();
133  LOG("DFRHadronicVtx", pINFO) << "xi = " << xi;
134 
135  double costheta = xi;
136  double sintheta = TMath::Sqrt(TMath::Max(0.,1.-xi*xi));
137 
138  SLOG("DFRHadronicVtx", pINFO) << "cos(pion, q) = " << costheta;
139 
140  // compute transverse and longitudinal ppi components along q
141  double ppiL = ppi*costheta;
142  double ppiT = ppi*sintheta;
143  double phi = 2*kPi* rnd->RndHadro().Rndm();
144 
145  TVector3 ppi3(0,ppiT,ppiL); // @ NRF' (NRF rotated so that q along z)
146 
147  ppi3.RotateZ(phi); // randomize transverse components
148  ppi3.RotateUz(qnrf.Vect().Unit()); // align longit. component with q in NRF
149 
150  TLorentzVector p4pi(ppi3,Epi);
151  p4pi.Boost(beta);
152 
153 
154  //-- Now figure out the f/s nucleon 4-p
155 
156  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - p4pi.Px();
157  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - p4pi.Py();
158  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - p4pi.Pz();
159  double ENf = nu->E() + Ni->E() - fsl->E() - p4pi.E();
160 
161  //-- Save the particles at the GHEP record
162 
163  int mom = evrec->HitNucleonPosition();
164  const TLorentzVector & vtx = *(nu->X4());
165 
166  const Target & tgt = interaction->InitState().Tgt();
167  GHepStatus_t ist = (tgt.IsNucleus()) ?
169 
170  evrec->AddParticle(nucl_pdgc, ist, mom,-1,-1,-1,
171  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
172 
173  evrec->AddParticle(pion_pdgc, ist, mom,-1,-1,-1,
174  p4pi.Px(), p4pi.Py(),p4pi.Pz(),p4pi.E(), vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
175 }
176 //___________________________________________________________________________
void ProcessEventRecord(GHepRecord *event_rec) const
double E(void) const
Get energy.
Definition: GHepParticle.h:91
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
double HitNucMass(void) const
Definition: Target.cxx:233
#define pFATAL
Definition: Messenger.h:56
bool IsNucleus(void) const
Definition: Target.cxx:272
int NPiMinus(void) const
Definition: XclsTag.h:60
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:418
double x(bool selected=false) const
Definition: Kinematics.cxx:99
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
An exception thrown by EventRecordVisitorI when the normal processing sequence has to be disrupted (f...
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:333
const Kinematics & Kine(void) const
Definition: Interaction.h:71
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int NPiPlus(void) const
Definition: XclsTag.h:59
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
#define pINFO
Definition: Messenger.h:62
int NPi0(void) const
Definition: XclsTag.h:58
#define pWARN
Definition: Messenger.h:60
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:53
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:117
const TLorentzVector * X4(void) const
Definition: GHepParticle.h:79
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:313
const int kPdgPiM
Definition: PDGCodes.h:159
virtual void AddParticle(const GHepParticle &p)
Definition: GHepRecord.cxx:499
const InitialState & InitState(void) const
Definition: Interaction.h:69
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
const Target & Tgt(void) const
Definition: InitialState.h:66
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:84
double ProbeE(RefFrame_t rf) const
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
enum genie::EGHepStatus GHepStatus_t
Initial State information.
Definition: InitialState.h:48
double Py(void) const
Get Py.
Definition: GHepParticle.h:89