CMS 3D CMS Logo

Cascade2Hadronizer.cc
Go to the documentation of this file.
2 
3 #include "HepMC/GenEvent.h"
4 #include "HepMC/PdfInfo.h"
5 #include "HepMC/PythiaWrapper6_4.h"
7 
8 #include "HepMC/HEPEVT_Wrapper.h"
9 #include "HepMC/IO_HEPEVT.h"
10 #include "HepMC/IO_GenEvent.h"
11 
16 
17 HepMC::IO_HEPEVT hepevtio;
18 
19 #include "HepPID/ParticleIDTranslations.hh"
20 
22 
23 //-- Pythia6 routines and functionalities to pass around Pythia6 params
24 
27 
29 
30 #include "CLHEP/Random/RandomEngine.h"
31 
32 using namespace edm;
33 using namespace std;
34 
35 #define debug 0
36 
37 static CLHEP::HepRandomEngine* cascade2RandomEngine;
38 
39 extern "C" {
40 
41  double dcasrn_(int *idummy) {
42 
43  static int call = 0;
44 
45  double rdm_nb = cascade2RandomEngine->flat();
46 
47  if(debug && ++call < 100) cout<<"dcasrn from c++, call: "<<call<<" random number: "<<rdm_nb<<endl;
48 
49  return rdm_nb;
50  }
51 }
52 
53 
54 namespace gen {
55 
57 
58  public:
60 
61  private:
62  void upInit() override {
63  FortranCallback::getInstance()->fillHeader();
64  }
65 
66  void upEvnt() override {
67  FortranCallback::getInstance()->fillEvent();
68  }
69 
70  bool upVeto() override {
71  bool veto = false;
72  if (!hepeup_.nup) veto = true; //-- LHE Common Blocks
73  return(veto);
74  }
75 
76  };
77 
78  static struct {
79  int n, npad, k[5][pyjets_maxn];
80  double p[5][pyjets_maxn], v[5][pyjets_maxn];
81  } pyjets_local;
82 
83  const std::vector<std::string> Cascade2Hadronizer::theSharedResources = { edm::SharedResourceNames::kPythia6,
85 
86  Cascade2Hadronizer::Cascade2Hadronizer(edm::ParameterSet const& pset)
87  : BaseHadronizer(pset),
88  fPy6Service(new Pythia6ServiceWithCallback(pset)), //-- this will store py6 parameters for further settings
89 
90  fComEnergy(pset.getParameter<double>("comEnergy")),
91  fextCrossSection(pset.getUntrackedParameter<double>("crossSection",-1.)),
92  fextCrossSectionError(pset.getUntrackedParameter<double>("crossSectionError",-1.)),
93  fFilterEfficiency(pset.getUntrackedParameter<double>("filterEfficiency",-1.)),
94 
95  fMaxEventsToPrint(pset.getUntrackedParameter<int>("maxEventsToPrint",0)),
96  fHepMCVerbosity(pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity",false)),
97  fPythiaListVerbosity(pset.getUntrackedParameter<int>("pythiaPylistVerbosity",0)),
98 
99  fDisplayPythiaBanner(pset.getUntrackedParameter<bool>("displayPythiaBanner",false)),
100  fDisplayPythiaCards(pset.getUntrackedParameter<bool>("displayPythiaCards",false)){
101 
102  fParameters = pset.getParameter<ParameterSet>("Cascade2Parameters");
103 
104  fConvertToPDG = false;
105  if(pset.exists("doPDGConvert"))
106  fConvertToPDG = pset.getParameter<bool>("doPDGConvert");
107 
108  //-- silence Pythia6 banner printout, unless display requested
109 
111  if(!call_pygive("MSTU(12)=12345")){
112  throw edm::Exception(edm::errors::Configuration,"PythiaError")
113  <<" Pythia did not accept MSTU(12)=12345";
114  }
115  }
116 
117  //-- silence printouts from PYGIVE, unless display requested
118 
119  if(!fDisplayPythiaCards){
120  if(!call_pygive("MSTU(13)=0")){
121  throw edm::Exception(edm::errors::Configuration,"PythiaError")
122  <<" Pythia did not accept MSTU(13)=0";
123  }
124  }
125 
126  //-- tmp stuff to deal with EvtGen corrupting pyjets
127  flushTmpStorage();
128  }
129 
131  if(fPy6Service != nullptr) delete fPy6Service;
132  }
133 
134  void Cascade2Hadronizer::doSetRandomEngine(CLHEP::HepRandomEngine* v) {
137  }
138 
140 
141  pyjets_local.n = 0 ;
142  pyjets_local.npad = 0 ;
143 
144  for(int ip=0; ip<pyjets_maxn; ip++){
145  for(int i=0; i<5; i++){
146  pyjets_local.k[i][ip] = 0;
147  pyjets_local.p[i][ip] = 0.;
148  pyjets_local.v[i][ip] = 0.;
149  }
150  }
151  return;
152  }
153 
155 
156  pyjets_local.n = pyjets.n;
157  pyjets_local.npad = pyjets.npad;
158 
159  for(int ip=0; ip<pyjets_maxn; ip++){
160  for(int i=0; i<5; i++){
161  pyjets_local.k[i][ip] = pyjets.k[i][ip];
162  pyjets_local.p[i][ip] = pyjets.p[i][ip];
163  pyjets_local.v[i][ip] = pyjets.v[i][ip];
164  }
165  }
166  return ;
167  }
168 
170 
171  HepMC::PdfInfo pdf;
172 
173  //-- filling factorization "Q scale" now! pthat moved to binningValues()
174 
175  if(event()->signal_process_id() <= 0) event()->set_signal_process_id(pypars.msti[0]);
176  if(event()->event_scale() <=0 ) event()->set_event_scale(pypars.pari[22]);
177  if(event()->alphaQED() <= 0) event()->set_alphaQED(pyint1.vint[56]);
178  if(event()->alphaQCD() <= 0) event()->set_alphaQCD(pyint1.vint[57]);
179 
180  //-- get pdf info directly from Pythia6 and set it up into HepMC::GenEvent
181  //-- S. Mrenna: Prefer vint block
182 
183  if(pdf.id1() <= 0) pdf.set_id1(pyint1.mint[14] == 21 ? 0 : pyint1.mint[14]);
184  if(pdf.id2() <= 0) pdf.set_id2(pyint1.mint[15] == 21 ? 0 : pyint1.mint[15]);
185  if(pdf.x1() <= 0) pdf.set_x1(pyint1.vint[40]);
186  if(pdf.x2() <= 0) pdf.set_x2(pyint1.vint[41]);
187  if(pdf.pdf1() <= 0) pdf.set_pdf1(pyint1.vint[38] / pyint1.vint[40]);
188  if(pdf.pdf2() <= 0) pdf.set_pdf2(pyint1.vint[39] / pyint1.vint[41]);
189  if(pdf.scalePDF() <= 0) pdf.set_scalePDF(pyint1.vint[50]);
190 
191  event()->set_pdf_info(pdf) ;
192 
193  if(debug) {
194  cout<<"standard Py6 event weight: pyint1.vint[96]: "<<pyint1.vint[96]<<endl;
195  cout<<"event weight returned by PYEVWT: 1./(pyint1.vint[98]): "<<1./(pyint1.vint[98])<<endl;
196  }
197 
198  //-- this is "standard" Py6 event weight (corresponds to PYINT1/VINT(97))
199  // event()->weights().push_back(pyint1.vint[96]);
200 
201  //-- this is event weight as 1./VINT(99) (PYINT1/VINT(99) is returned by the PYEVWT)
202  // event()->weights().push_back(1./(pyint1.vint[98]));
203 
204  //-- all cascade events have weight = 1
205  event()->weights().push_back(1.);
206 
207  //-- now create the GenEventInfo product from the GenEvent and fill the missing pieces
208 
209  eventInfo().reset(new GenEventInfoProduct(event().get()));
210 
211  //-- in Pythia6 pthat is used to subdivide samples into different bins
212  //-- in LHE mode the binning is done by the external ME generator
213  //-- which is likely not pthat, so only filling it for Py6 internal mode
214 
215  eventInfo()->setBinningValues(vector<double>(1,pypars.pari[16]));
216 
217  //-- here we treat long-lived particles
218 
219  if (pydat1.mstj[21]==3 || pydat1.mstj[21]==4) imposeProperTime();
220 
221  //-- convert particle IDs Py6 -> PDG (if requested)
222 
223  if(fConvertToPDG){
224  for(HepMC::GenEvent::particle_iterator part = event()->particles_begin();
225  part != event()->particles_end(); ++part){
226  (*part)->set_pdg_id(HepPID::translatePythiatoPDT((*part)->pdg_id()));
227  }
228  }
229 
230  //-- service printouts, if requested
231  if(fMaxEventsToPrint > 0){
234  if(fHepMCVerbosity){
235  cout <<"Event process = "<<pypars.msti[0]<<endl<<"----------------------"<<endl;
236  event()->print();
237  }
238  }
239 
240  return;
241  }
242 
244  return false;
245  }
246 
248 
249  //-- grab Py6 instance
251 
253 
254  //-- skip event if py6 considers it bad
255  if(pyint1.mint[50] != 0 ){
256  event().reset();
257  return false;
258  }
259 
260  //-- generation of the event with CASCADE
261  call_event();
262 
263  //-- pythia pyhepc routine converts common PYJETS in common HEPEVT
264  call_pyhepc(1);
265 
266  //-- delete the created event from memory
267  event().reset(hepevtio.read_next_event());
268 
269  //-- this is to deal with post-gen tools & residualDecay() that may reuse PYJETS
270  flushTmpStorage();
271  fillTmpStorage();
272 
273  return true;
274  }
275 
277  return true;
278  }
279 
281 
282  //-- grab Py6 instance
284 
285  int NPartsBeforeDecays = pyjets_local.n ;
286  int NPartsAfterDecays = event().get()->particles_size();
287  int barcode = NPartsAfterDecays;
288 
289  // JVY: well, in principle, it's not a 100% fair to go up to NPartsBeforeDecays,
290  // because Photos will attach gamma's to existing vertices, i.e. in the middle
291  // of the event rather than at the end; but this will only shift pointers down,
292  // so we'll be going again over a few "original" particle...
293  // in the alternative, we may go all the way up to the beginning of the event
294  // and re-check if anything remains to decay, that's fine even if it'll take
295  // some extra CPU...
296 
297  for(int ipart=NPartsAfterDecays; ipart>NPartsBeforeDecays; ipart-- ){
298 
299  HepMC::GenParticle* part = event().get()->barcode_to_particle( ipart );
300  int status = part->status();
301  if ( status != 1 ) continue; // check only "stable" particles,
302  // as some undecayed may still be marked as such
303  int pdgid = part->pdg_id();
304  int py6ptr = pycomp_( pdgid );
305 
306  if(pydat3.mdcy[0][py6ptr-1] != 1) continue; // particle is not expected to decay
307 
308  int py6id = HepPID::translatePDTtoPythia(pdgid);
309 
310  // first, will need to zero out, then fill up PYJETS
311  // I better do it directly (by hands) rather than via py1ent
312  // - to avoid (additional) mass smearing
313 
314  if(part->momentum().t() <= part->generated_mass()){
315  continue ; // e == m -> 0-momentum, nothing to decay...
316  }
317 
318  else{
319  pyjets.n = 0;
320 
321  for(int i=0; i<5; i++){
322  pyjets.k[i][0] = 0;
323  pyjets.p[i][0] = 0.;
324  pyjets.v[i][0] = 0.;
325  }
326 
327  pyjets.k[0][0] = 1;
328  pyjets.k[1][0] = py6id;
329  pyjets.p[4][0] = part->generated_mass();
330  pyjets.p[3][0] = part->momentum().t();
331  pyjets.p[0][0] = part->momentum().x();
332  pyjets.p[1][0] = part->momentum().y();
333  pyjets.p[2][0] = part->momentum().z();
334  HepMC::GenVertex* prod_vtx = part->production_vertex();
335  if ( !prod_vtx ) continue; // in principle, should never happen but...
336  pyjets.v[0][0] = prod_vtx->position().x();
337  pyjets.v[1][0] = prod_vtx->position().y();
338  pyjets.v[2][0] = prod_vtx->position().z();
339  pyjets.v[3][0] = prod_vtx->position().t();
340  pyjets.v[4][0] = 0.;
341  pyjets.n = 1;
342  pyjets.npad = pyjets_local.npad;
343  }
344 
345  //-- now call Py6 decay routine
346 
347  int parent = 1; // since we always pass to Py6 a single particle
348  pydecy_(parent);
349 
350  //-- now attach decay products to mother
351 
352  for(int iprt1=1; iprt1<pyjets.n; iprt1++){
353 
354  part->set_status( 2 );
355 
356  HepMC::GenVertex* DecVtx = new HepMC::GenVertex( HepMC::FourVector(pyjets.v[0][iprt1],
357  pyjets.v[1][iprt1],
358  pyjets.v[2][iprt1],
359  pyjets.v[3][iprt1]) );
360  DecVtx->add_particle_in( part ); // this will cleanup end_vertex if exists, replace with the new one
361  // I presume (vtx) barcode will be given automatically
362 
363  HepMC::FourVector pmom(pyjets.p[0][iprt1],pyjets.p[1][iprt1],
364  pyjets.p[2][iprt1],pyjets.p[3][iprt1] );
365 
366  int dstatus = 0;
367  if(pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10) dstatus = 1;
368  else if(pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20) dstatus = 2;
369  else if(pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30) dstatus = 3;
370  else if(pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100) dstatus = pyjets.k[0][iprt1];
371 
372  HepMC::GenParticle* daughter = new HepMC::GenParticle(pmom,
373  HepPID::translatePythiatoPDT(pyjets.k[1][iprt1]),
374  dstatus);
375  barcode++;
376  daughter->suggest_barcode( barcode );
377  DecVtx->add_particle_out( daughter );
378 
379  int iprt2;
380  for(iprt2=iprt1+1; iprt2<pyjets.n; iprt2++){ // the pointer is shifted by -1, c++ style
381 
382  if(pyjets.k[2][iprt2] != parent){
383  parent = pyjets.k[2][iprt2];
384  break; // another parent particle; reset & break the loop
385  }
386 
387  HepMC::FourVector pmomN(pyjets.p[0][iprt2],pyjets.p[1][iprt2],
388  pyjets.p[2][iprt2],pyjets.p[3][iprt2] );
389 
390  dstatus = 0;
391  if(pyjets.k[0][iprt2] >= 1 && pyjets.k[0][iprt2] <= 10) dstatus = 1;
392  else if(pyjets.k[0][iprt2] >= 11 && pyjets.k[0][iprt2] <= 20) dstatus = 2;
393  else if(pyjets.k[0][iprt2] >= 21 && pyjets.k[0][iprt2] <= 30) dstatus = 3;
394  else if(pyjets.k[0][iprt2] >= 31 && pyjets.k[0][iprt2] <= 100) dstatus = pyjets.k[0][iprt2];
395 
396  HepMC::GenParticle* daughterN = new HepMC::GenParticle(pmomN,
397  HepPID::translatePythiatoPDT( pyjets.k[1][iprt2] ),
398  dstatus);
399  barcode++;
400  daughterN->suggest_barcode( barcode );
401  DecVtx->add_particle_out( daughterN );
402  } //-- end iprt2 loop
403 
404  iprt1 = iprt2-1; // reset counter such that it doesn't go over the same child more than once
405  // don't forget to offset back into c++ counting, as it's already +1 forward
406 
407  event().get()->add_vertex( DecVtx );
408  } //-- end iprt1 loop
409  } //-- end loop over decay products
410 
411  // now restore the very original Py6 event record
412 
413  if(pyjets_local.n != pyjets.n){
414 
415  // restore pyjets to its state as it was before external decays -
416  // might have been jammed by action above or by py1ent calls in EvtGen
417 
418  pyjets.n = pyjets_local.n;
419  pyjets.npad = pyjets_local.npad;
420  for(int ip=0; ip<pyjets_local.n; ip++){
421  for(int i=0; i<5; i++){
422  pyjets.k[i][ip] = pyjets_local.k[i][ip];
423  pyjets.p[i][ip] = pyjets_local.p[i][ip];
424  pyjets.v[i][ip] = pyjets_local.v[i][ip];
425  }
426  }
427  }
428 
429  return true;
430  }
431 
433 
434  //-- grab Py6 instance
436 
438 
439  if (key == 0) fPy6Service->setCSAParams();
441 
442  return true;
443  }
444 
446  return false;
447  }
448 
450 
451  //-- grab Py6 instance
453 
454  //-- change standard parameters of JETSET/PYTHIA - replace call_pytcha()
455 
457 
461 
463 
464  //-- mstu(8) is set to NMXHEP in this dummy call (version >=6.404)
465  call_pyhepc(1);
466 
467  //-- initialise random number generator: has been changed to be CMSSW compliant
468  //-- dcasrn overloaded: call to rluxgo not needed anymore
469  //-- call_rluxgo(4,314159265,0,0);
470 
471  //-- initialise CASCADE parameters (default values)
472  call_casini();
473 
474  //-- Read the parameters and pass them to the common blocks
475  //-- call_steer();
476 
477  //-- initialise CASCADE parameters (user values)
478 
479  //-- retrieve all the different sets
480  vector<string> AllSets = fParameters.getParameter<vector<string> >("parameterSets");
481 
482  //-- loop over the different sets
483  for(unsigned i=0; i<AllSets.size();++i) {
484 
485  string Set = AllSets[i];
486  vector<string> Para_Set = fParameters.getParameter<vector<string> >(Set);
487 
488  //-- loop over all the parameters and stop in case of mistake
489  for(vector<string>::const_iterator itPara = Para_Set.begin(); itPara != Para_Set.end(); ++itPara) {
490 
491  if(!cascadeReadParameters(*itPara)) {
492  throw edm::Exception(edm::errors::Configuration,"CascadeError")
493  <<" Cascade did not accept the following parameter: \""<<*itPara<<"\""<<endl;
494  }
495  } //-- end loop over all the parameters
496  } //-- end loop over the different sets
497 
498  cainpu.plepin = -fComEnergy/2;
499  cainpu.ppin = fComEnergy/2;
500 
502 
503  //-- change standard parameters of CASCADE
504  call_cascha();
505 
506  //-- change standard parameters of JETSET/PYTHIA
507  //-- call_pytcha();
508 
509  //-- set up for running CASCADE (integration of the cross-section)
510  call_cascade();
511 
512  //-- print cross-section result from integration
513  call_caend(1);
514 
516 
518 
519  return true;
520  }
521 
522  bool Cascade2Hadronizer::declareStableParticles(const vector<int>& _pdg){
523  vector<int> pdg = _pdg;
524  for(size_t i=0; i<pdg.size(); i++){
525  int PyID = HepPID::translatePDTtoPythia(pdg[i]);
526  int pyCode = pycomp_( PyID );
527 
528  if(pyCode > 0){
529  ostringstream pyCard ;
530  pyCard << "MDCY(" << pyCode << ",1)=0";
531  //-- cout << "pdg= " << pdg[i] << " " << pyCard.str() << endl;
532 
533  call_pygive(pyCard.str());
534  }
535  }
536 
537  return true;
538  }
539 
541 
542  // this is practically a copy/paste of the original code by J.Alcaraz,
543  // taken directly from PythiaSource
544 
545  int dumm=0;
546  HepMC::GenEvent::vertex_const_iterator vbegin = event()->vertices_begin();
547  HepMC::GenEvent::vertex_const_iterator vend = event()->vertices_end();
548  HepMC::GenEvent::vertex_const_iterator vitr = vbegin;
549 
550  for(; vitr != vend; ++vitr){
551 
552  HepMC::GenVertex::particle_iterator pbegin = (*vitr)->particles_begin(HepMC::children);
553  HepMC::GenVertex::particle_iterator pend = (*vitr)->particles_end(HepMC::children);
554  HepMC::GenVertex::particle_iterator pitr = pbegin;
555 
556  for(; pitr != pend; ++pitr){
557 
558  if((*pitr)->end_vertex()) continue;
559  if((*pitr)->status()!=1) continue;
560 
561  int pdgcode= abs((*pitr)->pdg_id());
562  // Do nothing if the particle is not expected to decay
563  if ( pydat3.mdcy[0][pycomp_(pdgcode)-1] !=1 ) continue;
564 
565  double ctau = pydat2.pmas[3][pycomp_(pdgcode)-1];
566  HepMC::FourVector mom = (*pitr)->momentum();
567  HepMC::FourVector vin = (*vitr)->position();
568  double x = 0.;
569  double y = 0.;
570  double z = 0.;
571  double t = 0.;
572  bool decayInRange = false;
573  while(!decayInRange){
574 
575  double unif_rand = fPy6Service->call(pyr_, &dumm);
576  // Value of 0 is excluded, so following line is OK
577  double proper_length = - ctau * log(unif_rand);
578  double factor = proper_length/mom.m();
579  x = vin.x() + factor * mom.px();
580  y = vin.y() + factor * mom.py();
581  z = vin.z() + factor * mom.pz();
582  t = vin.t() + factor * mom.e();
583 
584  // Decay must be happen outside a cylindrical region
585  if(pydat1.mstj[21]==4){
586  if (sqrt(x*x+y*y)>pydat1.parj[72] || fabs(z)>pydat1.parj[73]) decayInRange = true;
587  // Decay must be happen outside a given sphere
588  }
589  else if (pydat1.mstj[21]==3){
590  if (sqrt(x*x+y*y+z*z)>pydat1.parj[71]) decayInRange = true;
591  }
592  // Decay is always OK otherwise
593  else {
594  decayInRange = true;
595  }
596  }
597 
598  HepMC::GenVertex* vdec = new HepMC::GenVertex(HepMC::FourVector(x,y,z,t));
599  event()->add_vertex(vdec);
600  vdec->add_particle_in((*pitr));
601  }
602  }
603 
604  return;
605  }
606 
608 
610 
612 
613  if(!runInfo().internalXSec()){
614  double sigma_CMS = 0.001 * caeffic.avgi;
615  double error_CMS = 0.001 * caeffic.sd;
616  runInfo().setInternalXSec(GenRunInfoProduct::XSec(sigma_CMS,error_CMS));
617  }
618 
619  //-- print out generated event summary
620  call_caend(2);
621 
622  //-- write out some information from Pythia to the screen
623  call_pystat(1);
624 
625  return;
626  }
627 
628  const char* Cascade2Hadronizer::classname() const {
629  return "gen::Cascade2Hadronizer";
630  }
631 
632  //-- Read the parameters and pass them to the common blocks
633 
634  bool Cascade2Hadronizer::cascadeReadParameters(const string& ParameterString) {
635 
636  bool accepted = true;
637 
638  if(!strncmp(ParameterString.c_str(),"KE",2))
639  caluco.ke = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
640 
641  else if(!strncmp(ParameterString.c_str(),"IRES(1)",7))
642  capar6.ires[0] = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
643 
644  else if(!strncmp(ParameterString.c_str(),"KP",2))
645  caluco.kp = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
646 
647  else if(!strncmp(ParameterString.c_str(),"IRES(2)",7))
648  capar6.ires[1] = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
649 
650  else if(!strncmp(ParameterString.c_str(),"NFRAG",5))
651  cainpu.nfrag = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
652 
653  else if(!strncmp(ParameterString.c_str(),"IPST",4))
654  cashower.ipst = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
655 
656  else if(!strncmp(ParameterString.c_str(),"IPSIPOL",7))
657  jpsi.ipsipol = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
658 
659  //-- from version 2.2.03 on
660  // else if(!strncmp(ParameterString.c_str(),"I23S",4))
661  // jpsi.i23s = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
662 
663  else if(!strncmp(ParameterString.c_str(),"IFPS",4))
664  cainpu.ifps = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
665 
666  else if(!strncmp(ParameterString.c_str(),"ITIMSHR",7))
667  casshwr.itimshr = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
668 
669  else if(!strncmp(ParameterString.c_str(),"IRUNAEM",7))
670  capar1.irunaem = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
671 
672  else if(!strncmp(ParameterString.c_str(),"IRUNA",5))
673  capar1.iruna = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
674 
675  else if(!strncmp(ParameterString.c_str(),"IQ2",3))
676  capar1.iq2 = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
677 
678  else if(!strncmp(ParameterString.c_str(),"IPRO",4))
679  capar1.ipro = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
680 
681  else if(!strncmp(ParameterString.c_str(),"NFLAV",5))
682  caluco.nflav = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
683 
684  else if(!strncmp(ParameterString.c_str(),"INTER",5))
685  cainpu.inter = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
686 
687  else if(!strncmp(ParameterString.c_str(),"IHFLA",5))
688  cahflav.ihfla = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
689 
690  else if(!strncmp(ParameterString.c_str(),"IRPA",4))
691  cascol.irpa = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
692 
693  else if(!strncmp(ParameterString.c_str(),"IRPB",4))
694  cascol.irpb = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
695 
696  else if(!strncmp(ParameterString.c_str(),"IRPC",4))
697  cascol.irpc = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
698 
699  else if(!strncmp(ParameterString.c_str(),"ICCFM",5))
700  casshwr.iccfm = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
701 
702  else if(!strncmp(ParameterString.c_str(),"IFINAL",6))
703  cainpu.ifinal = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
704 
705  else if(!strncmp(ParameterString.c_str(),"IGLU",4))
706  cagluon.iglu = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
707 
708  else if(!strncmp(ParameterString.c_str(),"IRspl",5))
709  casprre.irspl = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
710 
711  else if(!strncmp(ParameterString.c_str(),"PT2CUT",6))
712  captcut.pt2cut[capar1.ipro-1] = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
713 
714  else if(!strncmp(ParameterString.c_str(),"NCB",3))
715  integr.ncb = atoi(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
716 
717  else if(!strncmp(ParameterString.c_str(),"ACC1",4))
718  integr.acc1 = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
719 
720  else if(!strncmp(ParameterString.c_str(),"ACC2",4))
721  integr.acc2 = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
722 
723  else if(!strncmp(ParameterString.c_str(),"SCALFAF",7))
724  scalf.scalfaf = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
725 
726  else if(!strncmp(ParameterString.c_str(),"SCALFA",6))
727  scalf.scalfa = atof(&ParameterString[strcspn(ParameterString.c_str(),"=")+1]);
728 
729  else accepted = false;
730 
731  return accepted;
732  }
733 
735 
736  cout<<""<<endl;
737  cout<<"----------------------------------"<<endl;
738  cout<<"---- Cascade Parameters ----"<<endl;
739  cout<<"----------------------------------"<<endl;
740  cout<<""<<endl;
741 
742  cout<<"computation switch: "<<endl;
743  cout<<"number of calls per iteration for bases: "<<integr.ncb<<endl;
744  cout<<"relative precision for grid optimisation: "<<integr.acc1<<" %"<<endl;
745  cout<<"relative precision for integration: "<<integr.acc2<<" %"<<endl;
746  cout<<""<<endl;
747 
748  cout<<"kinematics: "<<endl;
749  cout<<"flavour code of beam 1: "<<caluco.ke<<endl;
750  cout<<"direct or resolved particle 1: "<<capar6.ires[0]<<endl;
751  cout<<"pz of incoming beam 1: "<<cainpu.plepin<<" GeV"<<endl;
752  cout<<"flavour code of beam 2: "<<caluco.kp<<endl;
753  cout<<"direct or resolved particle 2: "<<capar6.ires[1]<<endl;
754  cout<<"pz of incoming beam 2: "<<cainpu.ppin<<" GeV"<<endl;
755  cout<<"number of active flavors: "<<caluco.nflav<<endl;
756  cout<<""<<endl;
757 
758  cout<<"hard subprocess selection: "<<endl;
759  cout<<"hard subprocess number (IPRO): "<<capar1.ipro<<endl;
760  cout<<"IPRO = 10, switch to select QCD process g* g* -> q qbar: "<<cascol.irpa<<endl;
761  cout<<"IPRO = 10, switch to select QCD process g* g -> g g: "<<cascol.irpb<<endl;
762  cout<<"IPRO = 10, switch to select QCD process g* q -> g q: "<<cascol.irpc<<endl;
763  cout<<"flavor of heavy quark produced (IPRO = 11, 504, 514): "<<cahflav.ihfla<<endl;
764  // cout<<"select vector meson state: "<<jpsi.i23s<<endl;
765  cout<<"use matrix element including J/psi (Upsilon) polarisation: "<<jpsi.ipsipol<<endl;
766  cout<<"pt2 cut in matrix element for massless partons: "<<captcut.pt2cut[capar1.ipro-1]<<" GeV^2"<<endl;
767  cout<<""<<endl;
768 
769  cout<<"parton shower and fragmentation: "<<endl;
770  cout<<"switch for fragmentation: "<<cainpu.nfrag<<endl;
771  cout<<"switch for parton shower (0 off - 1 initial state - 2 final state - 3 initial & final state): "<<cainpu.ifps<<endl;
772  cout<<"switch for time like parton shower in initial state: "<<casshwr.itimshr<<endl;
773  cout<<"select CCFM (1) or DGLAP (0) evolution: "<<casshwr.iccfm<<endl;
774  cout<<"scale switch for final state parton shower: "<<cainpu.ifinal<<endl;
775  cout<<"scale factor for final state parton shower: "<<scalf.scalfaf<<endl;
776  cout<<"switch for proton remnant treatment: "<<casprre.irspl<<endl;
777  cout<<"keep track of intermediate state in parton shower: "<<cashower.ipst<<endl;
778  cout<<"mode of interaction for e p: "<<cainpu.inter<<endl;
779  cout<<""<<endl;
780 
781  cout<<"pdfs, couplings and scales: "<<endl;
782  cout<<"switch for running alphaem: "<<capar1.irunaem<<endl;
783  cout<<"switch for running alphas: "<<capar1.iruna<<endl;
784  cout<<"scale in alphas: "<<capar1.iq2<<endl;
785  cout<<"scale factor for scale in alphas: "<<scalf.scalfa<<endl;
786  cout<<"unintegrated pdf: "<<cagluon.iglu<<endl;
787  cout<<"path where updf are stored: "<<caspdf.pdfpath<<endl;
788  cout<<""<<endl;
789 
790  cout<<"center of mass energy, cross section, filter efficiency: "<<endl;
791  cout<<"center of mass energy: "<<fComEnergy<<" GeV"<<endl;
792  cout<<"external LO cross section: "<<fextCrossSection<<" +/- "<<fextCrossSectionError<<" pb"<<endl;
793  cout<<"filter efficiency: "<<fFilterEfficiency<<endl;
794  cout<<""<<endl;
795  }
796 
798 
799  cout<<""<<endl;
800  cout<<"----------------------------------"<<endl;
801  cout<<"---- Pythia6 Parameters ----"<<endl;
802  cout<<"----------------------------------"<<endl;
803  cout<<""<<endl;
804 
805  cout<<"charm mass: "<<pydat2.pmas[0][3]<<" GeV (default value = 1.5 GeV)"<<endl;
806  cout<<"bottom mass: "<<pydat2.pmas[0][4]<<" GeV (default value = 4.8 GeV)"<<endl;
807  cout<<"Higgs mass: "<<pydat2.pmas[0][24]<<" GeV (default value = 115 GeV)"<<endl;
808 
809  cout<<"lambda QCD: "<<pydat1.paru[111]<<" GeV (default value = 0.25 GeV)"<<endl;
810 
811  cout<<"alphas behaviour: "<<pydat1.mstu[110]<<" (default value = 1)"<<endl;
812  cout<<"nr of flavours wrt lambda QCD: "<<pydat1.mstu[111]<<" (default value = 5)"<<endl;
813  cout<<"min nr of flavours for alphas: "<<pydat1.mstu[112]<<" (default value = 3)"<<endl;
814  cout<<"max nr of flavours for alphas: "<<pydat1.mstu[113]<<" (default value = 5)"<<endl;
815 
816  cout<<"maximum angle settings: "<<pydat1.mstj[47]<<" (default value = 0)"<<endl;
817 
818  cout<<""<<endl;
819  }
820 
821 } //-- namespace gen
#define cainpu
#define pydat1
T getParameter(std::string const &) const
double dcasrn_(int *idummy)
void call(void(&fn)())
#define captcut
Pythia6Service * fPy6Service
#define cahflav
#define caluco
struct HEPEUP_ hepeup_
bool call_pygive(const std::string &line)
bool exists(std::string const &parameterName) const
checks if a parameter exists
void call_casini()
#define cagluon
#define integr
void setInternalXSec(const XSec &xsec)
#define vend()
Definition: vmac.h:43
double v[5][pyjets_maxn]
std::auto_ptr< HepMC::GenEvent > & event()
void doSetRandomEngine(CLHEP::HepRandomEngine *v) override
static const std::string kPythia6
#define cascol
GenRunInfoProduct & runInfo()
double p[5][pyjets_maxn]
#define caspdf
void call_pylist(int mode)
static struct gen::@586 pyjets_local
Pythia6ServiceWithCallback(const edm::ParameterSet &pset)
edm::ParameterSet fParameters
void call_cascade()
T sqrt(T t)
Definition: SSEVec.h:18
void call_caend(int mode)
#define vbegin()
Definition: vmac.h:36
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
#define scalf
bool declareStableParticles(const std::vector< int > &)
#define pyint1
void pydecy_(int &ip)
void setPYUPDAParams(bool afterPyinit)
const char * classname() const
static CLHEP::HepRandomEngine * cascade2RandomEngine
#define capar6
std::auto_ptr< GenEventInfoProduct > & eventInfo()
bool cascadeReadParameters(const std::string &ParameterString)
int k[5][pyjets_maxn]
static const std::string kFortranInstance
#define capar1
static FortranCallback * getInstance()
int pycomp_(int &)
#define pypars
part
Definition: HCALResponse.h:20
#define cashower
void call_event()
#define casshwr
#define debug
void setRandomEngine(CLHEP::HepRandomEngine *v)
void call_cascha()
return(e1-e2)*(e1-e2)+dp *dp
#define jpsi
HLT enums.
#define casprre
void setExternalXSecNLO(const XSec &xsec)
double pyr_(int *idummy)
void setExternalXSecLO(const XSec &xsec)
#define caeffic
HepMC::IO_HEPEVT hepevtio