CDF II
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 //--------------------------------------------------------------------------
002 // File and Version Information:
003 //
004 // Description:
005 //      Class WGRAD_Module
006 //
007 //      Interface WGRAD, a program for simulating W production with O(alpha) 
008 //      EW radiative corrections. Ref. : hep-ph/9807417
009 //                                       Fermilab-Pub-98/164-T
010 // Author List:
011 //      Dave Waters, Oxford University
012 //
013 //  revision history:
014 //  -----------------
015 //  8.7.99  Dave Waters : first attempt.
016 // 14.8.99  Dave Waters : set talk-to parameters that will be used to
017 //                        run WGRAD directly from simulation.
018 // 17.8.99  Dave Waters : complete the integration of WGRAD. No more
019 //                        need for ascii file input.
020 //  3.9.99  Dave Waters : Pythia interface menu. Choice of bypassing
021 //                        Pythia altogether, with particle 4-momenta
022 //                        loaded into HEPEVT by hand.
023 // aug 29 2001 lena: clean up, remove gEvent. 
024 // aug 30 2001 lena: make include in standard way
025 // jan 14 2001 lena: register as pythia user at beginJob
026 //--------------------------------------------------------------------------
027 
028 //-----------------------
029 // This Class's Header --
030 //-----------------------
031 #include "generatorMods/WGRAD_Module.hh"
032 
033 //---------------
034 // C++ Headers --
035 //---------------
036 #include <iostream>
037 
038 //-------------------------------
039 // Collaborating Class Headers --
040 //-------------------------------
041 #include "BaBar/Cdf.hh"
042 #include "inc/bcs.h"
043 #include "evt/evt.h"
044 //#include "ParticleDB/hepevt.hh"
045 #include "stdhep_i/CdfHepevt.hh"
046 //#include "ParticleDB/CdfParticleDatabase.hh"
047 #include "pythia_i/Pythia.hh"
048 
049 // extern "C" {
050 //  void pyupev_(int*, float*);
051 // }
052 
053 extern "C" {
054   void upinit_(int*, float*);
055   void upevnt_();
056 }
057 
058 //-----------------------------------------------------------------------
059 // Local Macros, Typedefs, Structures, Unions and Forward Declarations --
060 //-----------------------------------------------------------------------
061 #include "r_n/CdfRn.hh"
062 #include <sstream>
063 using std::ostringstream ;
064 
065 const long WGRAD_Module::_defaultRandomSeed1 = 9228;
066 const long WGRAD_Module::_defaultRandomSeed2 = 1097;
067 
068 const char* WGRAD_Module::genId = "WGRAD";
069 
070 //----------------
071 // Constructors --
072 //----------------
073 
074 WGRAD_Module::WGRAD_Module()
075   : AbsGenModule( WGRAD_Module::genId, 
076                   "W production including O(alpha) EW radiative corrections"),
077     _debug("debug",this,false),
078     _usePythia("usePythia",this,true),
079     _useLesHouchesModule("useLesHouchesModule",this,false),
080     _beamType("BeamType",this,1),
081     _W_charge("W_Charge",this,1),
082     _widthTreatment("widthTreatment",this,2),
083     _decayLepton("DecayLepton",this,2),
084     _QED_radType("QED_RadType",this,2),
085     _fsrType("fsrType",this,1),
086     _smearingRecombination("SmearingRecombination",this,0),
087     _order("Order",this,2),
088     _scale("Scale",this,0),
089     _fsrCollinearCut("fsrCollinearCut",this,0),
090     _photonEnergyCut("PhotonEnergyCut",this,0.0001),
091     _photonCollinearCut("PhotonCollinearCut",this,0.0001),
092     _pdflibGroup("pdflibGroup",this,3),
093     _pdflibSet("pdflibSet",this,54),
094     _wMass("W_Mass",this,80.3),
095     _wWidth("W_Width",this,2.1),
096     _CME("CME",this,2000.0),
097     _lepton_ptmin("Lepton_ptmin",this,0.),
098     _lepton_ymax("Lepton_ymax",this,100.),
099     _neutrino_ptmin("Neutrino_ptmin",this,0.),
100     _neutrino_ymax("Neutrino_ymax",this,100.),
101     _photon_ptmin("Photon_ptmin",this,0.),
102     _photon_ymax("Photon_ymax",this,100.),
103     _tmass_low("mT_low",this,0.),
104     _tmass_high("mT_high",this,2000.),
105     _itmaxGrid("itmaxGrid",this,1),
106     _ncallGrid("ncallGrid",this,2000),
107     _itmaxPrecision("itmaxPrecision",this,2),
108     _ncallPrecision("ncallPrecision",this,5000),
109     _rndmSeed("randomSeed",this,4567),
110     // No longer the correct way to setup Pythia : 
111     //_mstj_41("mstj(41)",this,-1),
112     //_mrpy_1("mrpy(1)",this,-1),
113     //_mstp_81("mstp(81)",this,-1),
114     //_mstp_91("mstp(91)",this,-1),
115     //_mstp_111("mstp(111)",this,-1),
116     //_parp_91("parp(91)",this,-1.0),
117     //_parp_92("parp(92)",this,-1.0),
118     //_parp_93("parp(93)",this,-1.0),
119     _randomSeed1("RandomSeed1",this,WGRAD_Module::_defaultRandomSeed1),
120     _randomSeed2("RandomSeed2",this,WGRAD_Module::_defaultRandomSeed2),
121     _pythiaMenu(new PythiaMenu(this,0,"PythiaMenu"))
122 {
123   ostringstream debugDesc;
124   debugDesc  << "\tPrint out debug information ?";
125   _debug.addDescription(debugDesc.str());
126   commands()->append(&_debug);
127 
128   ostringstream usePythiaDesc;
129   usePythiaDesc << "\tUse Pythia for PS and underlying event sim. ?";
130   _usePythia.addDescription(usePythiaDesc.str());
131   _useLesHouchesModule.addDescription("\t Use Les Houches Modules for all interfacing");
132   commands()->append(&_usePythia);
133   commands()->append(&_useLesHouchesModule);
134 
135   //---------------------------------------------------------------------
136   // Cards Menu :
137   //---------------------------------------------------------------------
138   _cardsMenu.initialize("cardsMenu",this);
139   _cardsMenu.initTitle("cardsMenu");
140   commands()->append(&_cardsMenu);
141 
142   _cardsMenu.commands()->append(&_beamType);
143   _cardsMenu.commands()->append(&_W_charge);
144   _cardsMenu.commands()->append(&_widthTreatment);
145   _cardsMenu.commands()->append(&_decayLepton);
146   _cardsMenu.commands()->append(&_QED_radType);
147   _cardsMenu.commands()->append(&_fsrType);
148   _cardsMenu.commands()->append(&_smearingRecombination);
149   _cardsMenu.commands()->append(&_order);
150   _cardsMenu.commands()->append(&_scale);
151   _cardsMenu.commands()->append(&_fsrCollinearCut);
152   _cardsMenu.commands()->append(&_photonEnergyCut);
153   _cardsMenu.commands()->append(&_photonCollinearCut);
154   _cardsMenu.commands()->append(&_pdflibGroup);
155   _cardsMenu.commands()->append(&_pdflibSet);
156   _cardsMenu.commands()->append(&_wMass);
157   _cardsMenu.commands()->append(&_wWidth);
158   _cardsMenu.commands()->append(&_CME);
159   _cardsMenu.commands()->append(&_lepton_ptmin);
160   _cardsMenu.commands()->append(&_lepton_ymax);
161   _cardsMenu.commands()->append(&_neutrino_ptmin);
162   _cardsMenu.commands()->append(&_neutrino_ymax);
163   _cardsMenu.commands()->append(&_photon_ptmin);
164   _cardsMenu.commands()->append(&_photon_ymax);
165   _cardsMenu.commands()->append(&_tmass_low);
166   _cardsMenu.commands()->append(&_tmass_high);
167   _cardsMenu.commands()->append(&_itmaxGrid);  
168   _cardsMenu.commands()->append(&_ncallGrid);
169   _cardsMenu.commands()->append(&_itmaxPrecision);
170   _cardsMenu.commands()->append(&_ncallPrecision);
171   _cardsMenu.commands()->append(&_rndmSeed);
172 
173   //---------------------------------------------------------------------
174   // Pythia Interface Menu :
175   //---------------------------------------------------------------------
176   // No longer the correct way to setup Pythia : 
177   //_pythiaInterfaceMenu.initialize("pythiaInterfaceMenu",this);
178   //_pythiaInterfaceMenu.initTitle("pythiaInterfaceMenu");
179   //commands()->append(&_pythiaInterfaceMenu);
180   //_pythiaInterfaceMenu.commands()->append(&_mstj_41);
181   //_pythiaInterfaceMenu.commands()->append(&_mrpy_1);
182   //_pythiaInterfaceMenu.commands()->append(&_mstp_81);
183   //_pythiaInterfaceMenu.commands()->append(&_mstp_91);
184   //_pythiaInterfaceMenu.commands()->append(&_mstp_111);
185   //_pythiaInterfaceMenu.commands()->append(&_parp_91);
186   //_pythiaInterfaceMenu.commands()->append(&_parp_92);
187   //_pythiaInterfaceMenu.commands()->append(&_parp_93);
188 
189   _pythiaMenu->initialize("PythiaMenu",this);
190   _pythiaMenu->initTitle("Pythia Menu");
191   commands()->append(_pythiaMenu);
192   
193   // Initialize the relevant submenu
194   _randomNumberMenu.initialize("RandomNumberMenu",this);
195   _randomNumberMenu.initTitle("Fake event random number menu");
196   commands()->append(&_randomNumberMenu);
197 
198   // Add the commands to the menu
199   _randomNumberMenu.commands()->append(&_randomSeed1);
200   _randomNumberMenu.commands()->append(&_randomSeed2);
201 
202   ostringstream tmpSstream1;
203   ostringstream tmpSstream2;
204 
205   // Describe them
206   tmpSstream1 << "      \t\t\tSeed #1 for the random number generator"
207               << "\n\t\t\t(default " << _randomSeed1.value() << ").";
208   _randomSeed1.addDescription(tmpSstream1.str());
209   tmpSstream2 << "      \t\t\tSeed #2 for the random number generator"
210               << "\n\t\t\t(default " << _randomSeed2.value() << ").";  
211   _randomSeed2.addDescription(tmpSstream2.str());
212 
213 }
214 
215 //--------------
216 // Destructor --
217 //--------------
218 
219 WGRAD_Module::~WGRAD_Module() 
220 { 
221   delete _pythiaMenu;
222 }
223 
224 //--------------
225 // Operations --
226 //--------------
227 
228 AppResult WGRAD_Module::genBeginJob( )
229 {
230   Pythia* pythia = Pythia::Instance();
231   // Set the random number seeds :
232     CdfRn* rn = CdfRn::Instance();
233     if ( !rn->isReadingFromFile() ) {
234       rn->SetEngineSeeds(_randomSeed1.value(), _randomSeed2.value(),"WGRAD");
235     }
236 
237   // Call WGRAD directly.
238   std::cout << "  WGRAD_Module::genBeginJob : passing cuts and parameters to WGRAD ..." << std::endl;
239 
240   // Set the WGRAD internal parameters from the talk-to parameters :
241   wgrad_set_external_();
242 
243   // Use Pythia or not ?
244   int wgrad_use_pythia[1];
245   if (_usePythia.value())
246     {
247       wgrad_use_pythia[0]=1;
248     }
249   else
250     {
251       wgrad_use_pythia[0]=0;
252     }
253   wgrad_use_pythia_(wgrad_use_pythia);
254 
255   // Set some basic flags and switches :
256   int wgrad_switches[10];
257   wgrad_switches[0] = _beamType.value();
258   wgrad_switches[1] = _W_charge.value();
259   wgrad_switches[2] = _decayLepton.value();
260   wgrad_switches[3] = _QED_radType.value();
261   wgrad_switches[4] = _fsrType.value();
262   wgrad_switches[5] = _smearingRecombination.value();
263   wgrad_switches[6] = _order.value();
264   wgrad_switches[7] = _scale.value();
265   wgrad_switches[8] = _fsrCollinearCut.value();
266   wgrad_switches[9] = _widthTreatment.value();
267   
268   wgrad_set_switches_(wgrad_switches);
269   
270   // Set soft and collinear photon cut-offs :
271   double wgrad_soft_collinear_cuts[2];
272   wgrad_soft_collinear_cuts[0] = _photonEnergyCut.value();
273   wgrad_soft_collinear_cuts[1] = _photonCollinearCut.value();
274   
275   wgrad_set_soft_collinear_cuts_(wgrad_soft_collinear_cuts);
276   
277   // Set choice of structure functions :
278   int wgrad_pdflib[2];
279   wgrad_pdflib[0] = _pdflibGroup.value();
280   wgrad_pdflib[1] = _pdflibSet.value();
281   
282   wgrad_set_pdflib_(wgrad_pdflib);
283   
284   // Set W mass, width and center of mass energy :
285   double wgrad_properties[3];
286   wgrad_properties[0]=_wMass.value();
287   wgrad_properties[1]=_wWidth.value();
288   wgrad_properties[2]=_CME.value();
289   
290   wgrad_set_properties_(wgrad_properties);
291   
292   // Set kinematic cuts :
293   double wgrad_kinematic_cuts[8];
294   wgrad_kinematic_cuts[0] = _lepton_ptmin.value();
295   wgrad_kinematic_cuts[1] = _lepton_ymax.value();
296   wgrad_kinematic_cuts[2] = _neutrino_ptmin.value();
297   wgrad_kinematic_cuts[3] = _neutrino_ymax.value();
298   wgrad_kinematic_cuts[4] = _photon_ptmin.value();
299   wgrad_kinematic_cuts[5] = _photon_ymax.value();
300   wgrad_kinematic_cuts[6] = _tmass_low.value();
301   wgrad_kinematic_cuts[7] = _tmass_high.value();
302 
303   wgrad_set_kinematic_cuts_(wgrad_kinematic_cuts);
304   
305   // Set parameters for grid search and precision calculation :
306   int wgrad_grid[4];
307   wgrad_grid[0] = _itmaxGrid.value();
308   wgrad_grid[1] = _ncallGrid.value();
309   wgrad_grid[2] = _itmaxPrecision.value();
310   wgrad_grid[3] = _ncallPrecision.value();
311   
312   wgrad_set_grid_pars_(wgrad_grid);
313   
314   // Set the random number seed :
315   int wgrad_rndm_seed[1];
316   wgrad_rndm_seed[0] = _rndmSeed.value();
317   wgrad_set_rndm_seed_(wgrad_rndm_seed);
318   
319   // Set the information in the Pythia common blocks for the WGRAD-Pythia interface :
320   // int wgrad_pythia_mstj[1];
321   // wgrad_pythia_mstj[0]=_mstj_41.value();
322   // wgrad_pythia_set_mstj_(wgrad_pythia_mstj);
323   // int wgrad_pythia_mrpy[1];
324   // wgrad_pythia_mrpy[0]=_mrpy_1.value();
325   // wgrad_pythia_set_mrpy_(wgrad_pythia_mrpy);
326   // int wgrad_pythia_mstp[3];
327   // wgrad_pythia_mstp[0]=_mstp_81.value();
328   // wgrad_pythia_mstp[1]=_mstp_91.value();
329   // wgrad_pythia_mstp[2]=_mstp_111.value();
330   // wgrad_pythia_set_mstp_(wgrad_pythia_mstp);
331   // double wgrad_pythia_parp[3];
332   // wgrad_pythia_parp[0]=_parp_91.value();
333   // wgrad_pythia_parp[1]=_parp_92.value();
334   // wgrad_pythia_parp[2]=_parp_93.value();
335   // wgrad_pythia_set_parp_(wgrad_pythia_parp);
336   
337   wgrad_begin_();
338 
339   std::cout << "  WGRAD_Module::genBeginJob : WGRAD initialised. " << std::endl;      
340   // AST, 17.03.2003. Only call PYTHIA directly if LesHouches is not used
341   if(!_useLesHouchesModule.value()) {
342   // DSW, 23.07.2002. Use new Pythia singleton interface.
343     pythia->setPythiaUid("WGRAD");
344 
345   // Call with special value for frame in order to trigger UPINIT :
346     std::string special_frame = "USER";
347     std::string target = "dummy";
348     double winval = 1;
349     std::cout << "  WGRAD_Module::genBeginJob : calling pythia->Pyinit " << std::endl;  
350     pythia->Pyinit(special_frame,target,&winval);
351   }
352 
353   return AppResult::OK;
354 }
355 
356 AppResult WGRAD_Module::genBeginRun( AbsEvent* run )
357 {
358   return AppResult::OK;
359 }
360 
361 AppResult WGRAD_Module::genEndJob()
362 {
363   // In case this gets called without WGRAD being initialised,
364   // set some flags which control termination :
365   wgrad_set_external_();
366   // Call WGRAD subroutine to print out summary information :
367   wgrad_end_();
368 
369   return AppResult::OK;
370 }
371 
372 AppResult WGRAD_Module::genEndRun( AbsEvent* run )
373 {
374   return AppResult::OK;
375 }
376 
377 int WGRAD_Module::callGenerator( AbsEvent* anEvent ) 
378 {
379   // Call WGRAD directly.
380   wgrad_evt_();
381 
382   // Clear up HEPEVT before filling it again :
383   CdfHepevt* hepevt = CdfHepevt::Instance();
384   hepevt->clear();
385   hepevt->clearCommon();
386 
387   // Extract the event weight from WGRAD :
388   double wgrad_event_weight[1];
389   wgrad_get_event_weight_(wgrad_event_weight);
390 
391   // AST, 17.03.2003. If using Les Houches Module let it do the rest       
392   if(_useLesHouchesModule.value()) return 0;
393 
394   // Pythia or not
395   if (!_usePythia.value()) 
396     {
397       // Extract the charged lepton, neutrino and photon four momenta from WGRAD :
398       int wgrad_particle_codes[3];
399       wgrad_get_particle_codes_(wgrad_particle_codes);
400       double wgrad_charged_lepton[4];
401       wgrad_get_charged_lepton_(wgrad_charged_lepton);
402       double wgrad_neutrino[4];
403       wgrad_get_neutrino_(wgrad_neutrino);
404       double wgrad_photon[4];
405       wgrad_get_photon_(wgrad_photon);
406 
407       if (_debug.value())
408         {
409           std::cout << " WGRAD_Module::callGenerator : Charged Lepton code, px, py, pz, E = " << wgrad_particle_codes[0] << " , " 
410                << wgrad_charged_lepton[0] << " , " << wgrad_charged_lepton[1] << " , " 
411                << wgrad_charged_lepton[2] << " , " << wgrad_charged_lepton[3] << std::endl;
412           std::cout << " WGRAD_Module::callGenerator : Neutrino code, px, py, pz, E       = " << wgrad_particle_codes[1] << " , " 
413                << wgrad_neutrino[0] << " , " << wgrad_neutrino[1] << " , " 
414                << wgrad_neutrino[2] << " , " << wgrad_neutrino[3] << std::endl;
415           std::cout << " WGRAD_Module::callGenerator : Photon code, px, py, pz, E         = " << wgrad_particle_codes[2] << " , " 
416                << wgrad_photon[0] << " , " << wgrad_photon[1] << " , " 
417                << wgrad_photon[2] << " , " << wgrad_photon[3] << std::endl;
418         }
419 
420 
421       // Determine whether 2 or 3 body final state :
422       int nPart = 3;
423       if (wgrad_photon[3] == 0.0) nPart = 2;
424       if (_debug.value()) std::cout << " WGRAD_Module::callGenerator : nPart = " << nPart << std::endl;
425 
426       CdfHepevt* hepevt = CdfHepevt::Instance();
427 
428       for (int iPart = 0; iPart < nPart; iPart++)
429         {         
430           hepevt->HepevtPtr()->IDHEP[iPart]     = wgrad_particle_codes[iPart];
431                   
432           hepevt->HepevtPtr()->ISTHEP[iPart]    = 1;   // status code for final state particles
433           hepevt->HepevtPtr()->JMOHEP[iPart][0] = 0;   // mother position 
434           hepevt->HepevtPtr()->JMOHEP[iPart][1] = 0;   // second mother position
435                   
436           hepevt->HepevtPtr()->JDAHEP[iPart][0] = 0;   // first daughter
437           hepevt->HepevtPtr()->JDAHEP[iPart][1] = 0;   // last daughter
438                   
439           if (iPart==0)
440             {
441               hepevt->HepevtPtr()->PHEP[iPart][0]   = wgrad_charged_lepton[0];
442               hepevt->HepevtPtr()->PHEP[iPart][1]   = wgrad_charged_lepton[1];
443               hepevt->HepevtPtr()->PHEP[iPart][2]   = wgrad_charged_lepton[2];
444               hepevt->HepevtPtr()->PHEP[iPart][3]   = wgrad_charged_lepton[3];
445             }
446           else if (iPart==1)
447             {
448               hepevt->HepevtPtr()->PHEP[iPart][0]   = wgrad_neutrino[0];
449               hepevt->HepevtPtr()->PHEP[iPart][1]   = wgrad_neutrino[1];
450               hepevt->HepevtPtr()->PHEP[iPart][2]   = wgrad_neutrino[2];
451               hepevt->HepevtPtr()->PHEP[iPart][3]   = wgrad_neutrino[3];
452             }
453           else if (iPart==2)
454             {
455               hepevt->HepevtPtr()->PHEP[iPart][0]   = wgrad_photon[0];
456               hepevt->HepevtPtr()->PHEP[iPart][1]   = wgrad_photon[1];
457               hepevt->HepevtPtr()->PHEP[iPart][2]   = wgrad_photon[2];
458               hepevt->HepevtPtr()->PHEP[iPart][3]   = wgrad_photon[3];
459             }
460 
461           hepevt->HepevtPtr()->VHEP[iPart][0]   = 0.;  // production vertex x-position
462           hepevt->HepevtPtr()->VHEP[iPart][1]   = 0.;  // production vertex y-position
463           hepevt->HepevtPtr()->VHEP[iPart][2]   = 0.;  // production vertex z-position
464           hepevt->HepevtPtr()->VHEP[iPart][3]   = 0.;  // production vertex time
465         }             
466       hepevt->HepevtPtr()->NHEP         = nPart;   // number of particles
467     }
468   else
469     {
470       // DSW, 23.07.2002. Use new Pythia singleton interface.
471       Pythia* pythia = Pythia::Instance();
472       // std::cout << "  WGRAD_Module::callGenerator : calling pythia->Pyevnt()" << std::endl;
473       // Let Pythia do it's stuff :
474       pythia->Pyevnt();
475       // Update the event counter :
476       pythia->addEvent();
477       int mevlist = pythia->EventListLevel();
478       if (pythia->events() >= _pythiaMenu->firstlistevent() &&
479           pythia->events() <= _pythiaMenu->lastlistevent())
480         {
481           std::cout << "PYEVNT event number " << pythia->events() << std::endl;
482           pythia->Pylist(&mevlist);
483         }
484       // Tell lunhep to convert from PYJETS to HEPEVT
485       int mconv = 1;
486       // Fill HEPEVT :
487       pythia->Lunhep(&mconv);
488     }
489   
490   return 0;
491 }
492 
493 void WGRAD_Module::fillHepevt() {
494   // this is done somewhere else, should be moved here
495 }
496
[ source navigation ] [ diff markup ] [ identifier search ] [ freetext search ] [ file search ]
This page was automatically generated by the LXR engine.
The LXR team
Valid HTML 4.01!
Send problems or questions to cdfcode@fnal.gov