Ekhara Class Reference

#include <Ekhara.h>

Inheritance diagram for Ekhara:

Public Member Functions
	Ekhara (const std::string &name, ISvcLocator *pSvcLocator)
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
StatusCode	checkAndPrintParameters ()

Detailed Description

Definition at line 15 of file Ekhara.h.

Constructor & Destructor Documentation

Ekhara()

Ekhara::Ekhara	(	const std::string &	name,
		ISvcLocator *	pSvcLocator )

• CMS-energy (GeV)

• # of events to determine the max xs

• handle UB underestim <0> = halt,
  • <positive> = use new UB if overshoot,
  • <negative> = ignore

• Enlarge factor for upper bound

• 0:pi0pi0, 1:pi+pi-, 2:pi0, 3:eta, 4:eta', 5:Chi_c0, 6:Chi_c1, 7:Chi_c2

• pi0,eta,eta': s(1); t(2); s+t(3); t:signal

• include vacuum polarization (true or false)

• NLO radiative corrections (true or false)

• return weighted events as output for NLO (true or false)

• phase space regulator eps_ph: mgamma=me*eps_ph

• pi0: WZWconst(1); rho pole(2);
  • LMD(3); LMD+V(4); LMD+Vnew(5);
  • 1-Octet(7); 2-Octet(8);
  • eta/eta': WZWconst(1); VMD(3);
  • 1-Octet(7); 2-Octet(8);

• e+ minimum angle

• e+ maximum angle

• e+ minimum Energy

• e+ maximum Energy

• e- minimum angle

• e- maximum angle

• e- minimum Energy

• e- maximum Energy

• pi0, eta, eta' minimum angle

• pi0, eta, eta' maximum angle

• pi0, eta, eta' minimum Energy

• pi0, eta, eta' maximum Energy

• the choice of generation algorithm: 0-v.1.0; 1-v.3.1

• pi+pi-: s(1); s+t(2); s+t+2g_Born(3); s+t+2g_Full(4)

• the choice of the pion form factor to be used: 0-KS; 1-GS; 2-GS new

• pion angle cut

• missing momentum angle cut

• Chi_cj minimum angle

• Chi_cj maximum angle

• Chi_cj minimum energy

• Chi_cj minimum energy

• Tagging angle (symmetric) in degrees

• minimum Q^2 of tagged lepton in GeV^2

• Tagging mode 0:no tagging; 1/11:untagged
  • 2/12: single tagging; 3/13:double tagging
    

Definition at line 30 of file Ekhara.cxx.

                                                              : Algorithm(name, pSvcLocator)
{
  declareProperty( "Ecm",                 m_cmsEnergy = 3.773 );      //!- CMS-energy (GeV) 
  declareProperty( "InitialEvents",       m_numUpperBound = 50000 );  //!- # of events to determine the max xs
  declareProperty( "IgnoreUpperBound",    m_ignoreUpperBound = 0.0 ); //!- handle UB underestim <0> = halt,
                                                                      //!- <positive> = use new UB if overshoot,
                                                                      //!- <negative> = ignore
  declareProperty( "UpperBoundEnlarge",   m_upperBoundEnlarge = 1.2 );//!- Enlarge factor for upper bound
  declareProperty( "FinalStateID",        m_finalState = 2 );         //!- 0:pi0pi0, 1:pi+pi-, 2:pi0, 3:eta, 4:eta', 5:Chi_c0, 6:Chi_c1, 7:Chi_c2 
 
  declareProperty( "MesonProdAmplitudes", m_mesonAmplitudes = 2 );    //!- pi0,eta,eta': s(1); t(2); s+t(3); t:signal
  declareProperty( "IncludeVacuumPolarization", m_switchVP = false ); //!- include vacuum polarization (true or false)
  declareProperty( "ConsiderNLO",         m_switchNLO = false );      //!- NLO radiative corrections (true or false)
  declareProperty( "OutputWeightedNLOEvents",  m_nloWithWeights = false );      //!- return weighted events as output for NLO (true or false)
  declareProperty( "PhaseSpaceRegulator", m_eps_ph = 0.01 );          //!- phase space regulator eps_ph: mgamma=me*eps_ph
  declareProperty( "MesonFormfactor",     m_TFF_ID = 8 );             //!- pi0: WZWconst(1); rho pole(2);
                                                                      //!-      LMD(3); LMD+V(4); LMD+Vnew(5);
                                                                      //!-      1-Octet(7); 2-Octet(8);
                                                                      //!- eta/eta': WZWconst(1); VMD(3);
                                                                      //!-           1-Octet(7); 2-Octet(8);
  declareProperty( "PositronThetaMin",    m_posiThetaMin   = 0.0 );   //!- e+ minimum angle 
  declareProperty( "PositronThetaMax",    m_posiThetaMax   = 180.0 ); //!- e+ maximum angle
  declareProperty( "PositronEnergyMin",   m_posiEnergyMin  = 0.0 );   //!- e+ minimum Energy
  declareProperty( "PositronEnergyMax",   m_posiEnergyMax  = 110.0 ); //!- e+ maximum Energy
  declareProperty( "ElectronThetaMin",    m_elecThetaMin   = 0.0 );   //!- e- minimum angle
  declareProperty( "ElectronThetaMax",    m_elecThetaMax   = 180.0 ); //!- e- maximum angle
  declareProperty( "ElectronEnergyMin",   m_elecEnergyMin  = 0.0 );   //!- e- minimum Energy
  declareProperty( "ElectronEnergyMax",   m_elecEnergyMax  = 110.0 ); //!- e- maximum Energy
 
  declareProperty( "MesonThetaMin",       m_mesonThetaMin  = 0.0 );   //!- pi0, eta, eta' minimum angle
  declareProperty( "MesonThetaMax",       m_mesonThetaMax  = 180.0 ); //!- pi0, eta, eta' maximum angle
  declareProperty( "MesonEnergyMin",      m_mesonEnergyMin = 0.0 );   //!- pi0, eta, eta' minimum Energy
  declareProperty( "MesonEnergyMax",      m_mesonEnergyMax = 100.0 ); //!- pi0, eta, eta' maximum Energy
 
  declareProperty( "TwoPionPhaseSpaceAlg",  m_twoPionPhspAlg = 1 );         //!- the choice of generation algorithm: 0-v.1.0; 1-v.3.1
  declareProperty( "TwoPionProdAmplitudes", m_twoPionAmplitudes = 4 );      //!- pi+pi-:  s(1); s+t(2); s+t+2g_Born(3); s+t+2g_Full(4)
  declareProperty( "TwoPionFormFactor"   ,  m_twoPionFormFactor = 2 );      //!- the choice of the pion form factor to be used: 0-KS; 1-GS; 2-GS new
  declareProperty( "TwoPionThetaMin",       m_twoPionThetaMin     = 0.0);   //!- pion angle cut
  declareProperty( "TwoPionThetaMax",       m_twoPionThetaMax     = 180.0);
  declareProperty( "TwoPionMissThetaMin",   m_twoPionMissThetaMin = 0.0);   //!- missing momentum angle cut
  declareProperty( "TwoPionMissThetaMax",   m_twoPionMissThetaMax = 180.0);
 
  declareProperty( "Chi_cjThetaMin",     m_chicjThetaMin   = 0.0);    //!- Chi_cj minimum angle
  declareProperty( "Chi_cjThetaMax",     m_chicjThetaMax   = 180.0);  //!- Chi_cj maximum angle
  declareProperty( "Chi_cjEnergyMin",    m_chicjEnergyMin  = 0.0);    //!- Chi_cj minimum energy
  declareProperty( "Chi_cjEnergyMax",    m_chicjEnergyMax  = 100.0);  //!- Chi_cj minimum energy
 
  declareProperty( "TaggingAngle",        m_taggingAngle   = 20. );   //!- Tagging angle (symmetric) in degrees 
  declareProperty( "TaggingQsquare",      m_taggingQsquare = 0.3 );   //!- minimum Q^2 of tagged lepton in GeV^2 
  declareProperty( "TaggingMode",         m_taggingMode    = 0 );     //!- Tagging mode 0:no tagging; 1/11:untagged
                                                                      //!-              2/12: single tagging; 3/13:double tagging  
 
 
}

Member Function Documentation

checkAndPrintParameters()

StatusCode Ekhara::checkAndPrintParameters

(

)

• pi0pi0

• Check switches for NLO

• pi+pi-

• Check switches for NLO

• 1 - s, 2 - s+t, 3 - s+t+g_Born, 4 - s+t+g_Full

• 0 - KS, 1 - GS, 2 - new GS

-Single Pseudoscalar Meson

• Check switches for NLO and VacPol

• 1 - s, 2 - t, 3 - s+t

• 1:Untagged , 2:Single, 3: double

• Chi_c0, Chi_c1, Chi_c2

• Check switches for NLO and VacPol

Definition at line 513 of file Ekhara.cxx.

{
  MsgStream log(msgSvc(), name());
 
  log << MSG::INFO << "=============================================================" << endreq;
  log << MSG::INFO << "This is EKHARA, Version 3.1" << endreq;
  switch(m_finalState) {
  case 0:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi0pi0" << endreq;
    break;
  case 1:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi+pi-" << endreq;
    break;
  case 2:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- pi0" << endreq;
    break;
  case 3:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- eta" << endreq;
    break;
  case 4:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- etaPRIME" << endreq;
    break;
  case 5:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c0" << endreq;
    break;
  case 6:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c1" << endreq;
    break;
  case 7:
    log << MSG::INFO << "\tSimulating the process:\te+e- -> e+e- Chi_c2" << endreq;
    break;
  default:
    log << MSG::ERROR << "WRONG channel ID: Process not implemented!" << endreq;
    return StatusCode::FAILURE;
    break;
  }
 
  
  if ( m_finalState == 0 ) {  //!- pi0pi0
    if( m_switchNLO ) {//!- Check switches for NLO
      log << MSG::ERROR << "NLO not implemented for this Process!" << endreq;
      return StatusCode::FAILURE;
    }
    if (m_twoPionAmplitudes!=4){
      log << MSG::ERROR << " Neutral Pion Pairs only with full TwoPionAmplitudes == 4 ! \t You used " << m_twoPionAmplitudes << endreq;
      return StatusCode::FAILURE;
    }
    if (m_twoPionPhspAlg!=1){
      log << MSG::ERROR << " Neutral Pion Pairs only with full TwoPionPhaseSpaceAlg == 1 ! \t You used " << m_twoPionPhspAlg << endreq;
      return StatusCode::FAILURE;
    }
    log << MSG::INFO << "\tPhase space generation using algorithm: " << m_twoPionPhspAlg << endreq;
    log << MSG::INFO << "\tMatrix element: \t|M_s + M_t + M_2g(full)|^2" << endreq;
 
    log << MSG::INFO << "The following conditions are applied:" << endreq;
    log << MSG::INFO << "\tPion Momentum:    " << m_twoPionThetaMin << " < Theta [deg] < "     << m_twoPionThetaMax << endreq;
    log << MSG::INFO << "\tMissing Momentum: " << m_twoPionMissThetaMin << " < Theta [deg] < " << m_twoPionMissThetaMax << endreq;
 
  } else if (m_finalState==1) {//!- pi+pi-
    if( m_switchNLO ) {//!- Check switches for NLO
      log << MSG::ERROR << "NLO not implemented for this Process!" << endreq;
      return StatusCode::FAILURE;
    }
    if (m_twoPionPhspAlg<0 || m_twoPionPhspAlg>1){
      log << MSG::ERROR << " Wrong choice of phase space generation algorithm: " << m_twoPionPhspAlg << endreq;
      return StatusCode::FAILURE;
    }
    
    log << MSG::INFO << "\tPhase space generation using algorithm: " << m_twoPionPhspAlg << endreq;
    log << MSG::INFO << "\tMatrix element: ";
    switch(m_twoPionAmplitudes){ //!-  1 - s, 2 - s+t, 3 - s+t+g_Born, 4 - s+t+g_Full
    case 1:
      log << MSG::INFO << "\t|M_s|^2" << endreq;
      break;
    case 2:
      log << MSG::INFO << "\t|M_s + M_t|^2" << endreq;
      break;
    case 3:
      log << MSG::INFO << "\t|M_s + M_t + M_2g(Born)|^2" << endreq;
      break;
    case 4:
      log << MSG::INFO << "\t|M_s + M_t + M_2g(full)|^2" << endreq;
      break;
    default:
      break;
    }
 
    switch(m_twoPionFormFactor){//!-  0 - KS, 1 - GS, 2 - new GS
      case 0:
        log << MSG::INFO << "\tUsing KS pion form factor" << endreq;
        break;
      case 1:
        log << MSG::INFO << "\tUsing GS pion form factor" << endreq;
        break;
      case 2:
        log << MSG::INFO << "\tUsing new GS pion form factor" << endreq;
        break;
      default:
        log << MSG::WARNING << "Undefined pion form factor switch! Using new GS pion form factor instead!" << endreq;
        m_twoPionFormFactor = 2;
        break;
    }
 
    log << MSG::INFO << "The following conditions are applied:" << endreq;
    log << MSG::INFO << "\tPion Momentum:    " << m_twoPionThetaMin << " < Theta [deg] < "  << m_twoPionThetaMax << endreq;
    log << MSG::INFO << "\tMissing Momentum: " << m_twoPionMissThetaMin << " < Theta [deg] < "  << m_twoPionMissThetaMax << endreq;
 
  } else if ( m_finalState >= 2 && m_finalState <= 4 ) { //!-Single Pseudoscalar Meson
    if( m_switchNLO ) {//!- Check switches for NLO and VacPol
      log << MSG::INFO << "\tWith NLO corrections: mgamma = " << m_eps_ph * 0.51099906E-3 << endreq;
      if( m_nloWithWeights ) {
        log << MSG::WARNING << "Output will contain weighted events!" << endreq;
      } else {
        log << MSG::WARNING << "Attention! Events will be unweighted for output!" << endreq;
      }
      if( m_switchVP ) {
        log << MSG::INFO << "\tVacuum polarization included:" << endreq;
        log << MSG::INFO << "\thttp://www-com.physik.hu-berlin.de/~fjeger/software.html" << endreq;
      } else  {
        log << MSG::INFO << "\tVacuum polarization NOT included:" << endreq;
      }
    } else {
      log << MSG::INFO << "\tWithout NLO corrections!" << endreq;
      if( m_switchVP ) {
        log << MSG::ERROR << "Vacuum polarization can ONLY be included for NLO!" << endreq;
        return StatusCode::FAILURE;
      }
    }
  
    log << MSG::INFO << "The cross section will be calculated" << endreq;
    log << MSG::INFO << "\tMatrix element: ";
    switch(m_mesonAmplitudes){ //!-  1 - s, 2 - t, 3 - s+t
    case 1:
      log << MSG::INFO << "\t|M_s|^2" << endreq;
      if ( m_switchNLO==1 ) {
          log << MSG::ERROR << " NLO not implemented for s-channel!" << endreq;
          return StatusCode::FAILURE;
      }
      if( m_TFF_ID == 7 || m_TFF_ID == 8 ) {
          log << MSG::WARNING << "TFF Models 7 and 8 were not compared to data in the only-s-channel configuration!" << endreq;
      }
      break;
    case 2:
      log << MSG::INFO << "\t|M_t|^2" << endreq;
      break;
    case 3:
      log << MSG::INFO << "\t|M_s + M_t|^2" << endreq;
      break;
    default:
      break;
    }
 
    log << MSG::INFO << "The following conditions are applied:" << endreq;
    log << MSG::INFO << "\te+   : " << m_posiThetaMin << " < Theta [deg] < " << m_posiThetaMax << endreq; 
    log << MSG::INFO << "\t       " << m_posiEnergyMin << " < Energy [GeV] < " << m_posiEnergyMin << endreq;
    log << MSG::INFO << "\te-   : " << m_elecThetaMin << " < Theta [deg] < " << m_elecThetaMax  << endreq;
    log << MSG::INFO << "\t       " << m_elecEnergyMin << " < Energy [GeV] < " << m_elecEnergyMin << endreq;
    
    double cosTagAngleRad = fabs(cos(m_taggingAngle*TMath::DegToRad()));
    switch( m_taggingMode ){ //!-  1:Untagged , 2:Single, 3: double
    case 1:
      log << MSG::INFO << "\t Generating Untagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with: |cos(Theta Lepton)| > " << cosTagAngleRad << "!" << endreq;
      break;
    case 2:
      log << MSG::INFO << "\t Generating Single Tagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with:" << endreq;
      log << MSG::INFO << "\t\t |cos(Theta e+/-)| > " << cosTagAngleRad  
      << " and |cos(Theta e-/+)| < " << cosTagAngleRad << endreq;
      break;
    case 3:
      log << MSG::INFO << "\t Generating Double Tagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with: |cos(Theta Lepton)| < " << cosTagAngleRad << "!" << endreq;
      break;
    case 11:
      log << MSG::INFO << "\t Generating Untagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with: Q^2 (Lepton)| < " << m_taggingQsquare << "!" << endreq;
      break;
    case 12:
      log << MSG::INFO << "\t Generating Single Tagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with:" << endreq;
      log << MSG::INFO << "\t\t Q^2 (e+/-) > " << m_taggingQsquare  
      << " and Q^2 (e-/+) < " << m_taggingQsquare << endreq;
      break;
    case 13:
      log << MSG::INFO << "\t Generating Double Tagged event configuration!" << endreq;
      log << MSG::INFO << "\t Accepting only events with: Q^2 (Lepton) > " << m_taggingQsquare << "!" << endreq;
      break;
    default:
      if( m_taggingMode>0 ) {
          log << MSG::INFO << "\t Unknown tagging mode selected!" << endreq;
          log << MSG::INFO << "\t Generating events without tagging configuration!" << endreq;
      }
          m_taggingMode = 0;
      break;
    }
   
    switch( m_finalState ) {
    case 2:
      log << MSG::INFO << "\tpi0  : ";
      break;
    case 3:
      log << MSG::INFO << "\teta  : ";
      break;
    case 4:
      log << MSG::INFO << "\teta' : ";
      break;
    default:
      break;
    }
    log << MSG::INFO << m_mesonThetaMin << " < Theta [deg] < " << m_mesonThetaMax  << endreq;
    log << MSG::INFO << "\t       " << m_mesonEnergyMin << " < Energy [GeV] < " << m_mesonEnergyMin << endreq;
  } else if(m_finalState >= 5 && m_finalState <= 7) { //!- Chi_c0, Chi_c1, Chi_c2
    if( m_switchNLO ) {//!- Check switches for NLO and VacPol
      log << MSG::ERROR << "NLO not implemented for this Process!" << endreq;
      return StatusCode::FAILURE;
    }
    log << MSG::INFO << "The following conditions are applied:" << endreq;
    log << MSG::INFO << "\te+   : " << m_posiThetaMin << " < Theta [deg] < " << m_posiThetaMax << endreq; 
    log << MSG::INFO << "\t       " << m_posiEnergyMin << " < Energy [GeV] < " << m_posiEnergyMin << endreq;
    log << MSG::INFO << "\te-   : " << m_elecThetaMin << " < Theta [deg] < " << m_elecThetaMax  << endreq;
    log << MSG::INFO << "\t       " << m_elecEnergyMin << " < Energy [GeV] < " << m_elecEnergyMin << endreq;
    switch(m_finalState) {
    case 5:
      log << MSG::INFO << "\tchi_c0  : ";
      break;
    case 6:
      log << MSG::INFO << "\tchi_c1  : ";
      break;
    case 7:
      log << MSG::INFO << "\tchi_c2 : ";
      break;
    default:
      break;
    }
    log << MSG::INFO << m_chicjThetaMin << " < Theta [deg] < " << m_chicjThetaMax  << endreq;
    log << MSG::INFO << "\t       " << m_chicjEnergyMin << " < Energy [GeV] < " << m_chicjEnergyMin << endreq;
 
  }
  log << MSG::INFO << "=============================================================" << endreq;
 
 
 
 
 
 
 
  return StatusCode::SUCCESS;
}

Referenced by initialize().

execute()

StatusCode Ekhara::execute

(

)

• Generate an event

-Store weigth as particle at rest:

Definition at line 166 of file Ekhara.cxx.

                           {
  MsgStream log(msgSvc(), name());
  log << MSG::DEBUG << "Ekhara in execute()" << endreq;
 
  //!- Generate an event
  if( m_finalState >= 2 && m_finalState <= 4 && m_switchNLO && m_nloWithWeights ) {
    EKHARA(2);
  } else {
    EKHARA(0);
  }
 
  double p1[4], p2[4], q1[4], q2[4];
  double pi1[4], pi2[4], qpion[4], qcj[4], kphp[4];
  
  for (int i=0;i<4;i++) {
    p1[i]    = 0.0;      //four vector of incoming positron
    p2[i]    = 0.0;      //four vector of incoming electron
    q1[i]    = 0.0;      //four vector of outgoung positron
    q2[i]    = 0.0;      //four vector of outgoing electron
    pi1[i]   = 0.0;      //four vector of produced pi0,pi+ (m_finalState == 0,1)
    pi2[i]   = 0.0;      //four vector of produced pi0,pi- (m_finalState == 0,1)
    qpion[i] = 0.0;      //four vector of produced pi0,eta,eta' (m_finalState == 2,3,4)
    qcj[i]   = 0.0;      //four vector of produced chi_cj (m_finalState == 5,6,7)
    kphp[i]  = 0.0;      //four vector of radiative photon
  }
 
  // Fill event information
  GenEvent* evt = new GenEvent(1,1);
 
  GenVertex* prod_vtx = new GenVertex();
  evt->add_vertex( prod_vtx );
 
  GET_FOURMOMENTA_LEPTONS(p1,p2,q1,q2);
 
  // incoming beam e+
  GenParticle* posi = new GenParticle( CLHEP::HepLorentzVector(p1[1],p1[2],p1[3],p1[0]), -11, 3);
  posi->suggest_barcode( 1 );
  prod_vtx->add_particle_in(posi);
 
  // incoming beam e-
  GenParticle *elec = new GenParticle(CLHEP::HepLorentzVector(p2[1],p2[2],p2[3],p2[0]), 11, 3);
  elec->suggest_barcode( 2 );
  prod_vtx->add_particle_in(elec);
 
  evt->set_beam_particles(posi,elec);
 
  // outcoming beam e+
  GenParticle *p = new GenParticle(CLHEP::HepLorentzVector(q1[1],q1[2],q1[3],q1[0]), -11, 1);
  p->suggest_barcode( 3 );
  prod_vtx->add_particle_out(p);
 
  // outcoming beam e-
  p = new GenParticle( CLHEP::HepLorentzVector(q2[1],q2[2],q2[3],q2[0] ), 11, 1); 
  p->suggest_barcode( 4 );
  prod_vtx->add_particle_out(p);
 
  switch (m_finalState) { 
  case 0: // e+e-pi0 pi0
    GET_FOURMOMENTA_TWOPI(pi1,pi2);
    p = new GenParticle( CLHEP::HepLorentzVector(pi1[1],pi1[2],pi1[3],pi1[0]), 111,1); 
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    p = new GenParticle( CLHEP::HepLorentzVector(pi2[1],pi2[2],pi2[3],pi2[0]), 111,1); 
    p->suggest_barcode( 6 );
    prod_vtx->add_particle_out(p);
    break;
  case 1: // e+e-pi+ pi-
    GET_FOURMOMENTA_TWOPI(pi1,pi2);
    p = new GenParticle( CLHEP::HepLorentzVector(pi1[1],pi1[2],pi1[3],pi1[0]), 211,1); 
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    p = new GenParticle( CLHEP::HepLorentzVector(pi2[1],pi2[2],pi2[3],pi2[0]),-211,1); 
    p->suggest_barcode( 6 );
    prod_vtx->add_particle_out(p);
    double weights[6]; //event weight + five weights for the various contributions to e+e- --> e+e-pi+pi-
    for (int i=0;i<6;i++) {
      weights[i] = 0.0;
    }
    GET_TWOPI_WEIGHTS(weights);
    for (int ai=0; ai<(sizeof(weights)/sizeof(double)); ai++) {
      evt->weights().push_back(weights[ai]);
    }
    p = new GenParticle( CLHEP::HepLorentzVector( 0, 0, 0, weights[0] ), 99, 1); 
    p->suggest_barcode( 7 );
    prod_vtx->add_particle_out(p);
 
    break;
  case 2: //e+e- pi0
    GET_FOURMOMENTA_PION(qpion);
    p = new GenParticle( CLHEP::HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 111, 1); 
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    break;
  case 3: //e+e- eta
    GET_FOURMOMENTA_PION(qpion);
    p = new GenParticle( CLHEP::HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 221, 1); 
    p->suggest_barcode(5);
    prod_vtx->add_particle_out(p);
    break;
  case 4: //e+e- etaPrime
    GET_FOURMOMENTA_PION(qpion);
    p = new GenParticle( CLHEP::HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 331, 1); 
    p->suggest_barcode(5);
    prod_vtx->add_particle_out(p);
    break;
  case 5: //e+e- chi_c0
    GET_FOURMOMENTA_CHICJ(qcj);
    p = new GenParticle( CLHEP::HepLorentzVector( qcj[1],qcj[2],qcj[3],qcj[0]), 10441, 1); 
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    break;
  case 6: //e+e- chi_c1
    GET_FOURMOMENTA_CHICJ(qcj);
    p = new GenParticle( CLHEP::HepLorentzVector( qcj[1],qcj[2],qcj[3],qcj[0]), 20443, 1);
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    break;
  case 7: //e+e- chi_c2
    GET_FOURMOMENTA_CHICJ(qcj);
    p = new GenParticle( CLHEP::HepLorentzVector( qcj[1],qcj[2],qcj[3],qcj[0]), 445, 1); 
    p->suggest_barcode( 5 );
    prod_vtx->add_particle_out(p);
    break;
  default:
    break;
  }
 
 
  if( m_finalState >= 2 && m_finalState <= 4 && m_switchNLO) {
    if(NLOTYPE.EvtPhTyp == 1) {
      GET_FOURMOMENTA_PHOTON(kphp);
      p = new GenParticle( CLHEP::HepLorentzVector(kphp[1],kphp[2],kphp[3],kphp[0]), 22, 1 ); 
      p->suggest_barcode( 6 );
      prod_vtx->add_particle_out(p);
    }
    if( m_nloWithWeights ) {
      double weight = GET_WEIGHT();
      evt->weights().push_back(weight);
      //!-Store weigth as particle at rest:
      if(NLOTYPE.EvtPhTyp == 1) {
        p = new GenParticle( CLHEP::HepLorentzVector( 0,0,0,weight ), 99, 1); 
        p->suggest_barcode( 7 );
        prod_vtx->add_particle_out(p);
      } else if(NLOTYPE.EvtPhTyp == 0) {
        p = new GenParticle( CLHEP::HepLorentzVector( 0,0,0,weight ), 98, 1); 
        p->suggest_barcode( 6 );
        prod_vtx->add_particle_out(p);
      }
    }
  }
 
 
  if( log.level() < MSG::INFO ) {
    evt->print();
  }
  
  // Check if the McCollection already exists
  SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
  if (anMcCol!=0) {
    log << MSG::DEBUG << "Adding McGenEvent to existing collection" << endreq;
    McGenEvent* mcEvent = new McGenEvent(evt);
    anMcCol->push_back(mcEvent);
  } else {
    McGenEventCol *mcColl = new McGenEventCol;
    McGenEvent* mcEvent = new McGenEvent(evt);
    mcColl->push_back(mcEvent);
    StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
    if (sc != StatusCode::SUCCESS) {
      log << MSG::ERROR << "Could not register McGenEvent" << endreq;
      delete mcColl;
      delete evt;
      delete mcEvent;
      return StatusCode::FAILURE;
    } else {
      log << MSG::INFO << "McGenEventCol created " << endreq;
    }
  }
 
  m_totalEvents++;
 
  log<<MSG::DEBUG<< "before execute() return"<<endreq;
  return StatusCode::SUCCESS;
}

finalize()

StatusCode Ekhara::finalize

(

)

e+e- --> e+e-pi0pi0

e+e- --> e+e-pi+pi-

Definition at line 352 of file Ekhara.cxx.

                            {
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "Ekhara in finalize()" << endreq;
 
  EKHARA(1);
  
  char xsect[100];
  if(m_finalState == 0) { //! e+e- --> e+e-pi0pi0
    double tnpfinpar[10];
    for(int i=0; i<10; i++) {
      tnpfinpar[i]=0.0;
    }
    
    GET_FINAL_MESON_INFO(0,tnpfinpar);
 
    log << MSG::INFO << "Total number of iterations = " << tnpfinpar[0] << endreq;
    log << MSG::INFO << "Number of MC trials = " << tnpfinpar[1] << endreq;
    log << MSG::INFO << "ContribMax_MCloop  = " << tnpfinpar[2] << endreq;
    log << MSG::INFO << "MC efficiency (evt/eff.iter) = " << m_totalEvents/tnpfinpar[1] << endreq;
 
    log << MSG::INFO << "================================" << endreq;
    log << MSG::INFO << "From all weights:" << endreq;
    double CS  = tnpfinpar[3]/tnpfinpar[0];
    double CSerr = sqrt( fabs(tnpfinpar[4]/tnpfinpar[0]-CS*CS) / (tnpfinpar[0]-1.) );
    log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [nb]" << endreq;
    log << MSG::INFO << endreq;
    log << MSG::INFO << "From UNWEIGHTED events:" << endreq;
    if (tnpfinpar[1] <= 0.) {
      CS = 0.;
      CSerr = 0.;
    } else {
      CS = tnpfinpar[5] * m_totalEvents / tnpfinpar[1];
      CSerr = sqrt(CS*(tnpfinpar[5] - CS)) / sqrt(tnpfinpar[1]);
    }
    log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [nb]" << endreq;
    log << MSG::INFO << "================================" << endreq;
    if (tnpfinpar[7] > 0.) {
      log << MSG::WARNING << "There were overshootings." << endreq;
      log << MSG::WARNING << "The events, corresponding to them were ignored." << endreq;
      log << MSG::WARNING << "The unweighted sample may be not reliable!" << endreq;
      log << MSG::WARNING << "\nContribution from overshooted events:" << endreq;
      sprintf(xsect," Integral Cross Section = %1.5E +- %1.5E  [GeV^-2]",tnpfinpar[8],tnpfinpar[9]);
      log << MSG::WARNING << xsect << endreq;
      sprintf(xsect,"                        = %1.5E +- %1.5E  [nb]",tnpfinpar[8]*tnpfinpar[6],tnpfinpar[9]*tnpfinpar[6]);
      log << MSG::WARNING << xsect << endreq;
      log << MSG::INFO << "================================" << endreq;
    }
  } else if(m_finalState==1) { //! e+e- --> e+e-pi+pi-
    double pipifinpar[2];
    for(int i=0; i<2; i++) {
      pipifinpar[i]=0.0;
    }
    GET_FINAL_TWOPI_INFO(pipifinpar);
 
    double CS = pipifinpar[0]/m_totalEvents;
    double CSerr = sqrt((pipifinpar[1]/m_totalEvents-CS*CS)/(m_totalEvents-1.));
    log << MSG::INFO << "================================" << endreq;
    log << MSG::INFO << "Cross Section = " << CS << " +- " << CSerr << " [nb]" << endreq;
    log << MSG::INFO << "================================" << endreq;
 
 
  } else if(m_finalState>1 && m_finalState<5) {
    double mfp[10];
    for(int i=0; i<10; i++) {
      mfp[i]=0.0;
    }
    GET_FINAL_MESON_INFO(0,mfp);
    if(m_switchNLO) {
      log << MSG::INFO << "================================" << endreq;
      log << MSG::INFO << "Events 0ph:" << endreq;
      log << MSG::INFO << "-------------------------------- " << endreq;
    }
    log << MSG::INFO << "Total number of iterations = " << mfp[0] << endreq;
    if(!m_switchNLO) {
    log << MSG::INFO << "Number of MC trials = " << mfp[1] << endreq;
    log << MSG::INFO << "ContribMax_MCloop  = " << mfp[2] << endreq;
    log << MSG::INFO << "MC efficiency (evt/eff.iter) = " << m_totalEvents/mfp[1] << endreq;
    }
    log << MSG::INFO << "================================" << endreq;
    log << MSG::INFO << "From all weights:" << endreq;
    double CS  = mfp[3]/mfp[0];
    double CSerr = sqrt( fabs(mfp[4]/mfp[0]-CS*CS) / (mfp[0]-1.) );
    double sumCS = CS;
    double sumCSerr = CSerr*CSerr;
    log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [GeV^-2]" << endreq;
    log << MSG::INFO << "                        = " << CS*mfp[6] << " +- " << CSerr*mfp[6] << " [nb]" << endreq;
    log << MSG::INFO << endreq;
    if(m_switchNLO) {
      for(int i=0; i<10; i++) {
        mfp[i]=0.0;
      }
      GET_FINAL_MESON_INFO(1,mfp);
      log << MSG::INFO << "================================" << endreq;
      log << MSG::INFO << "Events 1ph:" << endreq;
      log << MSG::INFO << "-------------------------------- " << endreq;
      log << MSG::INFO << "Total number of iterations = " << mfp[0] << endreq;
      log << MSG::INFO << "================================" << endreq;
      log << MSG::INFO << "From all weights:" << endreq;
      CS  = mfp[3]/mfp[0];
      CSerr = sqrt( fabs(mfp[4]/mfp[0]-CS*CS) / (mfp[0]-1.) );
      sumCS+=CS;
      sumCSerr+=CSerr*CSerr;
      log << MSG::INFO << " Integral Cross Section = " << CS << " +- " << CSerr << " [GeV^-2]" << endreq;
      log << MSG::INFO << "                        = " << CS*mfp[6] << " +- " << CSerr*mfp[6] << " [nb]" << endreq;
      log << MSG::INFO << endreq;
      log << MSG::INFO << " -------------------------------- " << endreq;
      log << MSG::INFO << " sigma (0ph+1ph events) = " << sumCS << " +- " << sqrt(sumCSerr) << " [GeV^-2]" << endreq;
      log << MSG::INFO << "                        = " << sumCS*mfp[6] << " +- " << sqrt(sumCSerr)*mfp[6] << " [nb]" << endreq;
      log << MSG::INFO << "================================" << endreq;
    } else {
      log << MSG::INFO << "From UNWEIGHTED events:" << endreq;
      if (mfp[1] <= 0.) {
        CS = 0.;
        CSerr = 0.;
      } else {
          CS = mfp[5] * m_totalEvents / mfp[1];
        CSerr = sqrt(CS*(mfp[5] - CS)) / sqrt(mfp[1]);
      }
      sprintf(xsect," Integral Cross Section = %1.5E +- %1.5E  [GeV^-2]",CS,CSerr);
      log << MSG::INFO << xsect << endreq;
      sprintf(xsect,"                        = %1.5E +- %1.5E  [nb]",CS*mfp[6],CSerr*mfp[6]);
      log << MSG::INFO << xsect << endreq;
      log << MSG::INFO << "================================" << endreq;
      
      if (mfp[7] > 0.) {
        log << MSG::WARNING << "There were overshootings." << endreq;
        log << MSG::WARNING << "The events, corresponding to them were ignored." << endreq;
        log << MSG::WARNING << "The unweighted sample may be not reliable!" << endreq;
        log << MSG::WARNING << "\nContribution from overshooted events:" << endreq;
        sprintf(xsect," Integral Cross Section = %1.5E +- %1.5E  [GeV^-2]",mfp[8],mfp[9]);
        log << MSG::WARNING << xsect << endreq;
        sprintf(xsect,"                        = %1.5E +- %1.5E  [nb]",mfp[8]*mfp[6],mfp[9]*mfp[6]);
        log << MSG::WARNING << xsect << endreq;
        log << MSG::INFO << "================================" << endreq;
      }
    }
  } else if (m_finalState>4 && m_finalState<8) {
    double chicjfinpar[3];
    for(int i=0; i<3; i++) {
      chicjfinpar[i]=0.0;
    }
    GET_FINAL_CHICJ_INFO(chicjfinpar);
 
    double CS = chicjfinpar[1]/chicjfinpar[0];
    double CSerr = sqrt((chicjfinpar[2]/chicjfinpar[0]-CS*CS)/(chicjfinpar[0]-1.));
    log << MSG::INFO << "================================" << endreq;
    sprintf(xsect," Cross Section = %1.5E +- %1.5E  [nb]",CS,CSerr);
    log << MSG::INFO << xsect << endreq;
    log << MSG::INFO << "================================" << endreq;
 
 
 
  }
 
  log << MSG::DEBUG << "Ekhara has terminated successfully" << endreq;
 
  return StatusCode::SUCCESS;
}

initialize()

StatusCode Ekhara::initialize

(

)

• Set Bes unified random engine

• Perform consitency checks of jobOptions and print screen messages

• Transfer integer parameters to FORTRAN common blocks

• Prepare floating point parameters to be transfered to quadruple precision

• Transfer floating point parameters to FORTAN routines
• and initialize FORTAN routines

• Initialize event counter

Definition at line 85 of file Ekhara.cxx.

                              {
  MsgStream log(msgSvc(), name());
  log << MSG::DEBUG << "Ekhara in initialize()" << endreq;
 
  //!- Set Bes unified random engine
  static const bool CREATEIFNOTTHERE(true);
  StatusCode RndmStatus = service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
  if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc) {
    log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
    return RndmStatus;
  }
  HepRandomEngine* engine  = p_BesRndmGenSvc->GetEngine("Ekhara");
  EkharaRandom::setRandomEngine(engine);
 
 
  //!- Perform consitency checks of jobOptions and print screen messages
  if(!checkAndPrintParameters()) return StatusCode::FAILURE;
 
 
  //!- Transfer integer parameters to FORTRAN common blocks
  CHANNELSEL.channel_id = m_finalState;
  CHANNELSEL.NLO1P = 0;
  if( m_switchNLO ) CHANNELSEL.NLO1P = 1;
  CHANNELSEL.VPSW = 0;
  if( m_switchVP ) CHANNELSEL.VPSW = 1;
  CHANNELSEL.phsp_2pi = m_twoPionPhspAlg;
  SWDIAG.sw_2pi = m_twoPionAmplitudes;
  SWDIAG.sw_1pi = m_mesonAmplitudes;
  PIONFFSW.piggFFsw = m_TFF_ID;
  TAGGINGMODE.tagmode = m_taggingMode;
  FFPARAMSET.FF_pion = m_twoPionFormFactor;
  NLOTYPE.weighted = 0;
  if( m_nloWithWeights) NLOTYPE.weighted = 1;
 
  //!- Prepare floating point parameters to be transfered to quadruple precision
  double xpar[27];
  xpar[0]  = m_cmsEnergy;
  xpar[1]  = m_numUpperBound;
  xpar[2]  = m_ignoreUpperBound;
  xpar[3]  = m_upperBoundEnlarge;
  xpar[4]  = m_eps_ph;
  xpar[5]  = m_posiThetaMin;
  xpar[6]  = m_posiThetaMax;
  xpar[7]  = m_posiEnergyMin;
  xpar[8]  = m_posiEnergyMax;
  xpar[9]  = m_elecThetaMin;
  xpar[10] = m_elecThetaMax;
  xpar[11] = m_elecEnergyMin;
  xpar[12] = m_elecEnergyMax;
  xpar[13] = m_mesonThetaMin;
  xpar[14] = m_mesonThetaMax;
  xpar[15] = m_mesonEnergyMin;
  xpar[16] = m_mesonEnergyMax;
  xpar[17] = m_twoPionThetaMin;
  xpar[18] = m_twoPionThetaMax;
  xpar[19] = m_twoPionMissThetaMin;
  xpar[20] = m_twoPionMissThetaMax;
  xpar[21] = m_chicjThetaMin;
  xpar[22] = m_chicjThetaMax;
  xpar[23] = m_chicjEnergyMin;
  xpar[24] = m_chicjEnergyMax;
  xpar[25] = m_taggingAngle;
  xpar[26] = m_taggingQsquare;
 
  //!- Transfer floating point parameters to FORTAN routines
  //!- and initialize FORTAN routines
  BOSS_INIT_EKHARA(xpar);
 
  if( log.level() < MSG::INFO ) {
    DIAGNOSE();
  }
 
  //!- Initialize event counter
  //should be replaced by request to ApplicationManager for EvtMax
  m_totalEvents = 0;
  
  return StatusCode::SUCCESS;
}

---

The documentation for this class was generated from the following files:

• Ekhara.h
• Ekhara.cxx