DFRHadronicSystemGenerator.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 - Feb 15, 2009
9 
10  For the class documentation see the corresponding header file.
11 
12  Important revisions after version 2.0.0 :
13  @ Feb 15, 2009 - CA
14  This class was first added in version 2.5.1.
15 
16 */
17 //____________________________________________________________________________
18 
19 #include <cstdlib>
20 
21 #include <TVector3.h>
22 
36 
37 using namespace genie;
38 using namespace genie::constants;
39 
40 //___________________________________________________________________________
42 HadronicSystemGenerator("genie::DFRHadronicSystemGenerator")
43 {
44 
45 }
46 //___________________________________________________________________________
48 HadronicSystemGenerator("genie::DFRHadronicSystemGenerator", config)
49 {
50 
51 }
52 //___________________________________________________________________________
54 {
55 
56 }
57 //___________________________________________________________________________
59 {
60 // This method generates the final state hadronic system (meson + nucleus) in
61 // diffractive scattering
62 //
64 
65  //-- Access neutrino, initial nucleus and final state prim. lepton entries
66  GHepParticle * nu = evrec->Probe();
67  GHepParticle * Ni = evrec->HitNucleon();
68  GHepParticle * fsl = evrec->FinalStatePrimaryLepton();
69  assert(nu);
70  assert(Ni);
71  assert(fsl);
72 
73  Interaction * interaction = evrec->Summary();
74 
75  //-- Determine the pdg code of the final state pion
76  const XclsTag & xcls_tag = interaction->ExclTag();
77  int pion_pdgc = 0;
78  if (xcls_tag.NPi0() == 1) pion_pdgc = kPdgPi0;
79  else if (xcls_tag.NPiPlus() == 1) pion_pdgc = kPdgPiP;
80  else if (xcls_tag.NPiMinus() == 1) pion_pdgc = kPdgPiM;
81  else {
82  LOG("DFRHadronicVtx", pFATAL)
83  << "No final state pion information in XclsTag!";
84  exit(1);
85  }
86 
87  //-- Determine the pdg code of the recoil nucleon
88  int nucl_pdgc = Ni->Pdg(); // same as the hit nucleon
89 
90 
91  const TLorentzVector & p4nu = *(nu ->P4());
92  const TLorentzVector & p4Ni = *(Ni ->P4());
93  const TLorentzVector & p4fsl = *(fsl->P4());
94 
95  // q at LAB
96  TLorentzVector q = p4nu - p4fsl;
97 
98  // q at NRF
99  TLorentzVector qnrf(q);
100  TVector3 beta = p4Ni.BoostVector();
101  qnrf.Boost(-beta);
102 
103  const InitialState & init_state = interaction -> InitState();
104  const Target & target = init_state.Tgt();
105 
106  double E = init_state.ProbeE(kRfHitNucRest); // neutrino energy
107  double M = target.HitNucMass();
108  double mpi = PDGLibrary::Instance()->Find(pion_pdgc)->Mass();
109  double mpi2 = TMath::Power(mpi,2);
110  double xo = interaction->Kine().x(true);
111  double yo = interaction->Kine().y(true);
112  double to = interaction->Kine().t(true);
113  double Epi = qnrf.E() - 0.5*to/M;
114  double Epi2 = TMath::Power(Epi,2);
115  double ppi2 = Epi2-mpi2;
116  double ppi = TMath::Sqrt(TMath::Max(0.,ppi2));
117 
118  SLOG("DFRHadronicVtx", pINFO)
119  << "Ev = "<< E
120  << ", xo = " << xo << ", yo = " << yo << ", to = " << to;
121  SLOG("DFRHadronicVtx", pINFO)
122  << "f/s pion E = " << Epi << ", |p| = " << ppi;
123 
124  // find angle theta between q and ppi (xi=costheta)
125 
126  double xi = -to - mpi2 - qnrf.M2() + 2*qnrf.E()*Epi;
127  xi /= (2 * TMath::Sqrt(Epi2-mpi2) * qnrf.Vect().Mag());
128 
129  if(xi < 0. || xi > 1. || Epi<= mpi) {
130  LOG("DFRKinematics", pWARN) << "Invalid selected kinematics; Attempt regenerating";
131  evrec->EventFlags()->SetBitNumber(kKineGenErr, true);
133  exception.SetReason("Invalid selected kinematics");
134  exception.SwitchOnStepBack();
135  exception.SetReturnStep(0);
136  throw exception;
137  }
138 
139  LOG("DFRHadronicVtx", pINFO) << "|to| = " << -to;
140  LOG("DFRHadronicVtx", pINFO) << "q2 = " << qnrf.M2();
141  LOG("DFRHadronicVtx", pINFO) << "xi = " << xi;
142 
143  double costheta = xi;
144  double sintheta = TMath::Sqrt(TMath::Max(0.,1.-xi*xi));
145 
146  SLOG("DFRHadronicVtx", pINFO) << "cos(pion, q) = " << costheta;
147 
148  // compute transverse and longitudinal ppi components along q
149  double ppiL = ppi*costheta;
150  double ppiT = ppi*sintheta;
151  double phi = 2*kPi* rnd->RndHadro().Rndm();
152 
153  TVector3 ppi3(0,ppiT,ppiL); // @ NRF' (NRF rotated so that q along z)
154 
155  ppi3.RotateZ(phi); // randomize transverse components
156  ppi3.RotateUz(qnrf.Vect().Unit()); // align longit. component with q in NRF
157 
158  TLorentzVector p4pi(ppi3,Epi);
159  p4pi.Boost(beta);
160 
161 
162  //-- Now figure out the f/s nucleon 4-p
163 
164  double pxNf = nu->Px() + Ni->Px() - fsl->Px() - p4pi.Px();
165  double pyNf = nu->Py() + Ni->Py() - fsl->Py() - p4pi.Py();
166  double pzNf = nu->Pz() + Ni->Pz() - fsl->Pz() - p4pi.Pz();
167  double ENf = nu->E() + Ni->E() - fsl->E() - p4pi.E();
168 
169  //-- Save the particles at the GHEP record
170 
171  int mom = evrec->HitNucleonPosition();
172  const TLorentzVector & vtx = *(nu->X4());
173 
174  const Target & tgt = interaction->InitState().Tgt();
175  GHepStatus_t ist = (tgt.IsNucleus()) ?
177 
178  evrec->AddParticle(nucl_pdgc, ist, mom,-1,-1,-1,
179  pxNf, pyNf, pzNf, ENf, 0, 0, 0, 0);
180 
181  evrec->AddParticle(pion_pdgc, ist, mom,-1,-1,-1,
182  p4pi.Px(), p4pi.Py(),p4pi.Pz(),p4pi.E(), vtx.X(), vtx.Y(), vtx.Z(), vtx.T());
183 }
184 //___________________________________________________________________________
185 
const double kPi
Basic constants.
const XML_Char * target
Definition: expat.h:268
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
void ProcessEventRecord(GHepRecord *event_rec) const
double E(void) const
Get energy.
Definition: GHepParticle.h:92
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
enum genie::EGHepStatus GHepStatus_t
double HitNucMass(void) const
Definition: Target.cxx:250
#define pFATAL
Definition: Messenger.h:57
bool IsNucleus(void) const
Definition: Target.cxx:289
int NPiMinus(void) const
Definition: XclsTag.h:58
virtual int HitNucleonPosition(void) const
Definition: GHepRecord.cxx:454
double x(bool selected=false) const
Definition: Kinematics.cxx:109
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
Definition: config.py:1
Double_t beta
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
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:97
Float_t E
Definition: plot.C:20
virtual GHepParticle * FinalStatePrimaryLepton(void) const
Definition: GHepRecord.cxx:370
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:41
int NPiPlus(void) const
Definition: XclsTag.h:57
const int kPdgPiP
Definition: PDGCodes.h:135
const int kPdgPi0
Definition: PDGCodes.h:137
#define pINFO
Definition: Messenger.h:63
int NPi0(void) const
Definition: XclsTag.h:56
#define pWARN
Definition: Messenger.h:61
TRandom3 & RndHadro(void) const
rnd number generator used by hadronization models
Definition: RandomGen.h:54
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:43
virtual TBits * EventFlags(void) const
Definition: GHepRecord.h:118
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
virtual GHepParticle * HitNucleon(void) const
Definition: GHepRecord.cxx:350
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 InitialState & InitState(void) const
Definition: Interaction.h:69
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...
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
double ProbeE(RefFrame_t rf) const
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
Initial State information.
Definition: InitialState.h:49
double Py(void) const
Get Py.
Definition: GHepParticle.h:90