CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/IOMC/ParticleGuns/src/FlatEGunASCIIWriter.cc

Go to the documentation of this file.
00001 // ----------------------------------------------------------------------
00002 // FlatEGunASCIIWriter.cc
00003 // Author: Julia Yarba
00004 //
00005 // This code has been molded after examples in HepMC and HepPDT, and
00006 // after single particle gun example (private contacts with Lynn Garren)
00007 //
00008 // Plus, it uses the ParameterSet funtionalities for "user interface"  
00009 //
00010 // It creates a ParticleDataTable from PDG-2004 data file.
00011 // It then creates one or several HepMC particle(s) and adds it(them) 
00012 // to the HepMC::GenEvent.
00013 // For particle definition it uses flat random gun in a given eta, phi 
00014 // and energy ranges to calculate particle kinematics.
00015 // Vertex smearing is currently off, can be easily added later.
00016 // After all, it writes out the event into ASCII output file.
00017 //
00018 // ----------------------------------------------------------------------
00019 
00020 #include <iostream>
00021 #include <fstream>
00022 // header
00023 #include "IOMC/ParticleGuns/interface/FlatEGunASCIIWriter.h"
00024 
00025 // essentials !!!
00026 #include "FWCore/Framework/interface/Event.h"
00027 #include "FWCore/Framework/interface/MakerMacros.h"
00028 #include "FWCore/Utilities/interface/Exception.h"
00029 
00030 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00031 
00032 #include "CLHEP/Random/RandFlat.h"
00033 
00034 using namespace edm;
00035 using namespace std;
00036       
00037 FlatEGunASCIIWriter::FlatEGunASCIIWriter( const ParameterSet& pset )
00038    : fEvt(0), 
00039      fOutFileName( pset.getUntrackedParameter<string>("OutFileName","FlatEGunHepMC.dat") ),
00040      fCurrentEvent(0)
00041 {
00042       
00043   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters") ;
00044   fPartIDs    = pgun_params.getParameter< vector<int> >("PartID");
00045   fMinEta     = pgun_params.getParameter<double>("MinEta");
00046   fMaxEta     = pgun_params.getParameter<double>("MaxEta");
00047   fMinPhi     = pgun_params.getParameter<double>("MinPhi");
00048   fMaxPhi     = pgun_params.getParameter<double>("MaxPhi");
00049   fMinE       = pgun_params.getParameter<double>("MinE");
00050   fMaxE       = pgun_params.getParameter<double>("MaxE");
00051             
00052 }
00053    
00054 FlatEGunASCIIWriter::~FlatEGunASCIIWriter()
00055 {
00056    
00057   if ( fEvt != NULL ) delete fEvt ;
00058   if( fOutStream) delete fOutStream;
00059 
00060 }
00061 
00062 void FlatEGunASCIIWriter::beginRun( edm::Run &r, const EventSetup& es)
00063 {
00064 
00065   es.getData( fPDGTable );
00066 
00067 }
00068 
00069 void FlatEGunASCIIWriter::beginJob()
00070 { 
00071   fOutStream = new HepMC::IO_GenEvent( fOutFileName.c_str() ); 
00072   if ( fOutStream->rdstate() == std::ios::failbit ) {
00073     throw cms::Exception("FileNotOpen", "FlatEGunASCIIWriter::beginJob()")
00074       << "File " << fOutFileName << " was not open.\n";
00075   }
00076   
00077   return ;
00078 
00079 }
00080 
00081 void FlatEGunASCIIWriter::analyze( const Event& , 
00082                                    const EventSetup& /* not used so far, thus is "empty" */ ) 
00083 {
00084          
00085    // clean up GenEvent memory : also deletes all vtx/part in it
00086    // 
00087    if ( fEvt != NULL ) delete fEvt ;
00088    
00089    // here re-create fEvt (memory)
00090    fEvt = new HepMC::GenEvent() ;
00091 
00092    // now actualy, cook up the event from PDGTable and gun parameters
00093    //
00094    //HepMC::GenVertex* Vtx = new HepMC::GenVertex( CLHEP::HepLorentzVector(0.,0.,0.) );
00095    HepMC::GenVertex* Vtx = new HepMC::GenVertex( HepMC::FourVector(0.,0.,0.) );
00096 
00097    // loop over particles
00098    //
00099    for ( unsigned int ip=0; ip<fPartIDs.size(); ip++ )
00100    {
00101        double energy = CLHEP::RandFlat::shoot( fMinE, fMaxE ) ;
00102        double eta    = CLHEP::RandFlat::shoot( fMinEta, fMaxEta ) ;
00103        double phi    = CLHEP::RandFlat::shoot( fMinPhi, fMaxPhi ) ;
00104        const HepPDT::ParticleData* 
00105           PData = fPDGTable->particle(HepPDT::ParticleID(abs(fPartIDs[ip]))) ;
00106        double mass   = PData->mass().value() ;
00107        double mom2   = energy*energy - mass*mass ;
00108        double mom    = 0. ;
00109        if ( mom2 > 0. ) mom = sqrt(mom2) ;
00110        double theta  = 2.*atan(exp(-eta)) ;
00111        double px     = mom*sin(theta)*cos(phi) ;
00112        double py     = mom*sin(theta)*sin(phi) ;
00113        double pz     = mom*cos(theta) ;
00114        //CLHEP::Hep3Vector p(px,py,pz) ;
00115        //HepMC::GenParticle* Part = 
00116        //    new HepMC::GenParticle(CLHEP::HepLorentzVector(p,energy),fPartIDs[ip],1);
00117        HepMC::FourVector p(px,py,pz,energy);
00118        HepMC::GenParticle* Part = 
00119            new HepMC::GenParticle(p,fPartIDs[ip],1);
00120        Vtx->add_particle_out(Part);
00121    }
00122    fEvt->add_vertex( Vtx ) ;
00123    fEvt->set_event_number( fCurrentEvent+1 ) ;
00124    fEvt->set_signal_process_id(20) ; 
00125    
00126    // for testing purpose only
00127    // fEvt->print() ;     
00128 
00129    fOutStream->write_event( fEvt );
00130 
00131    fCurrentEvent++ ;
00132 
00133    return ;
00134    
00135 }
00136