Functions | Variables
gtestHadronization.cxx File Reference
#include <sstream>
#include <cassert>
#include <cstdlib>
#include <string>
#include <vector>
#include <RVersion.h>
#include <TFile.h>
#include <TDirectory.h>
#include <TTree.h>
#include <TH1F.h>
#include <TLorentzVector.h>
#include <TMath.h>
#include <TClonesArray.h>
#include <TIterator.h>
#include "Framework/Algorithm/Algorithm.h"
#include "Framework/Algorithm/AlgFactory.h"
#include "Physics/Hadronization/HadronizationModelI.h"
#include "Framework/GHEP/GHepStatus.h"
#include "Framework/GHEP/GHepParticle.h"
#include "Framework/GHEP/GHepRecord.h"
#include "Framework/Interaction/ProcessInfo.h"
#include "Framework/Interaction/InitialState.h"
#include "Framework/Interaction/Interaction.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Physics/Hadronization/FragmRecUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"

Go to the source code of this file.

Functions

void pysphe_ (double *, double *)
 
void pythru_ (double *, double *)
 
void PrintSyntax (void)
 
void GetCommandLineArgs (int argc, char **argv)
 
void FillQrkArray (InteractionType_t it, int nu, int *QrkCode, bool *SeaQrk, int nmax, int &nqrk)
 
int main (int argc, char **argv)
 
hadnt Write ("hadnt")
 
gOutFile Close ()
 

Variables

TFile * gOutFile = 0
 
int gNEvents = -1
 
bool gSetHitQrk = false
 
string gHadAlg = ""
 
string gHadConfig = ""
 
 return
 

Function Documentation

gOutFile Close ( )
void FillQrkArray ( InteractionType_t  it,
int  nu,
int *  QrkCode,
bool *  SeaQrk,
int  nmax,
int &  nqrk 
)

Definition at line 338 of file gtestHadronization.cxx.

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 }
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
genie::kPdgUQuark
const int kPdgUQuark
Definition: PDGCodes.h:42
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
genie::kPdgDQuark
const int kPdgDQuark
Definition: PDGCodes.h:44
void GetCommandLineArgs ( int  argc,
char **  argv 
)

Definition at line 372 of file gtestHadronization.cxx.

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 }
gHadConfig
string gHadConfig
Definition: gtestHadronization.cxx:79
gHadAlg
string gHadAlg
Definition: gtestHadronization.cxx:78
pFATAL
#define pFATAL
Definition: Messenger.h:56
freeze_graph.argv
argv
Definition: freeze_graph.py:218
gNEvents
int gNEvents
Definition: gtestHadronization.cxx:76
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pINFO
#define pINFO
Definition: Messenger.h:62
gSetHitQrk
bool gSetHitQrk
Definition: gtestHadronization.cxx:77
genie::CmdLnArgParser
Command line argument parser.
Definition: CmdLnArgParser.h:29
PrintSyntax
void PrintSyntax(void)
Definition: gtestHadronization.cxx:418
int main ( int  argc,
char **  argv 
)

Definition at line 81 of file gtestHadronization.cxx.

82 {
83  GetCommandLineArgs(argc, argv);
84 
85  AlgFactory * algf = AlgFactory::Instance();
86 
87  const HadronizationModelI * model =
88  dynamic_cast<const HadronizationModelI *> (
89  algf->GetAlgorithm(gHadAlg, gHadConfig));
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
gHadConfig
string gHadConfig
Definition: gtestHadronization.cxx:79
thr
G4double thr[100]
Definition: G4S2Light.cc:59
gHadAlg
string gHadAlg
Definition: gtestHadronization.cxx:78
truth_ana.inu
int inu
Definition: truth_ana.py:132
genie::GHepParticle::FirstDaughter
int FirstDaughter(void) const
Definition: GHepParticle.h:68
pysphe_
void pysphe_(double *, double *)
genie::kPdgNuMu
const int kPdgNuMu
Definition: PDGCodes.h:30
model
Definition: model.py:1
genie::GHepParticle::Mass
double Mass(void) const
Mass that corresponds to the PDG code.
Definition: GHepParticle.cxx:143
train.init
init
Definition: train.py:42
genie::GHepParticle::Pz
double Pz(void) const
Get Pz.
Definition: GHepParticle.h:90
genie::GHepParticle::Status
GHepStatus_t Status(void) const
Definition: GHepParticle.h:64
genie::GHepParticle::Energy
double Energy(void) const
Get energy.
Definition: GHepParticle.h:92
freeze_graph.argv
argv
Definition: freeze_graph.py:218
genie::GHepParticle::Px
double Px(void) const
Get Px.
Definition: GHepParticle.h:88
genie::kPdgK0
const int kPdgK0
Definition: PDGCodes.h:174
genie::GHepParticle::Pdg
int Pdg(void) const
Definition: GHepParticle.h:63
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
gNEvents
int gNEvents
Definition: gtestHadronization.cxx:76
genie::GHepParticle::LastDaughter
int LastDaughter(void) const
Definition: GHepParticle.h:69
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::kPdgKM
const int kPdgKM
Definition: PDGCodes.h:173
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
dump_to_simple_cpp.W
W
Definition: dump_to_simple_cpp.py:43
genie::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
genie::kPdgIndep
const int kPdgIndep
Definition: PDGCodes.h:231
genie::kPdgKP
const int kPdgKP
Definition: PDGCodes.h:172
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::kPdgString
const int kPdgString
Definition: PDGCodes.h:230
pINFO
#define pINFO
Definition: Messenger.h:62
gSetHitQrk
bool gSetHitQrk
Definition: gtestHadronization.cxx:77
gOutFile
TFile * gOutFile
Definition: gtestHadronization.cxx:75
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
GetCommandLineArgs
void GetCommandLineArgs(int argc, char **argv)
Definition: gtestHadronization.cxx:372
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::kPdgCluster
const int kPdgCluster
Definition: PDGCodes.h:229
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
genie::AlgFactory
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::InteractionType_t
enum genie::EInteractionType InteractionType_t
genie::GHepParticle
STDHEP-like event record entry that can fit a particle or a nucleus.
Definition: GHepParticle.h:39
pythru_
void pythru_(double *, double *)
FillQrkArray
void FillQrkArray(InteractionType_t it, int nu, int *QrkCode, bool *SeaQrk, int nmax, int &nqrk)
Definition: gtestHadronization.cxx:338
genie::InitialState
Initial State information.
Definition: InitialState.h:48
genie::GHepParticle::Py
double Py(void) const
Get Py.
Definition: GHepParticle.h:89
void PrintSyntax ( void  )

Definition at line 418 of file gtestHadronization.cxx.

419 {
420  LOG("test", pNOTICE)
421  << "\n\n" << "Syntax:" << "\n"
422  << " testHadronization -n nevents -a hadronizer -c config [-q]\n";
423 }
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
pNOTICE
#define pNOTICE
Definition: Messenger.h:61
void pysphe_ ( double *  ,
double *   
)
void pythru_ ( double *  ,
double *   
)
hadnt Write ( "hadnt"  )

Variable Documentation

string gHadAlg = ""

Definition at line 78 of file gtestHadronization.cxx.

string gHadConfig = ""

Definition at line 79 of file gtestHadronization.cxx.

int gNEvents = -1

Definition at line 76 of file gtestHadronization.cxx.

TFile* gOutFile = 0

Definition at line 75 of file gtestHadronization.cxx.

bool gSetHitQrk = false

Definition at line 77 of file gtestHadronization.cxx.

return

Definition at line 335 of file gtestHadronization.cxx.

Generated by   doxygen 1.8.11