CosmicGenerator Class Reference

#include <CosmicGenerator.h>

Inheritance diagram for CosmicGenerator:

Public Member Functions
	CosmicGenerator (const std::string &name, ISvcLocator *pSvcLocator)
	~CosmicGenerator ()
StatusCode	initialize ()
StatusCode	execute ()
StatusCode	finalize ()
HepLorentzVector	generateVertex (void)

Static Public Attributes
static IBesRndmGenSvc *	p_BesRndmGenSvc = 0

Detailed Description

Definition at line 45 of file CosmicGenerator.h.

Constructor & Destructor Documentation

CosmicGenerator()

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

Definition at line 143 of file CosmicGenerator.cxx.

    : Algorithm( name, pSvcLocator )
//--------------------------------------------------------------------------
{
  //
  // Migration to MeV and mm units: all conversions are done in this interface
  // to the CosmicGun. The CosmicGun itself uses GeV units internally - to call
  // the fortran code.
  //
 
  m_GeV      = 1000;
  m_mm       = 10;
  m_readfile = false;
  m_exzCut   = false;
  m_tuple    = 0;
  m_tuple1   = 0;
  declareProperty( "eventfile", m_infile = "NONE" );
  declareProperty( "emin", m_emin = 0.1 );
  declareProperty( "emax", m_emax = 10. );
  declareProperty( "xvert_low", m_xlow = -110.7 ); // EMC BSCRmax2
  declareProperty( "xvert_hig", m_xhig = 110.7 );
  declareProperty( "zvert_low", m_zlow = -171.2 ); // EMC WorldZposition+0.5*CrystalLength
  declareProperty( "zvert_hig", m_zhig = 171.2 );
  declareProperty( "yvert_val", m_yval = 0 );
  declareProperty( "tmin", m_tmin = 0. );
  declareProperty( "tmax", m_tmax = 24. );
  declareProperty( "tcor", m_tcor = 0. );
  declareProperty( "IPx", m_IPx = 0. );
  declareProperty( "IPy", m_IPy = 0. );
  declareProperty( "IPz", m_IPz = 0. );
  declareProperty( "Radius", m_radius = 0. );
  declareProperty( "ExzCut", m_exzCut );
  declareProperty( "CubeProjection", m_cubeProj = true );
  declareProperty( "OptimizeSphere", m_sphereOpt = false );
 
  declareProperty( "SwapYZAxis", m_swapYZAxis = false );
  declareProperty( "ctcut", m_ctcut = 0.0 ); // according to sin(acos(0.93)) =0.368
  // declareProperty("ReadTTR",      m_readTTR =false );
  //  declareProperty("ReadTR",      m_readTTR =false );
  declareProperty( "PrintEvent", m_printEvent = 10 );
  declareProperty( "PrintMod", m_printMod = 100 );
  declareProperty( "RMax", m_rmax = 10000000. );
  declareProperty( "ThetaMin", m_thetamin = 0. );
  declareProperty( "ThetaMax", m_thetamax = 1. );
  declareProperty( "PhiMin", m_phimin = -1 * PI );
  declareProperty( "PhiMax", m_phimax = PI );
  declareProperty( "Xposition", m_xpos = 263.5 - 0.0001 ); // 263.5*cm,263.5*cm,287.5*cm
  declareProperty( "Yposition", m_ypos = 263.5 - 0.0001 );
  declareProperty( "Zposition", m_zpos = 287.5 - 0.0001 );
 
  declareProperty( "DoublePlaneTrigger", m_doublePlaneTrigger = false );
  declareProperty( "xMinTop", m_xposMin_top = -16 );       // in cm
  declareProperty( "xMaxTop", m_xposMax_top = 16 );        // in cm
  declareProperty( "yTop", m_ypos_top = 52 );              // in cm
  declareProperty( "zMinTop", m_zposMin_top = -25 );       // in cm
  declareProperty( "zMaxTop", m_zposMax_top = 25 );        // in cm
  declareProperty( "xMinBottom", m_xposMin_bottom = -8 );  // in cm
  declareProperty( "xMaxBottom", m_xposMax_bottom = 8 );   // in cm
  declareProperty( "yBottom", m_ypos_bottom = -26 );       // in cm
  declareProperty( "zMinBottom", m_zposMin_bottom = -25 ); // in cm
  declareProperty( "zMaxBottom", m_zposMax_bottom = 25 );  // in cm
 
  declareProperty( "DumpMC", m_dumpMC = 0 );
}

Referenced by CosmicGenerator().

~CosmicGenerator()

CosmicGenerator::~CosmicGenerator

(

)

Definition at line 209 of file CosmicGenerator.cxx.

{}

Member Function Documentation

execute()

StatusCode CosmicGenerator::execute

(

)

Definition at line 387 of file CosmicGenerator.cxx.

                                    {
  //---------------------------------------------------------------------------
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "CosmicGenerator executing" << endmsg;
 
  ++m_events;
  //  MsgStream log(msgSvc(), name());
  log << MSG::DEBUG << "Event #" << m_events << endmsg;
 
  // clear up the vectors
  m_fourPos.clear();
  m_fourMom.clear();
  m_polarization.clear();
  m_pdgCode.clear();
 
  if ( m_readfile )
  {
    if ( !m_ffile.eof() )
    {
      CosmicEventParser evt;
      m_ffile >> evt;
 
      std::cout << evt << std::endl;
 
      double       polx = 0;
      double       poly = 0;
      double       polz = 0;
      Polarization thePolarization;
      thePolarization.set_normal3d( HepNormal3D( polx, poly, polz ) );
      m_polarization.push_back( thePolarization );
 
      //
      // units are already converted to MeV's and mm.
      //
      m_fourPos.push_back( evt.Vertex() );
      m_fourMom.push_back( evt.Momentum() );
      m_pdgCode.push_back( evt.pdgID() );
    }
    else
    {
      log << MSG::FATAL << "End of file reached - stop " << endmsg;
      exit( 1 );
      return StatusCode::FAILURE;
    }
  }
 
  else
  {
 
    bool             accepted  = false;
    bool             projected = false;
    HepLorentzVector pp;
    CosmicGun*       gun = CosmicGun::GetCosmicGun();
    HepLorentzVector vert;
    HepLorentzVector vert_proj;
    while ( !accepted )
    {
      m_tried++;
      vert = generateVertex();
      Hep3Vector vert3( vert.x(), vert.y(), vert.z() );
 
      pp = gun->GenerateEvent();
      if ( m_dumpMC == 1 )
      {
        m_cosmicE      = pp.e() * m_GeV;
        m_cosmicTheta  = pp.theta();
        m_cosmicPhi    = pp.phi();
        m_cosmicCharge = gun->GetMuonCharge();
        m_tuple->write();
      }
      double theta1 = pp.theta();
      double phi1   = pp.phi();
      double mag1   = pp.vect().mag();
 
      Hep3Vector pp_corr( mag1 * sin( theta1 ) * cos( phi1 ), -mag1 * cos( theta1 ),
                          mag1 * sin( theta1 ) * sin( phi1 ) );
      Hep3Vector direction( pp_corr.x(), pp_corr.y(), pp_corr.z() );
      Hep3Vector center_dir = m_center - vert3;
      // if optimization activated, check for the direction of the generated muon
 
      if ( m_cubeProj )
      {
        double alpha, beta;
        if ( m_sphereOpt )
        {
 
          beta  = direction.angle( center_dir );
          alpha = asin( m_radius / center_dir.mag() );
          if ( fabs( beta ) > alpha )
          {
            log << MSG::DEBUG << alpha << ", " << beta << " rejected muon due to sphere cut! "
                << endmsg;
            m_rejected++;
            continue;
          }
        }
 
        double     l_fight0, l_fight1, l_fight2;
        double     coor_x0, coor_y0, coor_z0;
        double     coor_x1, coor_y1, coor_z1;
        double     coor_x2, coor_y2, coor_z2;
        bool       isIn0 = false;
        bool       isIn1 = false;
        bool       isIn2 = false;
        HepPoint3D vert_pro0;
        HepPoint3D vert_pro1;
        HepPoint3D vert_pro2;
        HepPoint3D vert_org( vert3 );
 
        coor_y0 = m_ypos; // defined in jobOpt.
        coor_x0 = direction.x() * ( coor_y0 - vert.y() ) / direction.y() + vert.x();
        coor_z0 = direction.z() * ( coor_y0 - vert.y() ) / direction.y() + vert.z();
 
        if ( fabs( coor_x0 ) <= m_xpos && fabs( coor_z0 ) <= m_zpos )
        {
          vert_pro0 = HepPoint3D( coor_x0, coor_y0, coor_z0 );
          l_fight0  = vert_org.distance( vert_pro0 );
          isIn0     = true;
        }
        else
        {
 
          coor_z1 = sign( coor_z0 - vert.z() ) * m_zpos; // defined in jobOpt.
          coor_x1 = direction.x() * ( coor_z1 - vert.z() ) / direction.z() + vert.x();
          coor_y1 = direction.y() * ( coor_z1 - vert.z() ) / direction.z() + vert.y();
 
          vert_pro1 = HepPoint3D( coor_x1, coor_y1, coor_z1 );
          l_fight1  = vert_org.distance( vert_pro1 );
 
          coor_x2   = sign( coor_x0 - vert.x() ) * m_xpos; // defined in jobOpt.
          coor_z2   = direction.z() * ( coor_x2 - vert.x() ) / direction.x() + vert.z();
          coor_y2   = direction.y() * ( coor_x2 - vert.x() ) / direction.x() + vert.y();
          vert_pro2 = HepPoint3D( coor_x2, coor_y2, coor_z2 );
          l_fight2  = vert_org.distance( vert_pro2 );
          if ( l_fight1 < l_fight2 ) isIn1 = true;
          else isIn2 = true;
        }
        // reset vert(x,y,z,t), after projection
        if ( isIn0 )
        {
          vert_proj =
              HepLorentzVector( vert_pro0.x(), vert_pro0.y(), vert_pro0.z(), l_fight0 );
 
          projected = true;
          accepted  = true;
          m_accepted++;
        }
        else if ( isIn1 )
        {
          vert_proj =
              HepLorentzVector( vert_pro1.x(), vert_pro1.y(), vert_pro1.z(), l_fight1 );
 
          projected = true;
          accepted  = true;
          m_accepted++;
        }
        else if ( isIn2 )
        {
          vert_proj =
              HepLorentzVector( vert_pro2.x(), vert_pro2.y(), vert_pro2.z(), l_fight2 );
 
          projected = true;
          accepted  = true;
          m_accepted++;
        }
        else
        {
          log << MSG::DEBUG << " rejected muon due to cube projection request! " << endmsg;
          m_rejected++;
        }
      }
      else if ( m_doublePlaneTrigger )
      {
        double coor_x0, coor_y0, coor_z0; // intersection with the top plane
        coor_y0 = m_ypos_top;
        coor_x0 = direction.x() * ( coor_y0 - vert.y() ) / direction.y() + vert.x();
        coor_z0 = direction.z() * ( coor_y0 - vert.y() ) / direction.y() + vert.z();
 
        double coor_x1, coor_y1, coor_z1; // intersection with the bottom plane
        coor_y1 = m_ypos_bottom;
        coor_x1 = direction.x() * ( coor_y1 - vert.y() ) / direction.y() + vert.x();
        coor_z1 = direction.z() * ( coor_y1 - vert.y() ) / direction.y() + vert.z();
 
        if ( coor_x0 > m_xposMin_top && coor_x0 < m_xposMax_top && coor_z0 > m_zposMin_top &&
             coor_z0 < m_zposMax_top && coor_x1 > m_xposMin_bottom &&
             coor_x1 < m_xposMax_bottom && coor_z1 > m_zposMin_bottom &&
             coor_z1 < m_zposMax_bottom )
        {
          accepted = true;
          m_accepted++;
        }
        else m_rejected++;
      }
      else accepted = true; // if no opt required accept the first muon
    }
 
    //      pp.setX(pp.x()*m_GeV);
    //      pp.setY(pp.y()*m_GeV);
    //      pp.setZ(pp.z()*m_GeV);
    //      pp.setT(pp.t()*m_GeV);
 
    pp.setX( pp.x() );
    pp.setY( pp.y() );
    pp.setZ( pp.z() );
    pp.setT( pp.t() );
    // Get the mass of the particle to be generated
    int charge = gun->GetMuonCharge();
    // m_pdgCode.push_back(charge*13);
    m_pdgCode.push_back( charge * -13 );
 
    //      const ParticleData* pdtparticle =
    //      m_particleTable->particle(ParticleID(abs(m_pdgCode.back()))); double mass =
    //      pdtparticle->mass().value(); double mass =105.5658;
 
    // Compute the kinematic values.  First, the vertex 4-vector:
 
    // Set the polarization.  Realistically, this is going to be zero
    // for most studies, but you never know...
    double polx = 0;
    double poly = 0;
    double polz = 0;
    // m_polarization.set_normal3d(HepNormal3D(polx,poly,polz));
    Polarization thePolarization;
 
    // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ?
    // if not...do nothing...if so, invert position of y- and z- coordinate
    //
    // well and don't forget about the direction of the incoming cosmic muon(s) either
    // that means:  y -> -y
    //
    if ( !m_swapYZAxis )
    {
      // thePolarization.set_normal3d(HepNormal3D(polx,-poly,polz));
      thePolarization.set_normal3d( HepNormal3D( polx, poly, polz ) );
    }
    else thePolarization.set_normal3d( HepNormal3D( polx, polz, -poly ) );
 
    m_polarization.push_back( thePolarization );
 
    // The method of calculating e, theta, and phi depends on the user's
    // commands.  Let the KinematicManager handle it.
    double e     = pp.e();
    double theta = pp.theta();
    double phi   = pp.phi();
 
    // At this point, we have e, theta, and phi.  Put them together to
    // get the four-momentum.
 
    double p2 = e * e - mass * mass;
    if ( p2 < 0 )
    {
      log << MSG::ERROR << "Event #" << m_events << " E=" << e << ", mass=" << mass
          << " -- you have generated a tachyon! Increase energy or change particle ID."
          << endmsg;
      return StatusCode::FAILURE;
    }
 
    double p  = sqrt( p2 );
    double px = p * sin( theta ) * cos( phi );
    double pz = p * sin( theta ) * sin( phi );
    double py = -p * cos( theta );
 
    // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ?
    // if not...do nothing...if so, invert position of y- and z- coordinate
    //
    // well and don't forget about the direction of the incoming cosmic muon(s) either
    // that means:  y -> -y
    //
    if ( !m_swapYZAxis )
    {
      // Line below corrupted py sign and forces muons to be upwards, not downwards.
      // m_fourMom.push_back(HepLorentzVector(px,-py,pz,pp.e()));
      m_fourMom.push_back( HepLorentzVector( px, py, pz, pp.e() ) );
    }
    else m_fourMom.push_back( HepLorentzVector( px, pz, -py, pp.e() ) );
    double x = vert.x();
    double y = vert.y();
    double z = vert.z();
    double t = vert.t();
    if ( projected )
    {
 
      x = vert_proj.x();
      y = vert_proj.y();
      z = vert_proj.z();
      t = vert.t() - vert_proj.t() / HepLorentzVector( px, py, pz, pp.e() ).beta() + m_tcor;
    }
    //      log << MSG::DEBUG
    //          << "proj:"<<m_cubeProj<<" ("<<x<<", "<<y<<", "<<z<<", "<<t<<")"
    //      <<endmsg;
    // Do we need to swap Y- and Z-axis for the PixelEndCap A Cosmic test ?
    // if not...do nothing...if so, invert position of y- and z- coordinate
    //
    // but not only that...change also the direction of the incoming cosmic muon(s),
    // they must go towards the pixel endcap A, i.e. y -> -y
    //
 
    if ( !m_swapYZAxis ) m_fourPos.push_back( HepLorentzVector( x, y, z, t ) );
    else m_fourPos.push_back( HepLorentzVector( x, z, y, t ) );
 
    log << MSG::DEBUG << "  (x,y,z,t) = (" << m_fourPos.back().x() << ","
        << m_fourPos.back().y() << "," << m_fourPos.back().z() << "," << m_fourPos.back().t()
        << "), (Px,Py,Pz,E) = (" << m_fourMom.back().px() << "," << m_fourMom.back().py()
        << "," << m_fourMom.back().pz() << "," << m_fourMom.back().e() << ")" << endmsg;
    log << MSG::DEBUG << "  (theta,phi) = (" << theta << "," << phi << "), "
        << "polarization(x,y,z) = (" << m_polarization.back().normal3d().x() << ","
        << m_polarization.back().normal3d().y() << "," << m_polarization.back().normal3d().z()
        << ")" << endmsg;
    if ( m_dumpMC == 1 )
    {
      //      m_cosmicE=pp.e()*m_GeV;
      //      m_cosmicTheta=pp.theta();
      //      m_cosmicPhi=pp.phi();
      //      m_cosmicCharge=gun->GetMuonCharge();
      mc_theta = m_fourMom.back().theta();
      mc_phi   = m_fourMom.back().phi();
      mc_px    = m_fourMom.back().px();
      mc_py    = m_fourMom.back().py();
      mc_pz    = m_fourMom.back().pz();
      //             m_cosmicE=pz;
      m_tuple1->write();
    }
  }
 
  // fill evt
 
  GenEvent* event = new GenEvent( 1, 1 );
 
  if ( m_fourMom.size() == m_fourPos.size() && m_fourMom.size() == m_polarization.size() )
  {
 
    for ( std::size_t v = 0; v < m_fourMom.size(); ++v )
    {
 
      // Note: The vertex and particle are owned by the event, so the
      // event is responsible for those pointers.
 
      // Create the particle, and specify its polarization.
 
      GenParticle* particle = new GenParticle( m_fourMom[v], m_pdgCode[v], 1 );
      particle->set_polarization( m_polarization[v] );
 
      // Create the vertex, and add the particle to the vertex.
      GenVertex* vertex = new GenVertex( m_fourPos[v] );
      vertex->add_particle_out( particle );
 
      // Add the vertex to the event.
      event->add_vertex( vertex );
    }
 
    event->set_event_number( m_events ); // Set the event number
  }
  else
  {
 
    log << MSG::ERROR << "Wrong different number of vertexes/momenta/polaritazions!" << endmsg;
    return StatusCode::FAILURE;
  }
  // Check if the McCollection already exists
  SmartDataPtr<McGenEventCol> anMcCol( eventSvc(), "/Event/Gen" );
  if ( anMcCol != 0 )
  {
    // Add event to existing collection
 
    log << MSG::INFO << "Add McGenEvent to existing collection" << endmsg;
    McGenEvent* mcEvent = new McGenEvent( event );
    anMcCol->push_back( mcEvent );
  }
  else
  {
    // Create Collection and add  to the transient store
    McGenEventCol* mcColl  = new McGenEventCol;
    McGenEvent*    mcEvent = new McGenEvent( event );
    mcColl->push_back( mcEvent );
    StatusCode sc = eventSvc()->registerObject( "/Event/Gen", mcColl );
    if ( sc != StatusCode::SUCCESS )
    {
      log << MSG::ERROR << "Could not register McGenEvent" << endmsg;
      delete mcColl;
      delete event;
      delete mcEvent;
      return StatusCode::FAILURE;
    }
  }
 
  return StatusCode::SUCCESS;
}

finalize()

StatusCode CosmicGenerator::finalize

(

)

Definition at line 776 of file CosmicGenerator.cxx.

                                     {
  //---------------------------------------------------------------------------
  // Get the KinematicManager.
 
  if ( m_cubeProj )
  {
 
    MsgStream log( msgSvc(), name() );
 
    log << MSG::INFO << "***************************************************" << endmsg;
    log << MSG::INFO << "** you have run CosmicGenerator with some " << endmsg;
    log << MSG::INFO << "** filters for cosmic muon simulation" << endmsg;
    log << MSG::INFO << "** " << m_tried << " muons were tried" << endmsg;
    log << MSG::INFO << "** " << m_accepted << " muons were accepted" << endmsg;
    log << MSG::INFO << "** " << m_rejected << " muons were rejected" << endmsg;
    log << MSG::INFO << "***************************************************" << endmsg;
    std::cout << "***************************************************" << std::endl;
    std::cout << "** you have run CosmicGenerator with some " << std::endl;
    std::cout << "** filters for cosmic muon simulation" << std::endl;
    std::cout << "** " << m_tried << " muons were tried" << std::endl;
    std::cout << "** " << m_accepted << " muons were accepted" << std::endl;
    std::cout << "** " << m_rejected << " muons were rejected" << std::endl;
    std::cout << "***************************************************" << std::endl;
  }
 
  return StatusCode::SUCCESS;
}

generateVertex()

HepLorentzVector CosmicGenerator::generateVertex

(

void

)

Definition at line 361 of file CosmicGenerator.cxx.

                                                       {
 
  MsgStream log( msgSvc(), name() );
 
  // Get the pointer to the engine of the stream named SINGLE. If the
  // stream does not exist is created automaticaly
  //  HepRandomEngine* engine = CosmicGenerator::p_AtRndmGenSvc->GetEngine("COSMICS");
  HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine( "PYTHIA" );
  // Generate a random number according to the distribution.
 
  float x_val = RandFlat::shoot( engine, m_xlow, m_xhig );
  float z_val = RandFlat::shoot( engine, m_zlow, m_zhig );
 
  // Generate a random number for time offset
  float t_val = 0;
  if ( m_tmin < m_tmax ) { t_val = RandFlat::shoot( engine, m_tmin, m_tmax ); }
  else if ( m_tmin == m_tmax ) { t_val = m_tmin; }
  else
    log << MSG::FATAL << " You specified m_tmin = " << m_tmin << " and m_tmax " << m_tmax
        << endmsg;
  HepLorentzVector p( x_val, m_yval, z_val, t_val * c_light );
 
  return p;
}

Referenced by execute().

initialize()

StatusCode CosmicGenerator::initialize

(

)

Definition at line 214 of file CosmicGenerator.cxx.

                                       {
  //---------------------------------------------------------------------------
 
  // Initialize event count.
 
  m_events   = 0;
  m_tried    = 0;
  m_accepted = 0;
  m_rejected = 0;
  if ( m_emin < 0.1 )
  {
    m_emin = 0.1;
    std::cout << "WARNING: set emin to 0.1 GeV" << std::endl;
  }
  m_emin   = m_emin * m_GeV;
  m_emax   = m_emax * m_GeV;
  m_xlow   = m_xlow * m_mm;
  m_xhig   = m_xhig * m_mm;
  m_zlow   = m_zlow * m_mm;
  m_zhig   = m_zhig * m_mm;
  m_yval   = m_yval * m_mm;
  m_xpos   = m_xpos * m_mm;
  m_ypos   = m_ypos * m_mm;
  m_zpos   = m_zpos * m_mm;
  m_radius = m_radius * m_mm;
  m_tcor   = m_tcor + sqrt( ( m_xpos - m_xlow ) * ( m_xpos - m_xlow ) +
                            ( m_zpos - m_zlow ) * ( m_zpos - m_zlow ) +
                            ( m_ypos - m_yval ) * ( m_ypos - m_yval ) ) /
                        ( m_emin / sqrt( m_emin * m_emin + mass * mass * m_GeV * m_GeV ) );
 
  m_xposMin_top *= m_mm;
  m_xposMax_top *= m_mm;
  m_ypos_top *= m_mm;
  m_zposMin_top *= m_mm;
  m_zposMax_top *= m_mm;
  m_xposMin_bottom *= m_mm;
  m_xposMax_bottom *= m_mm;
  m_ypos_bottom *= m_mm;
  m_zposMin_bottom *= m_mm;
  m_zposMax_bottom *= m_mm;
 
  ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap.h
  p_msgSvc                = 0;
  StatusCode result       = svcLocator->service( "MessageSvc", p_msgSvc );
 
  if ( !result.isSuccess() || p_msgSvc == 0 )
  {
    std::cerr << "VGenCommand(): could not initialize MessageSvc!" << std::endl;
    throw GaudiException( "Could not initalize MessageSvc", "CosmicGenerator",
                          StatusCode::FAILURE );
  }
 
  MsgStream log( p_msgSvc, "CosmicGenerator" );
  if ( m_dumpMC == 1 )
  {
    StatusCode status;
    NTuplePtr  nt( ntupleSvc(), "FILE1/comp" );
    if ( nt ) m_tuple = nt;
    else
    {
      m_tuple = ntupleSvc()->book( "FILE1/comp", CLID_ColumnWiseTuple, "ks N-Tuple example" );
      if ( m_tuple )
      {
        status = m_tuple->addItem( "cosmic_e", m_cosmicE );
        status = m_tuple->addItem( "cosmic_theta", m_cosmicTheta );
        status = m_tuple->addItem( "cosmic_phi", m_cosmicPhi );
        status = m_tuple->addItem( "cosmic_charge", m_cosmicCharge );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
    NTuplePtr nt1( ntupleSvc(), "FILE1/compp" );
    if ( nt1 ) m_tuple1 = nt1;
    else
    {
      m_tuple1 =
          ntupleSvc()->book( "FILE1/compp", CLID_ColumnWiseTuple, "ks N-Tuple example" );
      if ( m_tuple1 )
      {
        status = m_tuple1->addItem( "mc_theta", mc_theta );
        status = m_tuple1->addItem( "mc_phi", mc_phi );
        status = m_tuple1->addItem( "mc_px", mc_px );
        status = m_tuple1->addItem( "mc_py", mc_py );
        status = m_tuple1->addItem( "mc_pz", mc_pz );
      }
      else
      {
        log << MSG::ERROR << "    Cannot book N-tuple:" << long( m_tuple1 ) << endmsg;
        return StatusCode::FAILURE;
      }
    }
  }
  if ( m_infile == "NONE" )
 
  {
    // Initialize the BES random number services.
 
    static const bool CREATEIFNOTTHERE( true );
    StatusCode        RndmStatus =
        svcLocator->service( "BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE );
 
    if ( !RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc )
    {
      log << MSG::ERROR << " Could not initialize Random Number Service" << endmsg;
      return RndmStatus;
    }
 
    CosmicGun* gun = CosmicGun::GetCosmicGun();
 
    gun->SetEnergyRange( m_emin / m_GeV, m_emax / m_GeV );
    gun->SetCosCut( m_ctcut );
    gun->PrintLevel( m_printEvent, m_printMod );
    float flux_withCT = gun->InitializeGenerator();
 
    log << MSG::INFO << "Initialisation cosmic gun done." << endmsg;
    log << MSG::INFO << "Accepted diff  flux after E and cos(theta) cuts = " << flux_withCT
        << " /cm^2/s" << endmsg;
    log << MSG::INFO << "Accepted total flux after E and cos(theta) cuts = "
        << flux_withCT * ( m_xhig - m_xlow ) / m_mm * ( m_zhig - m_zlow ) / m_mm << " /s"
        << endmsg;
 
    std::cout << "Accepted diff  flux after E and cos(theta) cuts = " << flux_withCT
              << " /cm^2/s" << std::endl;
    std::cout << "Accepted total flux after E and cos(theta) cuts = "
              << flux_withCT * ( m_xhig - m_xlow ) / m_mm * ( m_zhig - m_zlow ) / m_mm << " /s"
              << std::endl;
  }
  else
  {
    log << MSG::INFO << "Cosmics are read from file " << m_infile << endmsg;
    m_ffile.open( m_infile.c_str() );
    if ( !m_ffile )
    {
      log << MSG::FATAL << "Could not open input file - stop! " << endmsg;
      return StatusCode::FAILURE;
    }
    m_readfile = true;
  }
 
  m_center = Hep3Vector( m_IPx, m_IPy, m_IPz );
  log << MSG::INFO << "Initialisation done" << endmsg;
  return StatusCode::SUCCESS;
}

Member Data Documentation

p_BesRndmGenSvc

IBesRndmGenSvc * CosmicGenerator::p_BesRndmGenSvc = 0

static

Definition at line 55 of file CosmicGenerator.h.

Referenced by cosmicrndm_(), generateVertex(), and initialize().

---

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

• CosmicGenerator.h
• CosmicGenerator.cxx