CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/IOMC/ParticleGuns/src/MultiParticleInConeGunSource.cc

Go to the documentation of this file.
00001 /*
00002  *  $Date: 2010/10/03 10:25:53 $
00003  *  $Revision: 1.5 $
00004  *  \author Jean-Roch Vlimant
00005  */
00006 
00007 #include <ostream>
00008 
00009 #include "IOMC/ParticleGuns/interface/MultiParticleInConeGunSource.h"
00010 
00011 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00012 
00013 #include "FWCore/Framework/interface/Event.h"
00014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00015 
00016 // #include "FWCore/Utilities/interface/Exception.h"
00017 
00018 // #include "CLHEP/Random/RandFlat.h"
00019 
00020 using namespace edm;
00021 using namespace std;
00022 
00023 MultiParticleInConeGunSource::MultiParticleInConeGunSource(const ParameterSet& pset,
00024                                              const InputSourceDescription& desc) : 
00025    BaseFlatGunSource(pset, desc)
00026 {
00027 
00028 
00029    ParameterSet defpset ;
00030    ParameterSet pgun_params = 
00031       pset.getUntrackedParameter<ParameterSet>("PGunParameters",defpset) ;
00032   
00033    fMinPt = pgun_params.getUntrackedParameter<double>("MinPt",0.99);
00034    fMaxPt = pgun_params.getUntrackedParameter<double>("MaxPt",1.01);
00035 
00036    fInConeIds = pgun_params.getUntrackedParameter< vector<int> >("InConeID",vector<int>());
00037    fMinDeltaR = pgun_params.getUntrackedParameter<double>("MinDeltaR",0.0);
00038    fMaxDeltaR = pgun_params.getUntrackedParameter<double>("MaxDeltaR",0.1);
00039    fMinMomRatio = pgun_params.getUntrackedParameter<double>("MinMomRatio",0.99);
00040    fMaxMomRatio = pgun_params.getUntrackedParameter<double>("MaxMomRatio",1.01);
00041 
00042    fInConeMinEta = pgun_params.getUntrackedParameter<double>("InConeMinEta",-5.5);
00043    fInConeMaxEta = pgun_params.getUntrackedParameter<double>("InConeMaxEta",5.5);
00044    fInConeMinPhi = pgun_params.getUntrackedParameter<double>("InConeMinPhi",-3.14159265358979323846);
00045    fInConeMaxPhi = pgun_params.getUntrackedParameter<double>("InConeMaxPhi",3.14159265358979323846);
00046    fInConeMaxTry = pgun_params.getUntrackedParameter<unsigned int>("InConeMaxTry",100);
00047    
00048    produces<HepMCProduct>();
00049 }
00050 
00051 MultiParticleInConeGunSource::~MultiParticleInConeGunSource()
00052 {
00053    // no need to cleanup GenEvent memory - done in HepMCProduct
00054 }
00055 
00056 bool MultiParticleInConeGunSource::produce(Event &e) 
00057 {
00058 
00059    if ( fVerbosity > 0 )
00060    {
00061       cout << " MultiParticleInConeGunSource : Begin New Event Generation" << endl ; 
00062    }
00063    // event loop (well, another step in it...)
00064           
00065    // no need to clean up GenEvent memory - done in HepMCProduct
00066    // 
00067    
00068    // here re-create fEvt (memory)
00069    //
00070    fEvt = new HepMC::GenEvent() ;
00071    
00072    // now actualy, cook up the event from PDGTable and gun parameters
00073    //
00074    // 1st, primary vertex
00075    //
00076    //HepMC::GenVertex* Vtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.,0.,0.));
00077    HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0.,0.,0.));
00078 
00079    // loop over particles
00080    //
00081    int barcode = 1 ;
00082    for (unsigned int ip=0; ip<fPartIDs.size(); ++ip)
00083    {
00084 
00085        double pt     = fRandomGenerator->fire(fMinPt, fMaxPt) ;
00086        double eta    = fRandomGenerator->fire(fMinEta, fMaxEta) ;
00087        double phi    = fRandomGenerator->fire(fMinPhi, fMaxPhi) ;
00088        int PartID = fPartIDs[ip] ;
00089        const HepPDT::ParticleData* 
00090           PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
00091        double mass   = PData->mass().value() ;
00092        double theta  = 2.*atan(exp(-eta)) ;
00093        double mom    = pt/sin(theta) ;
00094        double px     = pt*cos(phi) ;
00095        double py     = pt*sin(phi) ;
00096        double pz     = mom*cos(theta) ;
00097        double energy2= mom*mom + mass*mass ;
00098        double energy = sqrt(energy2) ; 
00099 
00100        HepMC::FourVector p(px,py,pz,energy) ;
00101        HepMC::GenParticle* Part = 
00102            new HepMC::GenParticle(p,PartID,1);
00103        Part->suggest_barcode( barcode ) ;
00104        barcode++ ;
00105        Vtx->add_particle_out(Part);
00106 
00107        if ( fAddAntiParticle ){}
00108 
00109        // now add the particles in the cone
00110        for (unsigned iPic=0; iPic!=fInConeIds.size();iPic++){
00111          unsigned int nTry=0;
00112          while(true){
00113            //shoot flat Deltar
00114            double dR = fRandomGenerator->fire(fMinDeltaR, fMaxDeltaR);
00115            //shoot flat eta/phi mixing
00116            double alpha = fRandomGenerator->fire(-3.14159265358979323846, 3.14159265358979323846);
00117            double dEta = dR*cos(alpha);
00118            double dPhi = dR*sin(alpha);
00119            
00120            /*
00121            //shoot Energy of associated particle         
00122            double energyIc = fRandomGenerator->fire(fMinEInCone, fMaxEInCone);
00123            if (mom2Ic>0){ momIC = sqrt(mom2Ic);}
00124            */
00125            //    get kinematics
00126            double etaIc = eta+dEta;
00127            double phiIc = phi+dPhi;
00128            //put it back in -Pi:Pi if necessary. multiple time might be necessary if dR > 3
00129            const unsigned int maxL=100;
00130            unsigned int iL=0;
00131            while(iL++<maxL){
00132              if (phiIc > 3.14159265358979323846) phiIc-=2*3.14159265358979323846;
00133              else if(phiIc <-3.14159265358979323846) phiIc+=2*3.14159265358979323846;
00134              
00135              if (abs(phiIc)<3.14159265358979323846) break;
00136            }
00137              
00138 
00139            //allow to skip it if you did not run out of possible drawing
00140            if (nTry++<=fInConeMaxTry){
00141              //draw another one if this one is not in acceptance
00142              if (etaIc<fInConeMinEta || etaIc > fInConeMaxEta) continue;
00143              if (phiIc<fInConeMinPhi || phiIc > fInConeMaxPhi) continue;
00144            }
00145            else{
00146              if ( fVerbosity > 0 )
00147                {
00148                  cout << " MultiParticleInConeGunSource : could not produce a particle "
00149                       <<  fInConeIds[iPic]<<" in cone "
00150                       <<  fMinDeltaR<<" to "<<fMaxDeltaR<<" within the "<<fInConeMaxTry<<" allowed tries."<<endl;
00151                }
00152              /*      edm::LogError("MultiParticleInConeGunSource")<< " MultiParticleInConeGunSource : could not produce a particle "<<
00153                fInConeIds[iPic]<<" in cone "<<
00154                fMinDeltaR<<" to "<<fMaxDeltaR<<" within the "<<fInConeMaxTry<<" allowed tries.";*/
00155            }
00156            int PartIDIc=fInConeIds[iPic];
00157            const HepPDT::ParticleData* 
00158              PDataIc = fPDGTable->particle(HepPDT::ParticleID(abs(PartIDIc)));
00159            
00160            //shoot momentum ratio
00161            double momR = fRandomGenerator->fire(fMinMomRatio, fMaxMomRatio);
00162            double massIc= PDataIc->mass().value() ;
00163            double momIc = momR * mom;
00164            double energyIc = sqrt(momIc*momIc + massIc*massIc);
00165 
00166            double thetaIc = 2.*atan(exp(-etaIc)) ;
00167            double pxIc = momIc*sin(thetaIc)*cos(phiIc);
00168            double pyIc = momIc*sin(thetaIc)*sin(phiIc);
00169            double pzIc = momIc*cos(thetaIc);
00170 
00171            HepMC::FourVector pIc(pxIc,pyIc,pzIc,energyIc) ;
00172            HepMC::GenParticle* PartIc = new HepMC::GenParticle(pIc, PartIDIc, 1);
00173            PartIc->suggest_barcode( barcode ) ;
00174            barcode++ ;
00175            Vtx->add_particle_out(PartIc);
00176            break;
00177          }//try many times while not in acceptance
00178        }//loop over the particle Ids in the cone
00179    }
00180 
00181    fEvt->add_vertex(Vtx) ;
00182    fEvt->set_event_number(event()) ;
00183    fEvt->set_signal_process_id(20) ; 
00184         
00185    if ( fVerbosity > 0 )
00186    {
00187       fEvt->print() ;  
00188    }
00189 
00190    auto_ptr<HepMCProduct> BProduct(new HepMCProduct()) ;
00191    BProduct->addHepMCData( fEvt );
00192    e.put(BProduct);
00193    
00194    if ( fVerbosity > 0 )
00195      {
00196        cout << " MultiParticleInConeGunSource : Event Generation Done " << endl;
00197      }
00198 
00199    return true;
00200 }
00201