CosmicGenerator.cxx

Go to the documentation of this file.

// -------------------------------------------------------------
// File: CosmicGenerator/CosmicGenerator.cxx
// Description:
//
// Non-simulation package CosmicGenerator
//
// A Generator providing cosmicmuon flux at ground level (X-Z plane)
//
// Generator based on fits to measured surface fluxes, used by ALEPH
//
// Ported From Atlas to BES
//                   by Beijiang Liu <liubj@mail.ihep.ac.cn>
//
// Modification to generate on a plane in the center of detector, then project back to Geant 4
// world
//                      Beijiang Liu 18-May-2008
//
// Modification for cleanup
//                      Beijiang Liu 28-Nov-2007
// /
// Modification for  increasing efficiency of muon hitting the detector:
//            m_sphereOpt, Point to a sphere around the interaction point
//                      Beijiang Liu 28-Nov-2007
//
// Modification to project the muons to the Geant4 world
//                      Beijiang Liu 18-Nov-2007
//
//
//
//
// -------------------------------------------------------------
// File: CosmicGenerator/CosmicGenerator.cxx
// Description:
 
//    The output will be stored in the transient event store so it can be
//    passed to the simulation.
//
// AuthorList:
//         W. Seligman: Initial Code 08-Nov-2002,
//         based on work by M. Shapiro and I. Hinchliffe
//
 
// Modification for increasing efficiency of muon hitting the detector:
//                     H. Ma.    March 17, 2006
//   Property: ExzCut:
//  if true, the method exzCut(...) will be called to apply a
//               energy dependent position cut on the surface.
//               This rejects low energy muons at large distance.
//   Property: RMax
//               Used by exzCut to reject non-projective muons, which are
//               too far out on the surface
 
// Modifications to accomodate Pixel EndCap C Cosmic Test needs
//      Marian Zdrazil   June 7, 2006   mzdrazil@lbl.gov
//
// Modifications to accomodate replacement of Pixel EndCap C by a Pixel EndCap A
//      Marian Zdrazil   November 24, 2006  mzdrazil@lbl.gov
//
//
// Description:
// ------------
// It is easier and actually more useful to leave the EndCap A
// in the vertical position (the way it is positioned in the ATLAS detector)
// instead of rotating it clockwise by 90deg which corresponds to the
// placement during the Pixel EndCap A cosmic test in SR1 in November 2006.
// This is why we will generate cosmic muons coming from the positive Z-axis
// direction better than rotating the whole setup in PixelGeoModel.
//
// Modifications July 3rd 2007, Rob McPherson
//     - Fix mu+/mu- bug (always present in Athena versions)
//     - Fix sign of Py (since tag CosmicGenerator-00-00-21, muons only upward-going)
 
#include "CosmicGenerator.h"
#include "CosmicEventParser.h"
#include "CosmicGun.h"
 
// Framework Related Headers:-
#include "GaudiKernel/MsgStream.h"
/// #include "StoreGate/StoreGateSvc.h"
 
// Other classes used by this class:-
#include "CLHEP/Vector/ThreeVector.h"
// #include "HepPDT/ParticleDataTable.hh"
 
////#include "CLHEP/HepPDT/ParticleData.hh"
#include "CLHEP/Units/PhysicalConstants.h"
 
#include "HepMC/GenEvent.h"
#include "HepMC/GenParticle.h"
#include "HepMC/GenVertex.h"
#include "HepMC/Polarization.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/GaudiException.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
 
#include "GaudiKernel/IDataProviderSvc.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/NTuple.h"
 
#include "GeneratorObject/McGenEvent.h"
 
// #include "TrackRecord/TimedTrackRecordCollection.h"
/// #include "TrackRecord/TrackRecordCollection.h"
///  For the Athena-based random numbers.
/// #include "CLHEP/Random/RandFlat.h"
// For BES random numnber service
#include "BesRndmGenSvc/IBesRndmGenSvc.h"
 
#include "CLHEP/Random/RandFlat.h"
#include "CLHEP/Random/Ranlux64Engine.h"
 
#include <cmath>
#include <fstream>
#include <string>
#include <vector>
using namespace CLHEP;
 
static inline int sign( double x ) { return ( x > 0 ? 1 : -1 ); }
 
// definition of PI
float  PI   = 3.1415927; // <== Should use M_PI (defined in standard cmath header)
double mass = 0.1055658;
 
//// Pointer On AtRndmGenSvc
// IAtRndmGenSvc*   CosmicGenerator::p_AtRndmGenSvc = 0;
IBesRndmGenSvc* CosmicGenerator::p_BesRndmGenSvc = 0;
 
extern "C" float cosmicrndm_( int* /*dummy*/ ) {
  //  HepRandomEngine* engine = CosmicGenerator::p_AtRndmGenSvc->GetEngine("COSMICS");
  HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine( "PYTHIA" );
  //  std::cout<<"seed: "<<engine->getSeed()<<", "<< RandFlat::shoot(engine);
  return RandFlat::shoot( engine );
}
 
DECLARE_COMPONENT( CosmicGenerator )
//--------------------------------------------------------------------------
CosmicGenerator::CosmicGenerator( const std::string& name, ISvcLocator* pSvcLocator )
    : 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 );
}
 
//--------------------------------------------------------------------------
CosmicGenerator::~CosmicGenerator()
//--------------------------------------------------------------------------
{}
 
//---------------------------------------------------------------------------
StatusCode CosmicGenerator::initialize() {
  //---------------------------------------------------------------------------
 
  // 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;
}
 
HepLorentzVector CosmicGenerator::generateVertex( void ) {
 
  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;
}
 
//---------------------------------------------------------------------------
StatusCode CosmicGenerator::execute() {
  //---------------------------------------------------------------------------
  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;
}
 
//---------------------------------------------------------------------------
StatusCode CosmicGenerator::finalize() {
  //---------------------------------------------------------------------------
  // 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;
}
 
// Energy dependent position cut on the surface.
bool CosmicGenerator::exzCut( const Hep3Vector& pos, const HepLorentzVector& p ) {
  // p is in GeV...
 
  double r   = 0;
  bool   cut = false;
  if ( pos.z() < 0 )
  {
    r        = sqrt( ( pow( pos.x(), 2 ) + pow( pos.z() + 28000, 2 ) ) );
    double e = 0.45238 * r + 5000;
    cut      = p.e() * m_GeV > e;
  }
  else
  {
    r = sqrt( ( pow( pos.x(), 2 ) + pow( pos.z() - 20000, 2 ) ) );
    if ( r < 15000 ) { cut = true; }
    else
    {
      double e = 0.461538 * ( r - 15000 ) + 10000;
      cut      = p.e() * m_GeV > e;
      //      std::cout<<"z>0 r , e, p.e = "<<r <<" " <<e <<" " <<p.e()*m_GeV<<std::endl;
    }
  }
 
  cut = cut && r < m_rmax;
 
  return cut;
}