Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

Phokhara Class Reference

#include <Phokhara.h>

List of all members.

Public Member Functions

StatusCode execute ()
StatusCode execute ()
StatusCode finalize ()
StatusCode finalize ()
StatusCode initialize ()
StatusCode initialize ()
 Phokhara (const string &name, ISvcLocator *pSvcLocator)
 Phokhara (const string &name, ISvcLocator *pSvcLocator)
StatusCode storeParticles ()
StatusCode storeParticles ()

Private Attributes

double Ar [14]
double Ar_r [14]
double cos1max
double cos1min
double cos2max
double cos2min
double cos3max
double cos3min
double dsigm
double dsigm1
double dsigm2
int ievent
double m_E
int m_f0_model
int m_FF_Pion
int m_fsr
int m_fsrnlo
double m_gmin
long int m_initSeed
int m_ivac
int m_nlo
int m_nm
double m_phot1cut
double m_phot2cut
double m_pi1cut
double m_pi2cut
int m_pion
double m_q2_max_c
double m_q2_min_c
double m_q2min
int m_tagged
double m_w
double qqmax
double qqmin
double sigma
double sigma1
double sigma2

Constructor & Destructor Documentation

Phokhara::Phokhara (  const string &  name,
ISvcLocator *  pSvcLocator
) 
  

00039                                                               :Algorithm( name, pSvcLocator )
00040 {
00041   m_initSeed = (0);
00042 
00043   declareProperty( "InitialSeed", m_initSeed = 0);
00044 
00045   declareProperty( "InitialEvents",    m_nm = 50000 );       // # of events to determine the maximum   
00046   declareProperty( "NLO",              m_nlo = 1 );          // Born(0), NLO(1)
00047   declareProperty( "SoftPhotonCutoff", m_w = 0.0001 );       // soft photon cutoff                
00048   declareProperty( "Channel",          m_pion = 0 );         // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
00049                                                              // 2pi+2pi-(3),ppbar(4),nnbar(5),
00050                                                              // K+K-(6),K0K0bar(7),pi+pi-pi0(8), 
00051                                                              // Lamb Lambbar->pi-pi+ppbar(9) 
00052   declareProperty( "FSR",                m_fsr = 1 );        // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
00053   declareProperty( "FSRNLO",             m_fsrnlo = 1 );     // yes(1), no(0)
00054   declareProperty( "VacuumPolarization", m_ivac = 0 );       // yes(1), no(0)
00055   declareProperty( "TaggedPhotons",      m_tagged = 0 );     // tagged photons(0), untagged photons(1)
00056   declareProperty( "PionFormfactor",     m_FF_Pion = 0 );    // KS Pionformfactor(0), GS Pionformfactor(1)
00057   declareProperty( "F0model",            m_f0_model = 0 );   // f0+f0(600): KK model(0), no structure(1), 
00058                                                              // no f0+f0(600)(2), f0 KLOE(3)
00059   declareProperty( "Ecm",             m_E = 3.097 );         // CMS-energy (GeV)                       
00060   declareProperty( "CutTotalInvMass", m_q2min = 0.0 );       // minimal  hadrons(muons)-gamma-inv mass squared (GeV^2)
00061   declareProperty( "CutHadInvMass",   m_q2_min_c = 0.0447 ); // minimal inv. mass squared of the hadrons(muons)(GeV^2)
00062   declareProperty( "MaxHadInvMass",   m_q2_max_c = 1.06 );   // maximal inv. mass squared of the hadrons(muons)(GeV^2)
00063   declareProperty( "MinPhotonEnergy", m_gmin = 0.05 );       // minimal photon energy/missing energy   (GeV)          
00064   declareProperty( "MinPhotonAngle",  m_phot1cut = 0.0 );    // minimal photon angle/missing momentum angle (grad)
00065   declareProperty( "MaxPhotonAngle",  m_phot2cut = 180.0 );  // maximal photon angle/missing momentum angle (grad)
00066   declareProperty( "MinHadronsAngle", m_pi1cut = 0.0 );      // minimal hadrons(muons) angle (grad)
00067   declareProperty( "MaxHadronsAngle", m_pi2cut = 180.0 );    // maximal hadrons(muons) angle (grad)
00068 
00069   ievent = 0;
00070 }

Phokhara::Phokhara (  const string &  name,
ISvcLocator *  pSvcLocator
) 
 

Member Function Documentation

StatusCode Phokhara::execute (   ) 
 

StatusCode Phokhara::execute (   ) 
  

00465 { 
00466   MsgStream log(messageService(), name());
00467 
00468   cout << "Execute " << ievent << endl;
00469 
00470   int ntrials = 0;
00471   int tr_old[2];
00472   tr_old[0] = (int)MAXIMA.tr[0];
00473   tr_old[1] = (int)MAXIMA.tr[1];
00474 
00475   while( ntrials < 1000)
00476     {
00477       ievent++;
00478       RANLXDF(Ar_r,1);
00479       Ar[1] = Ar_r[0];
00480       
00481       if (Ar[1] <= (MAXIMA.Mmax[0]/(MAXIMA.Mmax[0]+MAXIMA.Mmax[1]))) { 
00482         MAXIMA.count[0] = MAXIMA.count[0]+1.0;
00483         GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
00484       } 
00485       else {
00486         MAXIMA.count[1] = MAXIMA.count[1]+1.0;
00487         GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
00488       }
00489       
00490       if( ((int)MAXIMA.tr[0]+(int)MAXIMA.tr[1]) > (tr_old[0]+tr_old[1]) ) // event accepted after cuts
00491         {
00492           storeParticles();
00493           return StatusCode::SUCCESS; 
00494         }
00495       ntrials ++;
00496     }
00497  
00498   log << MSG::WARNING << "Could not satisfy cuts after " << ntrials << "trials. Continue ..." << endreq;
00499   
00500   return StatusCode::SUCCESS; 
00501 }

StatusCode Phokhara::finalize (   ) 
 

StatusCode Phokhara::finalize (   ) 
  

00727 {
00728   MsgStream log(messageService(), name());
00729   // =================================================
00730   if(FLAGS.pion == 9)
00731     MAXIMA.Mmax[1] = MAXIMA.Mmax[1] * (1.0 + LAMBDA_PAR.alpha_lamb)*(1.0 + LAMBDA_PAR.alpha_lamb) 
00732       * LAMBDA_PAR.ratio_lamb*LAMBDA_PAR.ratio_lamb;
00733   
00734   // --- value of the cross section ------------------
00735   if (FLAGS.nlo == 0) { 
00736     sigma = MAXIMA.Mmax[0]/MAXIMA.count[0]*MAXIMA.tr[0];
00737     dsigm = MAXIMA.Mmax[0]*sqrt((MAXIMA.tr[0]/MAXIMA.count[0]-(MAXIMA.tr[0]/MAXIMA.count[0])*(MAXIMA.tr[0]/MAXIMA.count[0]))/MAXIMA.count[0]);
00738   } else {
00739     sigma1 = MAXIMA.Mmax[0]/MAXIMA.count[0]*MAXIMA.tr[0];
00740     sigma2 = MAXIMA.Mmax[1]/MAXIMA.count[1]*MAXIMA.tr[1];
00741     dsigm1 = MAXIMA.Mmax[0]*sqrt((MAXIMA.tr[0]/MAXIMA.count[0]-(MAXIMA.tr[0]/MAXIMA.count[0])*(MAXIMA.tr[0]/MAXIMA.count[0]))/MAXIMA.count[0]);
00742     dsigm2 = MAXIMA.Mmax[1]*sqrt((MAXIMA.tr[1]/MAXIMA.count[1]-(MAXIMA.tr[1]/MAXIMA.count[1])*(MAXIMA.tr[1]/MAXIMA.count[1]))/MAXIMA.count[1]);
00743     
00744     sigma = sigma1+sigma2;
00745     dsigm = sqrt(dsigm1*dsigm1 + dsigm2*dsigm2);
00746   }
00747   
00748   // --- output --------------------------------------
00749   cout << "-------------------------------------------------------------" << endl;
00750   cout << "     PHOKHARA 6.0 Final Statistcs                            " << endl;
00751   cout << "-------------------------------------------------------------"  << endl;
00752   cout <<  int(MAXIMA.tr[0]+MAXIMA.tr[1]) << " total events accepted of " << endl;
00753   cout <<  int(ievent)            <<  " total events generated" << endl;
00754   cout <<  int(MAXIMA.tr[0])    <<  " one photon events accepted of " << endl;
00755   cout <<  int(MAXIMA.count[0]) <<  " events generated" << endl;
00756   cout <<  int(MAXIMA.tr[1])    <<  " two photon events accepted of " << endl;
00757   cout <<  int(MAXIMA.count[1]) <<  " events generated" << endl;
00758   cout << endl;
00759   if (FLAGS.nlo != 0) cout <<  "sigma1(nbarn) = " << sigma1 << " +- " << dsigm1 << endl;
00760   if (FLAGS.nlo != 0) cout <<  "sigma2(nbarn) = " << sigma2 << " +- " << dsigm2 << endl;
00761   cout <<  "sigma (nbarn) = " << sigma << " +- " <<dsigm << endl;
00762   cout << endl;
00763   cout <<  "maximum1 = " << MAXIMA.gross[0] << "  minimum1 = " << MAXIMA.klein[0] << endl;
00764   cout <<  "Mmax1    = " << MAXIMA.Mmax[0]  << endl;        
00765   cout <<  "maximum2 = " << MAXIMA.gross[1] << "  minimum2 = " << MAXIMA.klein[1] << endl;
00766   cout <<  "Mmax2    = " << MAXIMA.Mmax[1]  << endl;         
00767   cout <<  "-------------------------------------------------------------" << endl;
00768     
00769   log << MSG::INFO << "Phokhara finalized" << endreq;
00770 
00771   return StatusCode::SUCCESS;
00772 }

StatusCode Phokhara::initialize (   ) 
 

StatusCode Phokhara::initialize (   ) 
  

00072                                {
00073 
00074   MsgStream log(messageService(), name());
00075   
00076   log << MSG::INFO  << "Phokhara initialize" << endreq;
00077 
00078   int i,s_seed[105];
00079   long int k,j;
00080 
00081   if( m_initSeed == 0) // if seed is not set in the jobptions file, set it using BesRndmGenSvc
00082     {
00083       IBesRndmGenSvc* p_BesRndmGenSvc;
00084       static const bool CREATEIFNOTTHERE(true);
00085       StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
00086       if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
00087         {
00088           log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
00089           return RndmStatus;
00090         }
00091       
00092       // Save the random number seeds in the event
00093       HepRandomEngine* engine  = p_BesRndmGenSvc->GetEngine("PHOKHARA");
00094       const long*  s       = engine->getSeeds();
00095       m_initSeed = s[0];
00096     }
00097 
00098   log << MSG::INFO << "Set initial seed " << m_initSeed << endreq;
00099 
00100   // --- initialize the seed ------
00101   /*  ifstream seeds("seed.dat");
00102   if( seeds.is_open() )
00103     {
00104       int ii=0;
00105       while( ! seeds.eof() )
00106         seeds >> s_seed[ii++];
00107     }
00108   else
00109     cerr << "Cannot open file seed.dat for reading" << endl;
00110   */
00111 
00112   RLXDINIT(1, m_initSeed);
00113 
00114   //RLXDRESETF(s_seed);
00115   //rlxd_reset(s_seed);
00116   
00117   // --- input parameters ----------------------------
00118   MAXIMA.iprint = 0;
00119 
00120   FLAGS.nlo = m_nlo;    
00121   FLAGS.pion = m_pion;   
00122   FLAGS.fsr = m_fsr;    
00123   FLAGS.fsrnlo = m_fsrnlo; 
00124   FLAGS.ivac = m_ivac;    
00125   FLAGS.tagged = m_tagged;  
00126   FLAGS.FF_pion  = m_FF_Pion; 
00127   FLAGS.f0_model = m_f0_model;
00128   
00129   CTES.Sp = m_E*m_E;        ;
00130   
00131   CUTS.w = m_w;
00132   CUTS.q2min = m_q2min;  
00133   CUTS.q2_min_c = m_q2_min_c; 
00134   CUTS.q2_max_c = m_q2_max_c;
00135   CUTS.gmin = m_gmin;    
00136   CUTS.phot1cut = m_phot1cut;
00137   CUTS.phot2cut = m_phot2cut;
00138   CUTS.pi1cut = m_pi1cut;  
00139   CUTS.pi2cut = m_pi2cut;  
00140 
00141   INPUT();
00142   // --- print run data ------------------------------
00143   cout << "-------------------------------------------------------------" << endl;
00144   if (FLAGS.pion == 0) 
00145     cout << "     PHOKHARA 6.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
00146   else if (FLAGS.pion == 1) 
00147     cout << "     PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- gamma" << endl;    
00148   else if (FLAGS.pion == 2) 
00149     cout << "   PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;        
00150   else if (FLAGS.pion == 3) 
00151     cout << "   PHOKHARA 6.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl; 
00152   else if (FLAGS.pion == 4) 
00153     cout << "   PHOKHARA 6.0: e^+ e^- -> p pbar gamma" << endl; 
00154   else if (FLAGS.pion == 5) 
00155     cout << "   PHOKHARA 6.0: e^+ e^- -> n nbar gamma" << endl;                
00156   else if (FLAGS.pion == 6) 
00157     cout << "   PHOKHARA 6.0: e^+ e^- -> K^+ K^- gamma" << endl;                
00158   else if (FLAGS.pion == 7) 
00159     cout << "   PHOKHARA 6.0: e^+ e^- -> K_0 K_0bar gamma" << endl;                
00160   else if (FLAGS.pion == 8) 
00161     cout << "   PHOKHARA 6.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;                
00162   else if (FLAGS.pion == 9) {
00163     cout <<  "PHOKHARA 6.0 : e^+ e^- ->" << endl;
00164     cout << "  Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
00165   } else 
00166     cout <<  "     PHOKHARA 6.0: not yet implemented" << endl;
00167   
00168   // --------------------------------
00169   cout <<  "--------------------------------------------------------------" << endl;
00170   printf(" %s %f %s\n","cms total energy                       = ",sqrt(CTES.Sp)," GeV  "); 
00171   if (FLAGS.tagged == 0) { 
00172     if((CUTS.gmin/2.0/CTES.ebeam) < 0.0098){
00173       cerr << " minimal missing energy set to small" << endl;
00174       return 0;
00175     }
00176     printf(" %s %f %s\n","minimal tagged photon energy           = ",CUTS.gmin," GeV  ");       
00177     printf(" %s %f,%f\n","angular cuts on tagged photon          = ",CUTS.phot1cut,CUTS.phot2cut);
00178   }
00179   else {
00180     if((CUTS.gmin/2.0/CTES.ebeam) < 0.0098){
00181       cerr << " minimal missing energy set to small" << endl;
00182       return 0;
00183     }
00184     printf(" %s %f %s\n","minimal missing energy                 = ",CUTS.gmin," GeV  ");
00185     printf(" %s %f,%f\n","angular cuts on missing momentum       = ",CUTS.phot1cut, CUTS.phot2cut);
00186   }
00187 
00188  
00189   // --------------------------------
00190   if (FLAGS.pion == 0) 
00191     printf(" %s %f,%f\n","angular cuts on muons                  = ",CUTS.pi1cut,CUTS.pi2cut);
00192   else if (FLAGS.pion == 4)  
00193     printf(" %s %f,%f\n","angular cuts on protons                = ",CUTS.pi1cut,CUTS.pi2cut);
00194   else if (FLAGS.pion == 5) 
00195     printf(" %s %f,%f\n","angular cuts on neutrons               = ", CUTS.pi1cut,CUTS.pi2cut);
00196   else if ((FLAGS.pion == 6)||(FLAGS.pion == 7)) 
00197     printf(" %s %f,%f\n","angular cuts on kaons                  = ", CUTS.pi1cut,CUTS.pi2cut);
00198   else if (FLAGS.pion == 9) 
00199     printf(" %s %f,%f\n","angular cuts on pions and protons      = ", CUTS.pi1cut,CUTS.pi2cut);
00200   else
00201     printf(" %s %f,%f\n","angular cuts on pions                  = ", CUTS.pi1cut,CUTS.pi2cut);
00202     
00203   if (FLAGS.tagged == 0) {
00204     if (FLAGS.pion == 0) 
00205       printf(" %s %f %s\n","min. muons-tagged photon inv.mass^2    = ", CUTS.q2min," GeV^2");
00206     else if (FLAGS.pion == 4) 
00207       printf(" %s %f %s\n","min. protons-tagged photon inv.mass^2  = ", CUTS.q2min," GeV^2" );
00208     else if (FLAGS.pion == 5) 
00209       printf(" %s %f %s\n","min. neutrons-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
00210     else if ((FLAGS.pion == 6)||(FLAGS.pion == 7)) 
00211       printf(" %s %f %s\n","min. kaons-tagged photon inv.mass^2    = ", CUTS.q2min," GeV^2" );
00212     else if (FLAGS.pion == 9) 
00213       printf(" %s %f %s\n"," min. lambdas-tagged photon inv.mass^2 = ", CUTS.q2min," GeV^2" );
00214     else
00215       printf(" %s %f %s\n","min. pions-tagged photon inv.mass^2    = ", CUTS.q2min," GeV^2" );
00216   }
00217         
00218 // --- set cuts ------------------------------------
00219       if (FLAGS.tagged == 0) { 
00220         cos1min = cos(CUTS.phot2cut*CTES.pi/180.0);     // photon1 angle cuts in the 
00221         cos1max = cos(CUTS.phot1cut*CTES.pi/180.0);     // LAB rest frame            
00222       } else {
00223         cos1min = -1.0;                         
00224         cos1max =  1.0;
00225         CUTS.gmin = CUTS.gmin/2.0;
00226       }
00227       cos2min = -1.0;                        // photon2 angle limits      
00228       cos2max =  1.0;                        //                           
00229       cos3min = -1.0;                        // hadrons/muons angle limits    
00230       cos3max =  1.0;                        // in their rest frame            
00231       if (FLAGS.pion == 0)                   // virtual photon energy cut 
00232         qqmin = 4.0*CTES.mmu*CTES.mmu;
00233       else if (FLAGS.pion == 1) 
00234         qqmin = 4.0*CTES.mpi*CTES.mpi;
00235       else if (FLAGS.pion == 2)  
00236         qqmin = 4.0*(CTES.mpi+CTES.mpi0)*(CTES.mpi+CTES.mpi0);
00237       else if (FLAGS.pion == 3)  
00238         qqmin = 16.0*CTES.mpi*CTES.mpi;
00239       else if (FLAGS.pion == 4)  
00240         qqmin = 4.0*CTES.mp*CTES.mp;
00241       else if (FLAGS.pion == 5)  
00242         qqmin = 4.0*CTES.mnt*CTES.mnt;
00243       else if (FLAGS.pion == 6)  
00244         qqmin = 4.0*CTES.mKp*CTES.mKp;
00245       else if (FLAGS.pion == 7)  
00246         qqmin = 4.0*CTES.mKn*CTES.mKn;
00247       else if (FLAGS.pion == 8) 
00248         qqmin = (2.0*CTES.mpi+CTES.mpi0)*(2.0*CTES.mpi+CTES.mpi0);
00249       else if (FLAGS.pion == 9) 
00250         qqmin = 4.0*CTES.mlamb*CTES.mlamb;
00251       //      else
00252       //        continue;
00253 
00254       qqmax = CTES.Sp-2.0*sqrt(CTES.Sp)*CUTS.gmin;         // if only one photon 
00255       if (CUTS.q2_max_c < qqmax) 
00256         qqmax=CUTS.q2_max_c;      // external cuts      
00257 
00258       // -------------------
00259       if ( (CUTS.q2_min_c > qqmin) && (CUTS.q2_min_c < (CTES.Sp*(1.0-2.0*(CUTS.gmin/sqrt(CTES.Sp)+CUTS.w))) ) )
00260         qqmin = CUTS.q2_min_c;
00261       else {
00262         cerr << "------------------------------" << endl;
00263         cerr << " Q^2_min TOO SMALL" << endl;
00264         cerr << " Q^2_min CHANGE BY PHOKHARA = " << qqmin << " GeV^2" << endl;
00265         cerr << "------------------------------" << endl;
00266       }
00267       // -------------------
00268       if(qqmax <= qqmin){
00269         cerr << " Q^2_max to small " << endl;
00270         cerr << " Q^2_max = " << qqmax << endl;
00271         cerr << " Q^2_min = " << qqmin << endl;
00272         return 0;
00273       }
00274 
00275       // -------------------
00276       if (FLAGS.pion == 0) {
00277         printf(" %s %f %s\n", "minimal muon-pair invariant mass^2     = ", qqmin," GeV^2");
00278         printf(" %s %f %s\n", "maximal muon-pair invariant mass^2     = ", qqmax," GeV^2");
00279       } else if (FLAGS.pion == 1) {
00280         printf(" %s %f %s\n", "minimal pion-pair invariant mass^2     = ", qqmin," GeV^2");
00281         printf(" %s %f %s\n", "maximal pion-pair invariant mass^2     = ", qqmax," GeV^2");
00282       } else if (FLAGS.pion == 4) {
00283         printf(" %s %f %s\n", "minimal proton-pair invariant mass^2   = ", qqmin," GeV^2");
00284         printf(" %s %f %s\n", "maximal proton-pair invariant mass^2   = ", qqmax," GeV^2");
00285       } else if (FLAGS.pion == 5) {
00286         printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2  = ", qqmin," GeV^2");
00287         printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2  = ", qqmax," GeV^2");
00288       } else if ((FLAGS.pion == 6) || (FLAGS.pion == 7)) {
00289         printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2     = ", qqmin," GeV^2");
00290         printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2     = ", qqmax," GeV^2");
00291       } else if(FLAGS.pion == 8){
00292         printf(" %s %f %s\n", "minimal three-pion invariant mass^2    = ", qqmin," GeV^2");
00293         printf(" %s %f %s\n", "maximal three-pion invariant mass^2    = ", qqmax," GeV^2");
00294       } else if(FLAGS.pion == 9){
00295         printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2  = ", qqmin," GeV^2");
00296         printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2  = ", qqmax," GeV^2");
00297       } else {
00298         printf(" %s %f %s\n", "minimal four-pion invariant mass^2     = ", qqmin," GeV^2" );
00299         printf(" %s %f %s\n", "maximal four-pion invariant mass^2     = ", qqmax," GeV^2");
00300       }
00301       // -------------------
00302       if (FLAGS.nlo == 0) { 
00303         cout <<  "Born" << endl;
00304          if(FLAGS.fsrnlo != 0){
00305            cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
00306            return 0;
00307          }
00308       }
00309       // -------------------
00310       if((FLAGS.pion == 9) && (FLAGS.nlo != 0)) {
00311         cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
00312         cerr << "If you feel that you need NLO"<< endl;
00313         cerr << "please contact the authors"<< endl;
00314         return 0;
00315       }
00316       // -------------------
00317       if (FLAGS.nlo == 1) printf(" %s %f\n", "NLO:    soft photon cutoff w           = ",CUTS.w);
00318       if ((FLAGS.pion <= 1) || (FLAGS.pion == 6)) 
00319         {
00320           
00321           if( ! ((FLAGS.fsr == 1) || (FLAGS.fsr == 2) || (FLAGS.fsrnlo == 0) 
00322               || (FLAGS.fsr == 1) || (FLAGS.fsrnlo == 1)
00323                  || (FLAGS.fsr == 0) || (FLAGS.fsrnlo == 0))) {
00324             cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
00325             return 0;
00326           }
00327           
00328           // ------------------
00329           if (FLAGS.fsr == 0) 
00330             cout << "ISR only" << endl;
00331           else if (FLAGS.fsr == 1) 
00332             cout << "ISR+FSR" << endl;
00333           else if (FLAGS.fsr == 2) {
00334             if (FLAGS.nlo == 0) 
00335               cout << "ISR+INT+FSR" << endl;
00336             else {
00337               cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
00338               return 0;
00339             }
00340           }
00341           else {
00342             cerr << "WRONG FSR flag" << FLAGS.fsr << endl;
00343             return 0;
00344           }
00345       
00346           if(FLAGS.fsrnlo == 1) 
00347             cout << "IFSNLO included" << endl;
00348         }
00349       else
00350         {
00351           if((FLAGS.fsr == 0) && (FLAGS.fsrnlo == 0))
00352             cout << "ISR only" << endl;
00353           else {
00354             cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
00355             return 0;
00356           }
00357         }
00358 
00359       // ------------------
00360       if(FLAGS.ivac == 0){
00361         cout << "Vacuum polarization is NOT included" << endl;
00362       } else if(FLAGS.ivac == 1){
00363         cout << "Vacuum polarization is included"  << endl;
00364       } else {
00365         cout << "WRONG vacuum polarization switch" << endl;
00366         return 0;
00367       }
00368 // -----------------
00369       if(FLAGS.pion == 1){
00370         if(FLAGS.FF_pion == 0)
00371           cout << "Kuhn-Santamaria PionFormFactor" << endl;
00372         else if(FLAGS.FF_pion == 1)
00373           cout << "Gounaris-Sakurai PionFormFactor" << endl;
00374         else {
00375           cout << "WRONG PionFormFactor switch" << endl;
00376           return 0;
00377         }
00378 // ------
00379         if(FLAGS.fsr != 0){
00380          if(FLAGS.f0_model == 0)
00381            cout << "f0+f0(600): K+K- model" << endl;
00382          else if(FLAGS.f0_model == 1)
00383            cout << "f0+f0(600): \"no structure\" model" << endl;
00384          else if(FLAGS.f0_model == 2)
00385            cout << "NO f0+f0(600)" << endl;
00386          else if(FLAGS.f0_model == 3)
00387            cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
00388          else {
00389            cout << "WRONG f0+f0(600) switch" << endl;
00390            return 0;      
00391          }
00392         }
00393       }
00394 
00395   cout << "--------------------------------------------------------------" << endl;
00396 //
00397  // =================================================
00398 // --- finding the maximum -------------------------
00399       for( i = 0; i<2; i++)
00400         {
00401           MAXIMA.Mmax[i] = 1.0;
00402           MAXIMA.gross[i] = 0.0;
00403           MAXIMA.klein[i] = 0.0;      
00404         }
00405       
00406       if (FLAGS.nlo == 0) 
00407         MAXIMA.Mmax[1]=0.0;   // only 1 photon events generated
00408       
00409       MAXIMA.tr[0] = 0.0;      
00410       MAXIMA.tr[1] = 0.0;
00411       MAXIMA.count[0] = 0.0;
00412       MAXIMA.count[1] = 0.0;
00413       
00414       // =================================================
00415       // --- beginning the MC loop event generation ------
00416       for(j = 1; j <= m_nm; j++)
00417         {
00418           RANLXDF(Ar_r,1);
00419           Ar[1] = Ar_r[0];
00420           
00421           if (Ar[1] <= (MAXIMA.Mmax[0]/(MAXIMA.Mmax[0]+MAXIMA.Mmax[1]))) { 
00422             MAXIMA.count[0] = MAXIMA.count[0]+1.0;
00423             GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
00424           } 
00425           else {
00426             MAXIMA.count[1] = MAXIMA.count[1]+1.0;
00427             GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
00428           }
00429         }
00430       // --- end of the MC loop --------------------------
00431       // =================================================
00432       // --- for the second run ---
00433       
00434       MAXIMA.Mmax[0] = MAXIMA.gross[0]+.05*sqrt(MAXIMA.gross[0]*MAXIMA.gross[0]);
00435       MAXIMA.Mmax[1] = MAXIMA.gross[1]+(.03+.02*CTES.Sp)*sqrt(MAXIMA.gross[1]*MAXIMA.gross[1]); 
00436       
00437       if((FLAGS.pion == 1) && (FLAGS.fsrnlo == 1)) 
00438         MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.5;
00439       if((FLAGS.pion == 0) && (FLAGS.fsrnlo == 1)) 
00440         MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.5;
00441       
00442       if((FLAGS.pion == 0) && (FLAGS.fsr == 1) && (FLAGS.fsrnlo == 0)) 
00443         MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.2;
00444       
00445       if((FLAGS.pion == 2) || (FLAGS.pion == 3)){
00446         MAXIMA.Mmax[0]=MAXIMA.Mmax[0]*1.1;
00447         MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.1;
00448       }
00449       
00450       if(FLAGS.pion == 8){
00451         MAXIMA.Mmax[0]=MAXIMA.Mmax[0]*1.08;
00452         MAXIMA.Mmax[1]=MAXIMA.Mmax[1]*1.1;
00453       }
00454 // --- end of the second run -----------------------
00455       
00456       MAXIMA.tr[0] = 0.0;      
00457       MAXIMA.tr[1] = 0.0;
00458       MAXIMA.count[0] = 0.0;
00459       MAXIMA.count[1] = 0.0;
00460       
00461   return StatusCode::SUCCESS;
00462 }

StatusCode Phokhara::storeParticles (   ) 
 

StatusCode Phokhara::storeParticles (   ) 
  

00505 {
00506   MsgStream log(messageService(), name());
00507   int npart = 0;
00508   
00509   // Fill event information
00510   GenEvent* evt = new GenEvent(1,1);
00511 
00512   GenVertex* prod_vtx = new GenVertex();
00513 //  prod_vtx->add_particle_out( p );
00514   evt->add_vertex( prod_vtx );
00515 
00516   // incoming beam e+
00517   GenParticle* p = new GenParticle( HepLorentzVector( CTES.momenta[1][0], CTES.momenta[2][0], CTES.momenta[3][0], CTES.momenta[0][0] ), -11, 3); 
00518   p->suggest_barcode( ++npart );
00519   prod_vtx->add_particle_in(p);
00520 
00521   // incoming beam e-
00522   p = new GenParticle( HepLorentzVector( CTES.momenta[1][1], CTES.momenta[2][1], CTES.momenta[3][1], CTES.momenta[0][1] ), 11, 3); 
00523   p->suggest_barcode( ++npart );
00524   prod_vtx->add_particle_in(p);
00525 
00526 
00527   // gamma
00528   p = new GenParticle( HepLorentzVector( CTES.momenta[1][2], CTES.momenta[2][2], CTES.momenta[3][2], CTES.momenta[0][2] ), 22, 1); 
00529   p->suggest_barcode( ++npart );
00530   prod_vtx->add_particle_out(p);
00531   
00532   if( CTES.momenta[0][3] != 0 ) // second photon exists
00533     {
00534       p = new GenParticle( HepLorentzVector( CTES.momenta[1][3], CTES.momenta[2][3], CTES.momenta[3][3], CTES.momenta[0][3] ), 22, 1); 
00535       p->suggest_barcode( ++npart );
00536       prod_vtx->add_particle_out(p);
00537     }
00538   
00539   // other products depending on channel
00540   if (FLAGS.pion == 0) { // mu+ mu-
00541     // mu+
00542     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), -13, 1); 
00543     p->suggest_barcode( ++npart );
00544     prod_vtx->add_particle_out(p);
00545     // mu -
00546     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 13, 1); 
00547     p->suggest_barcode( ++npart );
00548     prod_vtx->add_particle_out(p);
00549   }    
00550     
00551   if (FLAGS.pion == 1) { // pi+ pi-
00552     // pi+
00553     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),211, 1); 
00554     p->suggest_barcode( ++npart );
00555     prod_vtx->add_particle_out(p);
00556     // pi-
00557     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -211, 1); 
00558     p->suggest_barcode( ++npart );
00559     prod_vtx->add_particle_out(p);
00560   } 
00561   
00562   if (FLAGS.pion == 2) { // pi+ pi- pi0 pi0
00563     // pi0
00564     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),111, 1); 
00565     p->suggest_barcode( ++npart );
00566     prod_vtx->add_particle_out(p);
00567     // pi0
00568     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 111, 1); 
00569     p->suggest_barcode( ++npart );
00570     prod_vtx->add_particle_out(p);
00571     // pi-
00572     p = new GenParticle( HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ),-211, 1); 
00573     p->suggest_barcode( ++npart );
00574     prod_vtx->add_particle_out(p);
00575     // pi+
00576     p = new GenParticle( HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), 211, 1); 
00577     p->suggest_barcode( ++npart );
00578     prod_vtx->add_particle_out(p);
00579   } 
00580   
00581   if (FLAGS.pion == 3) { // pi+ pi- pi0 pi0
00582     // pi0
00583     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),111, 1); 
00584     p->suggest_barcode( ++npart );
00585     prod_vtx->add_particle_out(p);
00586     // pi0
00587     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 111, 1); 
00588     p->suggest_barcode( ++npart );
00589     prod_vtx->add_particle_out(p);
00590     // pi-
00591     p = new GenParticle( HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ),-211, 1); 
00592     p->suggest_barcode( ++npart );
00593     prod_vtx->add_particle_out(p);
00594     // pi+
00595     p = new GenParticle( HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), 211, 1); 
00596     p->suggest_barcode( ++npart );
00597     prod_vtx->add_particle_out(p);
00598   } 
00599 
00600   if (FLAGS.pion == 4) { // p pbar
00601     // pbar
00602     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),-2212, 1); 
00603     p->suggest_barcode( ++npart );
00604     prod_vtx->add_particle_out(p);
00605     // p
00606     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 2212, 1); 
00607     p->suggest_barcode( ++npart );
00608     prod_vtx->add_particle_out(p);
00609   } 
00610 
00611   if (FLAGS.pion == 5) { // n nbar
00612     // nbar
00613     p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),-2112, 1); 
00614     p->suggest_barcode( ++npart );
00615     prod_vtx->add_particle_out(p);
00616     // n
00617     p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 2112, 1); 
00618     p->suggest_barcode( ++npart );
00619     prod_vtx->add_particle_out(p);
00620   } 
00621   
00622      if (FLAGS.pion == 6) { // K+ K-
00623        // pbar
00624        p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ),321, 1); 
00625        p->suggest_barcode( ++npart );
00626        prod_vtx->add_particle_out(p);
00627        // p
00628        p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -321, 1); 
00629        p->suggest_barcode( ++npart );
00630        prod_vtx->add_particle_out(p);
00631      } 
00632      
00633      if (FLAGS.pion == 7) { // K0 K0bar
00634        // pbar
00635        p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), 311, 1); 
00636        p->suggest_barcode( ++npart );
00637        prod_vtx->add_particle_out(p);
00638        // p
00639        p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -311, 1); 
00640        p->suggest_barcode( ++npart );
00641        prod_vtx->add_particle_out(p);
00642      } 
00643 
00644      if (FLAGS.pion == 8) { // pi+ pi- pi0
00645        // pi+
00646        p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), 211, 1); 
00647        p->suggest_barcode( ++npart );
00648        prod_vtx->add_particle_out(p);
00649        // pi-
00650        p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), -211, 1); 
00651        p->suggest_barcode( ++npart );
00652        prod_vtx->add_particle_out(p);
00653        // pi0
00654        p = new GenParticle( HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ), 111, 1); 
00655        p->suggest_barcode( ++npart );
00656        prod_vtx->add_particle_out(p);
00657      } 
00658 
00659      if (FLAGS.pion == 9) { // Lambda Lambdabar Pi+ Pi- P Pbar
00660        // Lambdabar
00661        p = new GenParticle( HepLorentzVector( CTES.momenta[1][5], CTES.momenta[2][5], CTES.momenta[3][5], CTES.momenta[0][5] ), -3122, 2); 
00662        p->suggest_barcode( ++npart );
00663        prod_vtx->add_particle_out(p);
00664        // Lambda
00665        p = new GenParticle( HepLorentzVector( CTES.momenta[1][6], CTES.momenta[2][6], CTES.momenta[3][6], CTES.momenta[0][6] ), 3122, 2); 
00666        p->suggest_barcode( ++npart );
00667        prod_vtx->add_particle_out(p);
00668        // pi+
00669        p = new GenParticle( HepLorentzVector( CTES.momenta[1][7], CTES.momenta[2][7], CTES.momenta[3][7], CTES.momenta[0][7] ), 211, 1); 
00670        p->suggest_barcode( ++npart );
00671        prod_vtx->add_particle_out(p);
00672        // pbar
00673        p = new GenParticle( HepLorentzVector( CTES.momenta[1][8], CTES.momenta[2][8], CTES.momenta[3][8], CTES.momenta[0][8] ), -2212, 1); 
00674        p->suggest_barcode( ++npart );
00675        prod_vtx->add_particle_out(p);
00676        // pi-
00677        p = new GenParticle( HepLorentzVector( CTES.momenta[1][9], CTES.momenta[2][9], CTES.momenta[3][9], CTES.momenta[0][9] ), -211, 1); 
00678        p->suggest_barcode( ++npart );
00679        prod_vtx->add_particle_out(p);
00680        // p
00681        p = new GenParticle( HepLorentzVector( CTES.momenta[1][10], CTES.momenta[2][10], CTES.momenta[3][10], CTES.momenta[0][10] ), 2212, 1); 
00682        p->suggest_barcode( ++npart );
00683        prod_vtx->add_particle_out(p);
00684      } 
00685 
00686   if( log.level() < MSG::INFO )
00687   {
00688     evt->print();  
00689   }
00690 
00691   // Check if the McCollection already exists
00692   SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
00693   if (anMcCol!=0) 
00694   {
00695     // Add event to existing collection
00696     MsgStream log(messageService(), name());
00697     log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
00698     McGenEvent* mcEvent = new McGenEvent(evt);
00699     anMcCol->push_back(mcEvent);
00700   }
00701   else 
00702   {
00703     // Create Collection and add  to the transient store
00704     McGenEventCol *mcColl = new McGenEventCol;
00705     McGenEvent* mcEvent = new McGenEvent(evt);
00706     mcColl->push_back(mcEvent);
00707     StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
00708     if (sc != StatusCode::SUCCESS) 
00709     {
00710       log << MSG::ERROR << "Could not register McGenEvent" << endreq;
00711       delete mcColl;
00712       delete evt;
00713       delete mcEvent;
00714       return StatusCode::FAILURE;
00715     }
00716     else 
00717     {
00718       log << MSG::INFO << "McGenEventCol created and " << npart <<" particles stored in McGenEvent" << endreq;
00719     }
00720   }
00721 
00722   return StatusCode::SUCCESS; 
00723 
00724 }

Member Data Documentation

double Phokhara::Ar [private]
 

double Phokhara::Ar_r [private]
 

double Phokhara::cos1max [private]
 

double Phokhara::cos1min [private]
 

double Phokhara::cos2max [private]
 

double Phokhara::cos2min [private]
 

double Phokhara::cos3max [private]
 

double Phokhara::cos3min [private]
 

double Phokhara::dsigm [private]
 

double Phokhara::dsigm1 [private]
 

double Phokhara::dsigm2 [private]
 

int Phokhara::ievent [private]
 

double Phokhara::m_E [private]
 

int Phokhara::m_f0_model [private]
 

int Phokhara::m_FF_Pion [private]
 

int Phokhara::m_fsr [private]
 

int Phokhara::m_fsrnlo [private]
 

double Phokhara::m_gmin [private]
 

long int Phokhara::m_initSeed [private]
 

int Phokhara::m_ivac [private]
 

int Phokhara::m_nlo [private]
 

int Phokhara::m_nm [private]
 

double Phokhara::m_phot1cut [private]
 

double Phokhara::m_phot2cut [private]
 

double Phokhara::m_pi1cut [private]
 

double Phokhara::m_pi2cut [private]
 

int Phokhara::m_pion [private]
 

double Phokhara::m_q2_max_c [private]
 

double Phokhara::m_q2_min_c [private]
 

double Phokhara::m_q2min [private]
 

int Phokhara::m_tagged [private]
 

double Phokhara::m_w [private]
 

double Phokhara::qqmax [private]
 

double Phokhara::qqmin [private]
 

double Phokhara::sigma [private]
 

double Phokhara::sigma1 [private]
 

double Phokhara::sigma2 [private]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:40:30 2011 for BOSS6.5.5 by  doxygen 1.3.9.1