gtestHadronization.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestHadronization
5 
6 \brief Program used for testing/debugging the KNO & PYTHIA hadronizers
7 
8  Syntax :
9  gtestHadronization -n nevents -t test -a hadronizer -c config [-q]
10 
11  Options :
12  -n number of events
13  -a hadronizer (algorithm name, eg genie::AGKYLowW2019)
14  -c hadronizer config set
15  -q set hit quark (needed for PYTHIA, not needed for KNO)
16 
17 \author Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
18  University of Liverpool & STFC Rutherford Appleton Laboratory
19 
20 \created June 20, 2004
21 
22 \cpright Copyright (c) 2003-2020, The GENIE Collaboration
23  For the full text of the license visit http://copyright.genie-mc.org
24 
25 */
26 //____________________________________________________________________________
27 
28 #include <sstream>
29 #include <cassert>
30 #include <cstdlib>
31 #include <string>
32 #include <vector>
33 
34 #include <RVersion.h>
35 #include <TFile.h>
36 #include <TDirectory.h>
37 #include <TTree.h>
38 #include <TH1F.h>
39 #include <TLorentzVector.h>
40 #include <TMath.h>
41 #include <TClonesArray.h>
42 #include <TIterator.h>
43 
46 #include "Physics/Hadronization/HadronizationModelI.h"
56 #include "Physics/Hadronization/FragmRecUtils.h"
58 
59 using std::string;
60 using std::vector;
61 using std::endl;
62 using std::ostringstream;
63 
64 using namespace genie;
65 
66 extern "C" void pysphe_(double *, double *);
67 extern "C" void pythru_(double *, double *);
68 
69 void PrintSyntax (void);
70 void GetCommandLineArgs (int argc, char ** argv);
71 
72 void FillQrkArray(InteractionType_t it, int nu,
73  int * QrkCode, bool * SeaQrk, int nmax, int & nqrk);
74 
75 TFile * gOutFile = 0;
76 int gNEvents = -1;
77 bool gSetHitQrk = false;
78 string gHadAlg = "";
79 string gHadConfig = "";
80 //____________________________________________________________________________
81 int main(int argc, char ** argv)
82 {
83  GetCommandLineArgs(argc, argv);
84 
86 
87  const HadronizationModelI * model =
88  dynamic_cast<const HadronizationModelI *> (
90  assert(model);
91 
92  gOutFile = new TFile("./genie-hadronization.root","recreate");
93 
94  const int npmax = 100; // max number of particles
95 /*
96  const int kCcNc = 2;
97  const int kNNu = 2;
98  const int kNNuc = 2;
99  const int kNQrkMax = 4;
100  const int kNW = 1;
101 
102  int CcNc[kCcNc] = { 1, 2 };
103  int NuCode[kNNu] = { kPdgNuMu, kPdgAntiNuMu };
104  int NucCode[kNNuc] = { kPdgProton, kPdgNeutron };
105  double W[kNW] = { 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 };
106 */
107  const int kCcNc = 1;
108  const int kNNu = 1;
109  const int kNNuc = 1;
110  const int kNQrkMax = 4;
111  const int kNW = 12;
112 
113  int CcNc[kCcNc] = { 1 };
114  int NuCode[kNNu] = { kPdgNuMu };
115  int NucCode[kNNuc] = { kPdgProton };
116 // double W[kNW] = { 1.5, 1.8, 2.0, 2.5, 3.0, 4.0, 5.0 };
117  double W[kNW] = { 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 4.0, 6.0, 8.0, 10.0, 15.0, 20.0 };
118 
119  int QrkCode[kNQrkMax];
120  bool SeaQrk [kNQrkMax];
121 
122  gOutFile->cd();
123 
124  TTree * hadnt = new TTree("hadnt","hadronizer multiplicities");
125 
126  int br_iev; // event number
127  int br_nuc; // hit nucleon PDG code
128  int br_neut; // neutrino PDG code
129  int br_qrk; // hit quark PDG code
130  int br_sea; // hit quark is from the sea=1 (valence=0)
131  int br_ccnc; // CC=1, NC=2
132  float br_W; // hadronic invariant mass
133  int br_np; // number of generated p
134  int br_nn; // number of generated n
135  int br_npip; // number of generated pi+
136  int br_npim; // number of generated pi-
137  int br_npi0; // number of generated pi0
138  int br_nKp; // number of generated K+
139  int br_nKm; // number of generated K-
140  int br_nK0; // number of generated K0
141  int br_n; // total number of generated particles
142  int br_nstrst; // string daughters
143  int br_model; // jetset model info (1:string,2:cluster,3:indep), 0: for kno
144  float br_sphericity; // jetset sphericity
145  float br_aplanarity; // jetset aplanarity
146  float br_thrust; // jetset thrust
147  float br_oblateness; // jetset oblateness
148  int br_pdg[npmax]; // PDG code of each particle
149  int br_ist[npmax]; // Status code of each particle
150  int br_dec[npmax]; // Decay of a previous hadronization product (Delta,...)?
151  float br_px [npmax]; // px of each particle
152  float br_py [npmax]; // py of each particle
153  float br_pz [npmax]; // pz of each particle
154  float br_KE [npmax]; // kinetic energy of each particle
155  float br_E [npmax]; // energy of each particle
156  float br_M [npmax]; // mass of each particle
157  float br_pL [npmax]; // pL of each particle
158  float br_pT2[npmax]; // pT2 of each particle
159  float br_xF [npmax]; // xF of each particle
160  float br_z [npmax]; // xF of each particle
161 
162  hadnt->Branch("iev", &br_iev, "iev/I");
163  hadnt->Branch("nuc", &br_nuc, "nuc/I");
164  hadnt->Branch("neut", &br_neut, "neut/I");
165  hadnt->Branch("qrk", &br_qrk, "qrk/I");
166  hadnt->Branch("sea", &br_sea, "sea/I");
167  hadnt->Branch("ccnc", &br_ccnc, "ccnc/I");
168  hadnt->Branch("W", &br_W, "W/F");
169  hadnt->Branch("np", &br_np, "np/I");
170  hadnt->Branch("nn", &br_nn, "nn/I");
171  hadnt->Branch("npip", &br_npip, "npip/I");
172  hadnt->Branch("npim", &br_npim, "npim/I");
173  hadnt->Branch("npi0", &br_npi0, "npi0/I");
174  hadnt->Branch("nKp", &br_nKp, "nKp/I");
175  hadnt->Branch("nKm", &br_nKm, "nKm/I");
176  hadnt->Branch("nK0", &br_nK0, "nK0/I");
177  hadnt->Branch("n", &br_n, "n/I");
178  hadnt->Branch("jmod", &br_model, "jmod/I");
179  hadnt->Branch("nstrst",&br_nstrst,"nstrst/I");
180  hadnt->Branch("sph", &br_sphericity, "sph/F");
181  hadnt->Branch("apl", &br_aplanarity, "apl/F");
182  hadnt->Branch("thr", &br_thrust, "thr/F");
183  hadnt->Branch("obl", &br_oblateness, "obl/F");
184  hadnt->Branch("pdg", br_pdg, "pdg[n]/I");
185  hadnt->Branch("ist", br_ist, "ist[n]/I");
186  hadnt->Branch("dec", br_dec, "dec[n]/I");
187  hadnt->Branch("px", br_px, "px[n]/F");
188  hadnt->Branch("py", br_py, "py[n]/F");
189  hadnt->Branch("pz", br_pz, "pz[n]/F");
190  hadnt->Branch("KE", br_KE, "KE[n]/F");
191  hadnt->Branch("E", br_E, "E[n]/F");
192  hadnt->Branch("M", br_M, "M[n]/F");
193  hadnt->Branch("pL", br_pL, "pL[n]/F");
194  hadnt->Branch("pT2", br_pT2, "pT2[n]/F");
195  hadnt->Branch("xF", br_xF, "xF[n]/F");
196  hadnt->Branch("z", br_z, "z[n]/F");
197 
198  const int nnull_max=100;
199  int nnull=0;
200 
201  // CC/NC loop
202  for(int iccnc=0; iccnc<kCcNc; iccnc++) {
203  InteractionType_t it = (CcNc[iccnc]==1) ? kIntWeakCC : kIntWeakNC;
204 
205  // neutrino & hit nucleon loops
206  for(int inu=0; inu<kNNu; inu++) {
207  for(int inuc=0; inuc<kNNuc; inuc++) {
208 
209  InitialState init (26,56,NuCode[inu]);
210  ProcessInfo proc (kScDeepInelastic, it);
211  Interaction intr (init, proc);
212 
213  intr.InitStatePtr()->TgtPtr()->SetHitNucPdg(NucCode[inuc]);
214 
215  // hit quark loop (if requested)
216  int nqrk=1;
217  if(gSetHitQrk) {
218  FillQrkArray(it, NuCode[inu], QrkCode, SeaQrk, kNQrkMax, nqrk);
219  }
220  for(int iqrk=0; iqrk<nqrk; iqrk++) {
221  if(gSetHitQrk) {
222  intr.InitStatePtr()->TgtPtr()->SetHitQrkPdg(QrkCode[iqrk]);
223  intr.InitStatePtr()->TgtPtr()->SetHitSeaQrk(SeaQrk[iqrk]);
224  }
225 
226  LOG("test",pNOTICE) << "hadronizing: " << intr.AsString();
227 
228  for(int iw=0; iw<kNW; iw++) {
229  intr.KinePtr()->SetW(W[iw]);
230 
231  for(int in=0; in<gNEvents; in++) {
232  TClonesArray * plist = model->Hadronize(&intr);
233 
234  if(!plist) {
235  // don't count the current event and repeat
236  in--;
237  nnull++;
238  if(nnull>nnull_max) exit(1);
239  continue;
240  }
241 
242  utils::fragmrec::Print(plist);
243 
244  br_iev = in;
245  br_nuc = NucCode[inuc];
246  br_neut = NuCode[inu];
247  br_qrk = ((gSetHitQrk) ? QrkCode[iqrk] : 0);
248  br_sea = ((gSetHitQrk) ? SeaQrk[iqrk] : 0);
249  br_ccnc = CcNc[iccnc];
250  br_W = W[iw];
251  br_np = utils::fragmrec::NParticles(kPdgProton, plist);
252  br_nn = utils::fragmrec::NParticles(kPdgNeutron, plist);
253  br_npip = utils::fragmrec::NParticles(kPdgPiP, plist);
254  br_npim = utils::fragmrec::NParticles(kPdgPiM, plist);
255  br_npi0 = utils::fragmrec::NParticles(kPdgPi0, plist);
256  br_nKp = utils::fragmrec::NParticles(kPdgKP, plist);
257  br_nKm = utils::fragmrec::NParticles(kPdgKM, plist);
258  br_nK0 = utils::fragmrec::NParticles(kPdgK0, plist);
259  br_n = plist->GetEntries();
260 
261  br_model=0;
262  br_nstrst=0;
263 
264  GHepParticle * particle = 0;
265  TIter particle_iter(plist);
266 
267  unsigned int i=0;
268  unsigned int daughter1=0, daughter2=0;
269  bool model_set=false;
270 
271  while( (particle = (GHepParticle *) particle_iter.Next()) ) {
272  br_pdg[i] = particle->Pdg();
273  br_ist[i] = particle->Status();
274  br_px[i] = particle->Px();
275  br_py[i] = particle->Py();
276  br_pz[i] = particle->Pz();
277  br_KE[i] = particle->Energy() - particle->Mass();
278  br_E[i] = particle->Energy();
279  br_M[i] = particle->Mass();
280  br_pL[i] = particle->Pz();
281  br_pT2[i] = TMath::Power(particle->Px(),2) +
282  TMath::Power(particle->Py(),2);
283  br_xF[i] = particle->Pz() / (W[iw]/2);
284  br_z[i] = particle->Energy() / W[iw];
285 
286  if(particle->Pdg() == kPdgString || particle->Pdg() == kPdgCluster || particle->Pdg() == kPdgIndep)
287  if(model_set) exit(1);
288  model_set = true;
289  if (particle->KF() == kPdgString ) br_model=1;
290  else if (particle->KF() == kPdgCluster) br_model=2;
291  else if (particle->KF() == kPdgIndep ) br_model=3;
292 
293  daughter1 = particle->FirstDaughter();
294  daughter2 = particle->LastDaughter();
295  br_nstrst = daughter2-daughter1+1;
296  }
297 
298  if(model_set) {
299  if(i>=daughter1 && i<=daughter2) br_dec[i] = 1;
300  else br_dec[i] = 0;
301  } else {
302  br_dec[i] = 0;
303  }
304 
305  i++;
306  } // particle-iterator
307 
308  double sph=0, apl=0, thr=0, obl=0;
309  if(br_model!=0) {
310  pysphe_(&sph,&apl);
311  pythru_(&thr,&obl);
312  LOG("test", pINFO) << "Sphericity = " << sph << ", aplanarity = " << apl;
313  }
314 
315  br_sphericity = sph;
316  br_aplanarity = apl;
317  br_thrust = thr;
318  br_oblateness = obl;
319 
320  hadnt->Fill();
321 
322  plist->Delete();
323  delete plist;
324  }//n
325  }//w
326  }//qrk
327  }//inuc
328  }//inu
329  }//cc/nc
330 
331  hadnt->Write("hadnt");
332 
333  gOutFile->Close();
334 
335  return 0;
336 }
337 //____________________________________________________________________________
339  int * QrkCode, bool * SeaQrk, int nmax, int & nqrk)
340 {
341 // utility method: create/fill array with all possible hit quarks
342 //
343 
344  for(int i=0; i<nmax; i++) {
345  QrkCode[i] = -1;
346  }
347  if(it==kIntWeakNC) {
348  nqrk=4;
349  assert(nqrk<=nmax);
350  QrkCode[0] = kPdgUQuark; SeaQrk[0] = false;
351  QrkCode[1] = kPdgUQuark; SeaQrk[1] = true;
352  QrkCode[2] = kPdgDQuark; SeaQrk[2] = false;
353  QrkCode[3] = kPdgDQuark; SeaQrk[3] = true;
354  }
355  else if (it==kIntWeakCC) {
356  nqrk=2;
357  assert(nqrk<=nmax);
358  if(pdg::IsNeutrino(nu)) {
359  QrkCode[0] = kPdgDQuark; SeaQrk[0] = false;
360  QrkCode[1] = kPdgDQuark; SeaQrk[1] = true;
361  } else if(pdg::IsAntiNeutrino(nu)) {
362  QrkCode[0] = kPdgUQuark; SeaQrk[0] = false;
363  QrkCode[1] = kPdgUQuark; SeaQrk[1] = true;
364  } else {
365  exit(1);
366  }
367  } else {
368  exit(1);
369  }
370 }
371 //____________________________________________________________________________
372 void GetCommandLineArgs(int argc, char ** argv)
373 {
374 // Parse the command line arguments
375 
376  CmdLnArgParser parser(argc,argv);
377 
378  // number of events:
379  if( parser.OptionExists('n') ) {
380  LOG("test", pINFO) << "Reading number of events to generate";
381  gNEvents = parser.ArgAsInt('n');
382  } else {
383  LOG("test", pFATAL) << "Number of events was not specified";
384  PrintSyntax();
385  exit(1);
386  }
387 
388  // hadronizer:
389  if( parser.OptionExists('a') ) {
390  LOG("test", pINFO) << "Reading hadronization algorithm name";
391  gHadAlg = parser.ArgAsString('a');
392  } else {
393  LOG("test", pINFO) << "No hadronization algorithm was specified";
394  PrintSyntax();
395  exit(1);
396  }
397 
398  // hadronizer config:
399  if( parser.OptionExists('c') ) {
400  LOG("test", pINFO) << "Reading hadronization algorithm config name";
401  gHadConfig = parser.ArgAsString('c');
402  } else {
403  LOG("test", pINFO) << "No hadronization algorithm config was specified";
404  PrintSyntax();
405  exit(1);
406  }
407 
408  // set struck quark?
409  if( parser.OptionExists('q') ) {
410  LOG("test", pINFO) << "reading struck quark option";
411  gSetHitQrk = true;
412  } else {
413  LOG("test", pINFO) << "Using default option for setting hit quark";
414  gSetHitQrk = false;
415  }
416 }
417 //____________________________________________________________________________
418 void PrintSyntax(void)
419 {
420  LOG("test", pNOTICE)
421  << "\n\n" << "Syntax:" << "\n"
422  << " testHadronization -n nevents -a hadronizer -c config [-q]\n";
423 }
424 //____________________________________________________________________________
425 
426 
427 
string gHadConfig
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
string ArgAsString(char opt)
THE MAIN GENIE PROJECT NAMESPACE
Definition: AlgCmp.h:25
G4double thr[100]
Definition: G4S2Light.cc:59
string gHadAlg
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
std::string string
Definition: nybbler.cc:12
int FirstDaughter(void) const
Definition: GHepParticle.h:68
void pysphe_(double *, double *)
#define pFATAL
Definition: Messenger.h:56
const int kPdgUQuark
Definition: PDGCodes.h:42
const int kPdgNuMu
Definition: PDGCodes.h:30
struct vector vector
void SetHitQrkPdg(int pdgc)
Definition: Target.cxx:184
Definition: model.py:1
double Mass(void) const
Mass that corresponds to the PDG code.
init
Definition: train.py:42
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
string AsString(void) const
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
double Energy(void) const
Get energy.
Definition: GHepParticle.h:92
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
const int kPdgK0
Definition: PDGCodes.h:174
int Pdg(void) const
Definition: GHepParticle.h:63
Summary information for an interaction.
Definition: Interaction.h:56
int gNEvents
int LastDaughter(void) const
Definition: GHepParticle.h:69
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const int kPdgKM
Definition: PDGCodes.h:173
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
const int kPdgIndep
Definition: PDGCodes.h:231
const int kPdgKP
Definition: PDGCodes.h:172
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
const int kPdgString
Definition: PDGCodes.h:230
const int kPdgDQuark
Definition: PDGCodes.h:44
#define pINFO
Definition: Messenger.h:62
bool gSetHitQrk
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:279
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
void SetHitNucPdg(int pdgc)
Definition: Target.cxx:171
Target * TgtPtr(void) const
Definition: InitialState.h:67
TFile * gOutFile
const int kPdgPiM
Definition: PDGCodes.h:159
InitialState * InitStatePtr(void) const
Definition: Interaction.h:74
void GetCommandLineArgs(int argc, char **argv)
const int kPdgProton
Definition: PDGCodes.h:81
const int kPdgCluster
Definition: PDGCodes.h:229
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
void SetHitSeaQrk(bool tf)
Definition: Target.cxx:195
const int kPdgNeutron
Definition: PDGCodes.h:83
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
enum genie::EInteractionType InteractionType_t
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
bool OptionExists(char opt)
was option set?
void pythru_(double *, double *)
QTextStream & endl(QTextStream &s)
void FillQrkArray(InteractionType_t it, int nu, int *QrkCode, bool *SeaQrk, int nmax, int &nqrk)
Initial State information.
Definition: InitialState.h:48
void PrintSyntax(void)
int main(int argc, char **argv)
double Py(void) const
Get Py.
Definition: GHepParticle.h:89