Project:/generatorMods/src/GrappaModule.cc

CDF KITS
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
Architecture:	[ i386 ]
Version:	[ 4.10.4 ] [ 4.10.5 ] [ 4.8.4 ] [ 4.8.4l3s ] [ 4.8.5 ] [ 4.9.0 ] [ 4.9.1 ] [ 4.9.1.hpt3 ] [ 4.9.1hpt3 ] [ 4.9.1top1 ] [ 5.0.0 ] [ 5.1.0 ] [ 5.1.0beamonly ] [ 5.1.1 ] [ 5.2.0 ] [ 5.3.0 ] [ 5.3.1 ] [ 5.3.1dsp ] [ 5.3.3 ] [ 5.3.3_nt ] [ 5.3.4 ] [ 6.1.1 ] [ 6.1.1b ] [ 6.1.2 ] [ 6.1.3 ] [ 6.1.4 ] [ 6.1.4int3 ] [ 6.1.4mc ] [ 6.1.4mc_a ] [ 6.1.6 ] [ development ]
001 // Grappa interface module to Pythia.
002 // (GRace At Prton-Prton/Anti-proton)
003 //
004 // created                  02/14/2002    Soushi Tsuno
005 // add TAUOLA               03/15/2002    Soushi Tsuno
006 // upgraded to PYTHIA 6.2   03/26/2002    Soushi Tsuno
007 // add W + 0,1,2 jets       06/14/2002    Soushi Tsuno
008 //
009 
010 #include "generatorMods/GrappaModule.hh"
011 #include "pythia_i/Pythia.hh"
012 #include "stdhep_i/CdfHepevt.hh"
013 
014 //---------------
015 // C++ Headers -- MUST use defects !!!!
016 //---------------
017 #include <iostream>
018 #include <iomanip>
019 #include <math.h>
020 using std::cout;
021 using std::endl;
022 using std::setw;
023 using std::string;
024 
025 #include "AbsEnv/AbsEnv.hh"
026 
027 extern "C" {
028   //  void pyupev_(int*, double*);
029 
030   void upinit_();
031   void upevnt_();
032   void grcgproc_(int*,int*,int*);
033 
034   //  void gtauini_(int*);
035   //  void gtaumain_();
036 }
037 
038 const char* GrappaModule::genId="Grappa";
039 const long GrappaModule::_defaultRandomSeed1 = 922883591;
040 const long GrappaModule::_defaultRandomSeed2 = 109735476;
041 
042 //----------------
043 // Constructors --
044 //----------------
045 GrappaModule::GrappaModule() : 
046   AbsGenModule( GrappaModule::genId, 
047                 "AC++ Grappa module"),
048   _grcecm("CMEnergy",  this, 1960.),
049   _grcpap("Beam",      this, "PAP"),
050   _gptcut("gPtcut",    this,   5.0),
051   _getacut("gEtacut",  this,   3.0),
052   _grconcut("gRconcut",this,   0.1),
053   _hmas("hmass",       this, 120.0),
054   _hwid("hwidth",      this,   6.54e-3),
055   _tmas("tmass",       this, 175.0),
056   _twid("twidth",      this,   1.59),
057   _bmas("bmass",       this,   4.8),
058   _cmas("cmass",       this,   1.5),
059   _ijetflv("JetFlv",   this,0, 8,    1,  0,   1),
060   _igsub("iGsub",      this,0,50,  401,400, 550),
061   _incall("iNcall",    this,0,50,100000,0,9999999),
062   _ibswrt("iBswrt",    this,0,50,    0,  0,  12),
063   _inpfl("iNpfl",      this,0,50,    4,  1,   5),
064   _icoup("iCoup",      this,0,50,    1,  1,   6),
065   _ifact("iFact",      this,0,50,    0,  0,   6),
066   _isdep("iSdep",      this,0,50,    0,  0,   1),
067   _grcq("GrcQ",        this,0,50,  0.0,0.0,9999.0),
068   _grcfaq("GrcfaQ",    this,0,50,  0.0,0.0,9999.0),
069   _grcfile("GrcFile",  this,0,50,"bases"),
070 //*  _useTAUOLA("useTAUOLA",this,false),
071   _pythiaMenu( new PythiaMenu( this, 0, "PythiaMenu") ),
072   _randomNumberMenu(),
073   _randomSeed1("RandomSeed1",this,GrappaModule::_defaultRandomSeed1),
074   _randomSeed2("RandomSeed2",this,GrappaModule::_defaultRandomSeed2)
075 {
076   _initializeTalkTo();
077 }
078 
079 //--------------
080 // Destructor --
081 //--------------
082 
083 GrappaModule::~GrappaModule()
084 {
085   delete _pythiaMenu;
086 }
087 
088 AppResult GrappaModule::genBeginRun(AbsEvent* aRun) {
089   return AppResult::OK;
090 }
091 
092 AppResult GrappaModule::genBeginJob() {
093 
094   const int ProcNum = 33;
095 
096   int iGrcSetParm[500],iGrcSetJet[8];
097   double rGrcSetParm[500];
098   char cGrcSetParm[ProcNum][60];
099   double rGrcSetCut[8],rGrcSetMass[8];
100 
101   const char cGrcPar[9][9] = {
102     "iNcall","iBswrt","iNpfl","iCoup","iFact","iSdep",
103     "GrcQ","GrcfaQ","GrcFile"
104   };
105 
106   //============
107   // Initialize
108   //============
109   for (int i = 0; i < 8; ++i) { 
110     rGrcSetCut[i]  = 0.0;
111     rGrcSetMass[i] = 0.0;
112     iGrcSetJet[i] = 0;
113   }
114 
115   //======================
116   //  Default Parameters
117   //======================
118   const int iGrcParDef[ProcNum][7] = {
119     {400,  300000, 0, 4, 5, 0, 0},
120     {401,   10000, 0, 4, 5, 0, 0},
121     {402, 1000000, 0, 4, 5, 0, 0},
122     {403,  200000, 0, 4, 5, 0, 0},
123     {404,   20000, 0, 4, 5, 0, 0},
124     {405,  100000, 0, 4, 2, 0, 0},
125     {406,   10000, 0, 4, 2, 0, 0},
126     {407,  100000, 0, 4, 5, 0, 0},
127     {408,   10000, 0, 4, 5, 0, 0},
128     {409,  500000, 0, 4, 5, 0, 0},
129     {410,  200000, 0, 4, 5, 0, 1},
130     {411,   10000, 0, 4, 5, 0, 1},
131     {412,   10000, 0, 4, 5, 0, 0},
132     {413,  300000, 0, 4, 5, 0, 0},
133     {414,  200000, 0, 4, 5, 0, 1},
134     {415,   10000, 0, 4, 5, 0, 1},
135     {416,  100000, 0, 4, 5, 0, 0},
136     {417,   50000, 0, 4, 5, 0, 0},
137     {418,  500000, 0, 4, 5, 0, 0},
138     {419,   10000, 0, 4, 6, 0, 0},
139     {420,   10000, 0, 4, 6, 0, 0},
140     {421,    5000, 0, 5, 5, 0, 0},
141     {422,   50000, 0, 5, 5, 0, 0},
142     {423,   50000, 0, 5, 5, 0, 0},
143     {424,  200000, 0, 5, 5, 0, 0},
144     {425, 1000000, 0, 5, 5, 0, 0},
145     {426,  200000, 0, 5, 5, 0, 0},
146     {427,  300000, 0, 5, 5, 0, 0},
147     {428, 1000000, 0, 5, 5, 0, 0},
148     {429,  200000, 0, 5, 5, 0, 0},
149     {430,  300000, 0, 5, 5, 0, 0},
150     {500,  100000, 0, 5, 5, 0, 0},
151     {501, 1000000, 0, 5, 5, 0, 0}};
152 
153   const double rGrcParDef[ProcNum][3] = {
154     {400, 120.0,   0},
155     {401, 120.0,   0},
156     {402,  91.188, 0},
157     {403,  91.188, 0},
158     {404, sqrt(120.0*120.0 + 91.188*91.188),   0},
159     {405,   0.0,   0},
160     {406,   0.0,   0},
161     {407, sqrt(2.0)*91.188, 0},
162     {408,  80.45,  0},
163     {409, sqrt(175.0*175.0+80.45*80.45),  0},
164     {410,  91.188, 0},
165     {411,  91.188, 0},
166     {412,  80.45,  0},
167     {413,  80.45,  0},
168     {414,  91.188, 0},
169     {415,  91.188, 0},
170     {416,  91.188, 0},
171     {417,  91.188, 0},
172     {418,  80.45,  0},
173     {419,   0.0,   0},
174     {420,   0.0,   0},
175     {421,  80.45,  0},
176     {422,  80.45,  0},
177     {423,  80.45,  0},
178     {424,  80.45,  0},
179     {425,  80.45,  0},
180     {426,  80.45,  0},
181     {427,  80.45,  0},
182     {428,  80.45,  0},
183     {429,  80.45,  0},
184     {430,  80.45,  0},
185     {500,  80.45,  0},
186     {501,  80.45,  0}};
187       
188   const char cGrcParDef[ProcNum][10] 
189     = {"bases_400","bases_401","bases_402","bases_403","bases_404",
190        "bases_405","bases_406","bases_407","bases_408","bases_409",
191        "bases_410","bases_411","bases_412","bases_413","bases_414",
192        "bases_415","bases_416","bases_417","bases_418","bases_419",
193        "bases_420","bases_421","bases_422","bases_423","bases_424",
194        "bases_425","bases_426","bases_427","bases_428","bases_429",
195        "bases_430","bases_w1j","bases_w2j"};
196 
197   //=================================
198   //  Set CM Energy, Beam Particles
199   //=================================
200   double GrcEcmval = _grcecm.value();
201   const string& GrcBeam = _grcpap.value();
202   rGrcSetParm[0] = GrcEcmval;
203 
204   rGrcSetCut[0] = _gptcut.value();
205   rGrcSetCut[1] = _getacut.value();
206   rGrcSetCut[2] = _grconcut.value();
207 
208   rGrcSetMass[0] = _hmas.value();
209   rGrcSetMass[1] = _hwid.value();
210   rGrcSetMass[2] = _tmas.value();
211   rGrcSetMass[3] = _twid.value();
212   rGrcSetMass[4] = _bmas.value();
213   rGrcSetMass[5] = _cmas.value();
214 
215   //========================
216   //  Set input parameters
217   //========================
218   AbsParmList<int>::ConstIterator igsub  = _igsub.begin();
219   AbsParmList<int>::ConstIterator incall = _incall.begin();
220   AbsParmList<int>::ConstIterator ibswrt = _ibswrt.begin();
221   AbsParmList<int>::ConstIterator inpfl  = _inpfl.begin();
222   AbsParmList<int>::ConstIterator icoup  = _icoup.begin();
223   AbsParmList<int>::ConstIterator ifact  = _ifact.begin();
224   AbsParmList<int>::ConstIterator isdep  = _isdep.begin();
225 
226   AbsParmList<double>::ConstIterator grcq     = _grcq.begin();
227   AbsParmList<double>::ConstIterator grcfaq   = _grcfaq.begin();
228 
229   int nproc = 0;
230   for (AbsParmList<int>::ConstIterator it = _igsub.begin(); 
231        it != _igsub.end(); ++it) {
232     ++nproc;
233   }
234 
235   for (int i = 0; i < nproc; ++i, ++igsub, ++incall, ++ibswrt, 
236          ++inpfl, ++icoup, ++ifact, ++isdep, ++grcq, ++grcfaq) {
237     iGrcSetParm[i]         = *igsub;
238     iGrcSetParm[i+nproc]   = *incall;
239     iGrcSetParm[i+2*nproc] = *ibswrt;
240     iGrcSetParm[i+3*nproc] = *inpfl;
241     iGrcSetParm[i+4*nproc] = *icoup;
242     iGrcSetParm[i+5*nproc] = *ifact;
243     iGrcSetParm[i+6*nproc] = *isdep;
244 
245     rGrcSetParm[i+1]         = *grcq;
246     rGrcSetParm[i+nproc+1]   = *grcfaq;
247   }
248 
249   int i = 0;
250   for (AbsParmList<std::string>::ConstIterator grcfile = _grcfile.begin();
251        grcfile != _grcfile.end(); ++grcfile, ++i) {
252     std::string tmpbsfile = *grcfile;
253     strcpy(cGrcSetParm[iGrcSetParm[i]-400],tmpbsfile.c_str());
254   }
255 
256   //============
257   //  Debug...
258   //============
259   //Set the size of each parameter.
260   int iGrcParSize[12], isubtmp;
261 
262   iGrcParSize[0]  = _incall.size();
263   iGrcParSize[1]  = _ibswrt.size();
264   iGrcParSize[2]  = _inpfl.size();
265   iGrcParSize[3]  = _icoup.size();
266   iGrcParSize[4]  = _ifact.size();
267   iGrcParSize[5]  = _isdep.size();
268   iGrcParSize[6]  = _grcq.size();
269   iGrcParSize[7]  = _grcfaq.size();
270   iGrcParSize[8]  = _grcfile.size();
271 
272   std::cout << "Nr. of Process " << nproc << std::endl;
273 
274   for (int j = 0; j < 9; ++j) {
275     if (iGrcParSize[j] > nproc) {
276       std::cout << " Grappa::Error: different Nr. of " << cGrcPar[j] 
277                 << " = " << iGrcParSize[j] << std::endl;
278       abort();
279     }
280     else if (iGrcParSize[j] < nproc) {
281       for (int i = iGrcParSize[j]; i < nproc; ++i) {
282 
283         if (iGrcSetParm[0] < 500) isubtmp = iGrcSetParm[i]-400;
284         else                      isubtmp = iGrcSetParm[i]-469;
285 
286         if (j < 6) {
287           std::cout << " Use default parameter of " << cGrcPar[j] 
288                     << " ( = " << iGrcParDef[(isubtmp)][j+1] 
289                     << " ) in process = " << iGrcSetParm[i] << std::endl;
290 
291           cout << iGrcParDef[(isubtmp)][j+1] << " " 
292                << iGrcParDef[isubtmp][j+1] << " " << isubtmp << endl;
293           
294           iGrcSetParm[i+(j+1)*nproc] = iGrcParDef[isubtmp][j+1];
295         }
296         else if (j < 8) {
297           std::cout << " Use default parameter of " << cGrcPar[j] 
298                     << " ( = " << rGrcParDef[isubtmp][j-5]
299                     << " ) in process = " << iGrcSetParm[i] << std::endl;
300           
301           rGrcSetParm[i+(j-6)*nproc+1] = rGrcParDef[isubtmp][j-5];
302         }
303         else if (j == 8) {
304           std::cout << " Use default filename of " << cGrcPar[j] 
305                     << " ( = " << cGrcParDef[isubtmp]
306                     << " ) in process = " << iGrcSetParm[i] << std::endl;
307 
308           strcpy(cGrcSetParm[isubtmp],cGrcParDef[isubtmp]);
309         }
310       }
311     }
312   }
313   // check jet flavor
314   if (_ijetflv.size() > 8) {
315     std::cout << " Grappa: too large of Nr. of jet flavors... "
316               << "Execution stopped." << std::endl;
317     abort();
318   }
319   
320   int itmp = 0;
321   for (AbsParmList<int>::ConstIterator ijetflv = _ijetflv.begin();
322        ijetflv != _ijetflv.end(); ++ijetflv) {
323     iGrcSetJet[itmp] = *ijetflv;
324     ++itmp;
325   }
326 
327   // check NOT arrowed flavor
328   if (iGrcSetJet[5] == 1) {
329     std::cout << " Currently, Grappa does not support t-quark jets." 
330               << " ...proceeded as JetFlv[5] = 0..." << std::endl;
331     iGrcSetJet[5] = 0;
332   }
333 
334   // check process
335   if (nproc > 1) {
336     for (int j = 0; j < nproc; ++j) {
337       if (iGrcSetParm[j] > 499) {
338         std::cout << " Process::" << iGrcSetParm[j] 
339                   << " is only arrowed for single process." << std::endl;
340         std::cout << " Shoud not use other process numbers." << std::endl;
341         abort();
342       }
343     }
344   }
345 
346   //================
347   //  To GR@PPA...
348   //================
349   grcgive_(&nproc,iGrcSetParm,iGrcSetJet,rGrcSetParm,rGrcSetCut,rGrcSetMass);
350 
351   for (int i = 0; i < nproc; ++i) {
352     char* bsfile;
353     if (iGrcSetParm[0] < 500) bsfile = cGrcSetParm[iGrcSetParm[i]-400];
354     else                      bsfile = cGrcSetParm[iGrcSetParm[i]-469];
355 
356     int cGrcLen = strlen(bsfile);
357     grcgivec_(bsfile,&cGrcLen,&iGrcSetParm[i]);
358   }
359 
360   std::cout << "OK grcgive" << std::endl;
361 
362   //================
363   //  To Pythia...
364   //================
365   const char* grccol = GrcBeam.c_str();
366 
367   Pythia* pythia = Pythia::Instance();
368 
369   // duplicate some stuff
370 
371   pythia->setEventListLevel( _pythiaMenu->eventlistlevel());
372   pythia->setFirstListEvent( _pythiaMenu->firstlistevent());
373   pythia->setLastListEvent (  _pythiaMenu->lastlistevent());
374 
375   // Read a decay table:
376   const string& fileMode = _pythiaMenu->decayFileMode();
377   const string& fileName = _pythiaMenu->decayFile();
378   if (fileMode=="READ" || fileMode=="read") {
379     if (fileMode.empty()) {
380       ERRLOG(ELerror,"[PYTHIA_NO_DCYTBL]")
381         << "GrappaModule: decay table read selected, but no filename set"
382         << endmsg;
383     } else {
384       int mode=2;
385       int lun=37;
386       if (pythia->Opdcay(fileName.c_str(),&lun,"R")) {
387         ERRLOG(ELerror,"[PYTHIA_NO_DCYTBL]")
388           << "GrappaModule: decay table file " 
389           << "\"" << fileName << "\"\n"
390           << "cannot be opened for read!" << endmsg;
391       } else {
392         pythia->Pyupda(&mode,&lun);
393         pythia->Cldcay(&lun);
394         ERRLOG(ELsuccess,"[PYTHIA_DCYTBL_R]")
395           << "GrappaModule: read decay table file "
396           << "\"" << fileName << "\"" << endmsg;
397       }
398     }
399   }
400   
401   const char* frame = pythia->frame().c_str();
402   const char* beam = pythia->beam().c_str();
403   double winval = _pythiaMenu->win();
404   std::string target = _pythiaMenu->getTarget();
405   int minlist   = _pythiaMenu->initlistlevel();
406   pythia->pysubs().msel() = _pythiaMenu->msel();
407 
408   std::cout << "call pyinit" << std::endl;
409   // Note the `extra' arguments to pyinit on the end of the argument
410   // list which do *not* explicitly occur in the FORTRAN
411   // implementation of the function -- one for each character
412   // variable, in the same order (reference Randy Herber,
413   // herber@fnal.gov). Note also that the `extern "C"' function
414   // declaration matches the C usage rather than the FORTRAN
415   // definition of the PYINIT function and therefore includes the 3
416   // extra arguments.
417 
418   //=========================
419   //  Inintialize of Pythia
420   //=========================
421   //*  std::cout << "Grappa: trying to force PYUPEV " << std::endl;
422   pythia->setPythiaUid("Grappa");
423 
424   // Call Grappa directly.
425   std::cout << "  GrappaModule::genBeginJob : passing cuts and parameters to Grappa ..." << std::endl;
426 
427   //===============
428   //  Integration
429   //===============
430   for (int i = 0; i < nproc; ++i) {
431     int newsub = iGrcSetParm[i];
432     grcgproc_(&nproc,&i,&newsub);
433   }
434 
435   const char* frame62 = "USER";
436   const char* dummy62 = " ";
437   double winval62 = 0.0;
438   pythia->Pyinit(frame62,dummy62,&winval62);
439 
440   pythia -> Pylist(&minlist);
441   // Temporary dummy call to pyofsh to force private version to be used
442   pythia -> Pyofsh();
443   pythia -> setEvents(0);
444   
445   // Write a new decay table:
446   if (fileMode=="WRITE" || fileMode=="write") {
447     if (fileMode.empty()) {
448       ERRLOG(ELerror,"[PYTHIA_NO_DCYTBL]")
449         << "GrappaModule: decay table write selected, but no filename set"
450         << endmsg;
451     } else {
452       int mode=1;
453       int lun=37;
454       if (pythia->Opdcay(fileName.c_str(),&lun,"W")) {
455         ERRLOG(ELerror,"[PYTHIA_NO_DCYTBL]")
456           << "GrappaModule: decay table file " 
457           << "\"" << fileName << "\"\n"
458           << "cannot be opened for write!" << endmsg;
459       } else {
460         pythia->Pyupda(&mode,&lun);
461         pythia->Cldcay(&lun);
462         ERRLOG(ELsuccess,"[PYTHIA_DCYTBL_W]")
463           << "GrappaModule: wrote decay table file "
464           << "\"" << fileName << "\"" << endmsg;
465       }
466     }
467   }
468 
469 
470   //===================
471   // Initialize TAUOLA
472   //===================
473 
474 //*  if (_useTAUOLA.value()) {
475 //*    std::cout << "call gtauini" << std::endl;
476 //*    int newtmp = 100;
477 //*    gtauini_(&newtmp);
478 //*  }
479 
480   std::cout << "End beginjob..." << std::endl;
481 
482   return AppResult::OK;
483 }
484 
485 
486 int GrappaModule::callGenerator(AbsEvent* anEvent) {
487 
488   Pythia* pythia = Pythia::Instance();
489 
490   // Generate the event
491   pythia -> Pyevnt();
492   // Update our event counter 
493   pythia -> addEvent();
494 
495   // See if we are supposed to list this event  
496   int mevlist= pythia->EventListLevel();
497   if ( pythia -> events() >= pythia -> FirstListEvent() && 
498        pythia -> events() <= pythia -> LastListEvent() ) {
499     std::cout << "PYEVNT event no. " << pythia->events() << std::endl;
500     pythia->Pylist(&mevlist); // List this event if required
501   }
502   
503   return 0;
504 }
505 
506 void GrappaModule::fillHepevt() {
507 
508   Pythia* pythia = Pythia::Instance();
509   // Tell lunhep to convert from PYJETS to HEPEVT
510   int mconv=1;                          
511   // Do the conversion
512   //*  pythia->Pylist(&mconv);
513 
514   pythia->Lunhep(&mconv); // Use STDHEP translation routine 1999/01/27 CG
515   CdfHepevt::Instance()->HepevtPtr()->NEVHEP = pythia->events();
516 
517 }
518 
519 AppResult GrappaModule::genEndRun(AbsEvent* aRun) {
520 
521   Pythia* pythia = Pythia::Instance();
522   int mstat = _pythiaMenu->statlistlevel();
523   std::cout << "Call PYSTAT at endRun" << std::endl;
524   pythia->Pystat(&mstat);
525 
526   return AppResult::OK;
527 }
528 
529 AppResult GrappaModule::genEndJob() {
530 
531   for (AbsParmList<int>::ConstIterator it = _ibswrt.begin(); 
532        it != _ibswrt.end(); ++it) {
533     int ibspar = *it;
534     grcendjob_(&ibspar);
535   }
536 
537   Pythia* pythia = Pythia::Instance();
538   int mstat = _pythiaMenu->statlistlevel();
539   std::cout << "Call PYSTAT at endJob" << std::endl;
540   pythia->Pystat(&mstat);
541   
542   return AppResult::OK;
543 }
544
[ source navigation ]
[ diff markup ]
[ identifier search ]
[ freetext search ]
[ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Send problems or questions to cdfcode@fnal.gov