BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Generator/Bhwide/Bhwide-00-00-11/src/Bhwide.cxx Source File

00001 //*****************************************************************************
00002 //
00003 // Bhwide.cxx
00004 //
00005 // Algorithm runs large angle Bhabha event generator Bhwide 
00006 // and stores output to transient store
00007 //
00008 // Aug 2007 A. Zhemchugov, initial version written for BES3
00009 //*****************************************************************************
00010 
00011 // Header for this module:-
00012 #include "Bhwide/Bhwide.h"
00013 #include "Bhwide/BhwideRandom.h"
00014 
00015 // Framework Related Headers:-
00016 #include "HepMC/GenEvent.h"
00017 #include "HepMC/GenVertex.h"
00018 #include "HepMC/GenParticle.h"
00019 #include "CLHEP/Vector/LorentzVector.h"
00020 
00021 #include "GaudiKernel/MsgStream.h"
00022 #include "GaudiKernel/ISvcLocator.h"
00023 #include "GaudiKernel/AlgFactory.h"
00024 #include "GaudiKernel/DataSvc.h"
00025 #include "GaudiKernel/SmartDataPtr.h"
00026 
00027 #include "GeneratorObject/McGenEvent.h"
00028 #include "BesKernel/IBesRndmGenSvc.h"
00029 
00030 #include "cfortran/cfortran.h"
00031 
00032 #include <stdlib.h>
00033 
00034 using namespace std;
00035 //      COMMON / momset / p1(4),q1(4),p2(4),q2(4),phot(100,4),nphot
00036 //      * COMMON / wgtall / wtmod,wtcru1,wtcru2,wtset(300)
00037 
00038 typedef struct { 
00039  double p1[4]; 
00040  double q1[4]; 
00041  double p2[4]; 
00042  double q2[4]; 
00043  double phot[4][100]; 
00044  int nphot; } MOMSET_DEF;
00045 
00046 #define MOMSET COMMON_BLOCK(MOMSET_DEF, momset)
00047 COMMON_BLOCK_DEF(MOMSET_DEF,MOMSET);
00048 //MOMSET_DEF MOMSET; 
00049 
00050 //-------------------------
00051 
00052 PROTOCCALLSFSUB1(GLIMIT,glimit,INT)
00053 #define GLIMIT(LENMX) CCALLSFSUB1(GLIMIT,glimit,INT,LENMX)
00054                                                                                             
00055 PROTOCCALLSFSUB1(DUMPS,dumps,INT)
00056 #define DUMPS(NOUT) CCALLSFSUB1(DUMPS,dumps,INT,NOUT)
00057                                                                                             
00058 PROTOCCALLSFSUB3(BHWIDE,bhwide,INT,DOUBLEV,INTV)
00059 #define BHWIDE(MODE,XPAR,NPAR) CCALLSFSUB3(BHWIDE,bhwide,INT,DOUBLEV,INTV,MODE,XPAR,NPAR)
00060 
00061 PROTOCCALLSFSUB3(MARINI, marini, INT, INT, INT)
00062 #define MARINI(IJKLIN, NTOTIN, NTOT2N) CCALLSFSUB3(MARINI, marini, INT, INT, INT, IJKLIN, NTOTIN, NTOT2N)
00063 
00064 extern "C" {                                                                                           
00065   extern double ranmarr_();
00066   extern void marran_(float* rvec, int* nma);
00067   extern void ecuran_(double* rvec, int* nma);
00068   extern void carran_(double* rvec, int* nma);
00069 }
00070 
00071 double ranmarr_()
00072 {
00073   return BhwideRandom::random();
00074 }
00075 
00076 void marran_(float* rvec, int* nma)
00077 {
00078   int nmax = *nma;
00079   assert(nmax<100);
00080   double rvecd[100];
00081   BhwideRandom::FlatArray(rvecd, nmax);
00082   for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00083 }
00084 
00085 void carran_(double* rvec, int* nma)
00086 {
00087   int nmax = *nma;
00088   assert(nmax<100);
00089   double rvecd[100];
00090   BhwideRandom::FlatArray(rvecd, nmax);
00091   for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00092 }
00093 
00094 void ecuran_(double* rvec, int* nma)
00095 {
00096   int nmax = *nma;
00097   assert(nmax<100);
00098   double rvecd[100];
00099   BhwideRandom::FlatArray(rvecd, nmax);
00100   for(int i=0; i<nmax; i++) rvec[i]=rvecd[i];
00101 }
00102 
00103 
00104 
00105 
00106 //float ranmar_(){
00107 //    return BhwideRandom::random();                                                              
00108 //}                                                                                                      
00109                                                                                                        
00110 Bhwide::Bhwide(const string& name, ISvcLocator* pSvcLocator):Algorithm( name, pSvcLocator )
00111 {
00112   declareProperty("CMEnergy", m_cmEnergy = 3.097); // 2*Ebeam [GeV]
00113   declareProperty("EleMinThetaAngle", m_ThMine = 22); // [deg]
00114   declareProperty("EleMaxThetaAngle", m_ThMaxe = 158); // [deg]
00115   declareProperty("PosMinThetaAngle", m_ThMinp = 22); // [deg]
00116   declareProperty("PosMaxThetaAngle", m_ThMaxp = 158); // [deg]
00117   declareProperty("EleMinEnergy", m_EnMine = 0.01); // [GeV]
00118   declareProperty("PosMinEnergy", m_EnMinp = 0.01); // [GeV]
00119   declareProperty("Acollinearity", m_Acolli = 10);   // [deg]
00120   declareProperty("SoftPhotonCut", m_infraredCut   = 1E-5); // dimensionless, Ephoton > m_infraredCut*sqrt(s)/2
00121   declareProperty("VacuumPolarization", m_keyPia   = 3); // 0 - OFF, 1 - ON, Burkhardt et.al. 1989, as in BHLUMI 2.0x, 2 - ON, S. Eidelman, F. Jegerlehner, Z.Phys.C(1995), 3 - ON, Burkhardt and Pietrzyk 1995 (Moriond).
00122   declareProperty("KeyMod", m_keyMod   = 2); // type of MODEL subprogram and QED matrix element for hard bremsstrahlung: 1 - obtained by the authors (helicity amplitudes), 2 - from CALKUL, Nucl. Phys. B206 (1982) 61. Checked to be in a very good agreement!
00123   declareProperty("KeyLib", m_keyLib   = 2); // option for ElectroWeak Corrections Library: 1 - ElectroWeak Corr. from BABAMC (obsolete), 2 - ElectroWeak Corr. from ALIBABA, RECOMMENDED 
00124   declareProperty("EWC", m_keyEwc   = 1); // switching ON/OFF weak corrections: 0 - only QED corrections included (here both KeyLib =1,2 should be equivalent), 1 - all ElectroWeak Corrections included 
00125   m_initSeed.clear();
00126 //  m_initSeed.push_back(54217137); m_initSeed.push_back(0); m_initSeed.push_back(0);
00127 //  declareProperty("InitializedSeed", m_initSeed);
00128 }
00129 
00130 StatusCode Bhwide::initialize(){
00131 
00132   MsgStream log(messageService(), name());
00133   
00134   log << MSG::INFO << "Bhwide initialize" << endreq;
00135 
00136   static const bool CREATEIFNOTTHERE(true);                                                            
00137   StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);                 
00138   if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)                                                 
00139   {                                                                                                    
00140     log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;                      
00141     return RndmStatus;                                                                                 
00142   }                                                                                                    
00143   CLHEP::HepRandomEngine* engine  = p_BesRndmGenSvc->GetEngine("Bhwide");                            
00144   engine->showStatus();
00145   BhwideRandom::setRandomEngine(engine);                                                             
00146 
00147   GLIMIT(50000);
00148 
00152   double WTMAX  =  3.0; //    ! Maximum Weight for rejection
00153   double AMAZ   =  91.1882; // ! Z mass
00154   double GAMMZ  =  2.4952; //  ! Z width (may be recalculated by EW library)
00155   double SINW2  =  0.22225; // ! sin^2(theta_W) (may be recalculated by EW library)   
00156   double AMTOP  = 174.3; //    ! top quark mass
00157   double AMHIG  = 115.0; //    ! Higgs mass      
00159   int KeyWgt =   0;     // ! unweighted (WT=1) events, for detector simulation
00161   int    KeyRnd =   1;  // ! RANMAR random numbers
00162   int    KeyCha =   0;  // ! Channel choice: all/s-only/t-only: =0/1/2
00163   int    KeyZof =   0;  // ! Z-contribution ON/OFF: =0/1
00164   int    KeyOpt = 1000*KeyZof +100*KeyCha +10*KeyWgt + KeyRnd;
00165   // !#      KeyEWC = 0   ! QED corrections only
00166   int    KeyEWC =   m_keyEwc;  // ! Total O(alpha) ElectroWeak corr. included
00167   // !#      KeyLib = 1   ! ElectroWeak corrections from BABAMC (obsolete!)
00168   int   KeyLib =  m_keyLib;  // ! ElectroWeak corrections from ALIBABA
00169   // !#      KeyMod = 1   ! Hard bremsstr. matrix element from MODEL1
00170   int    KeyMod =   m_keyMod;  // ! Hard bremsstr. matrix alement from MODEL2
00171   int    KeyPia =   m_keyPia;  // ! Vacuum polarization option (0/1/2/3)
00172   int    KeyRad = 1000*KeyEWC + 100*KeyLib + 10*KeyMod + KeyPia;
00173   
00175   npar[0]=   KeyOpt;
00176   npar[1]=   KeyRad;
00177 
00178   xpar[0] = m_cmEnergy;  
00179   xpar[1] = m_ThMinp;
00180   xpar[2] = m_ThMaxp;
00181   xpar[3] = m_ThMine;
00182   xpar[4] = m_ThMaxe;
00183   xpar[5] = m_EnMinp;
00184   xpar[6] = m_EnMine;
00185   xpar[7] = m_Acolli;
00186   xpar[8] = m_infraredCut;  //        ! Infrared cut on photon energy
00187   xpar[9]  = WTMAX;
00188   xpar[10] = AMAZ;
00189   xpar[11] = GAMMZ;
00190   xpar[12] = SINW2;
00191   xpar[13] = AMTOP;
00192   xpar[14] = AMHIG;
00193   
00194   //MARINI(m_initSeed[0], m_initSeed[1], m_initSeed[2]);
00195   
00196   BHWIDE(-1,xpar,npar);
00197   
00198   return StatusCode::SUCCESS;
00199 }
00200 
00201 StatusCode Bhwide::execute() 
00202 { 
00203   MsgStream log(messageService(), name());
00204   log << MSG::INFO << "Bhwide executing" << endreq;
00205 
00206 
00207   BHWIDE( 0,xpar,npar);
00208 
00209   if( log.level() < MSG::INFO )
00210   {
00211     DUMPS(6);
00212    // dump output to file
00213    //  DUMPS(16);
00214   }
00215 
00216   int npart = 0;
00217   
00218   // Fill event information
00219   GenEvent* evt = new GenEvent(1,1);
00220 
00221   GenVertex* prod_vtx = new GenVertex();
00222 //  prod_vtx->add_particle_out( p );
00223   evt->add_vertex( prod_vtx );
00224                           
00225   // incoming beam e+
00226   GenParticle* p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p1[0],  MOMSET.p1[1], 
00227                                                       MOMSET.p1[2],  MOMSET.p1[3]),
00228                                                       -11, 3); 
00229   p->suggest_barcode( ++npart );
00230   prod_vtx->add_particle_in(p);
00231 
00232   // incoming beam e-
00233   p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q1[0],  MOMSET.q1[1], MOMSET.q1[2],  MOMSET.q1[3]),
00234                        11, 3); 
00235   p->suggest_barcode( ++npart );
00236   prod_vtx->add_particle_in(p);
00237   
00238   // scattered e+
00239   p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.p2[0],  MOMSET.p2[1], 
00240                                          MOMSET.p2[2],  MOMSET.p2[3]),
00241                                          -11, 1); 
00242   p->suggest_barcode( ++npart );
00243   prod_vtx->add_particle_out(p);
00244 
00245   // scattered e-
00246   p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.q2[0],  MOMSET.q2[1], MOMSET.q2[2],  MOMSET.q2[3]),
00247                        11, 1); 
00248   p->suggest_barcode( ++npart );
00249   prod_vtx->add_particle_out(p);
00250 
00251   int iphot=0;
00252   for (iphot=0; iphot<MOMSET.nphot; iphot++)
00253   {
00254     // gamma
00255     p = new GenParticle( CLHEP::HepLorentzVector( MOMSET.phot[0][iphot],  MOMSET.phot[1][iphot], 
00256                                            MOMSET.phot[2][iphot],  MOMSET.phot[3][iphot]),
00257                                            22, 1); 
00258     p->suggest_barcode( ++npart );
00259     prod_vtx->add_particle_out(p);
00260   }
00261   
00262   if( log.level() < MSG::INFO )
00263   {
00264     evt->print();  
00265   }
00266 
00267   // Check if the McCollection already exists
00268   SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00269   if (anMcCol!=0) 
00270   {
00271     // Add event to existing collection
00272     MsgStream log(messageService(), name());
00273     log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00274     McGenEvent* mcEvent = new McGenEvent(evt);
00275     anMcCol->push_back(mcEvent);
00276   }
00277   else 
00278   {
00279     // Create Collection and add  to the transient store
00280     McGenEventCol *mcColl = new McGenEventCol;
00281     McGenEvent* mcEvent = new McGenEvent(evt);
00282     mcColl->push_back(mcEvent);
00283     StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00284     if (sc != StatusCode::SUCCESS) 
00285     {
00286       log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00287       delete mcColl;
00288       delete evt;
00289       delete mcEvent;
00290       return StatusCode::FAILURE;
00291     }
00292     else 
00293     {
00294       log << MSG::INFO << "McGenEventCol created and " << npart <<" particles stored in McGenEvent" << endreq;
00295     }
00296   }
00297 
00298   return StatusCode::SUCCESS; 
00299 }
00300 
00301 StatusCode Bhwide::finalize() 
00302 {
00303   MsgStream log(messageService(), name());
00304 
00305   BHWIDE( 2,xpar,npar);
00306 
00307   log << MSG::INFO << "Bhwide finalized" << endreq;
00308 
00309   return StatusCode::SUCCESS;
00310 }