Reference Manual

Go to the documentation of this file.
00001 #include <CLHEP/Random/RandomEngine.h>
00002 
00003 #include "FWCore/Framework/interface/Run.h"
00004 #include "FWCore/ServiceRegistry/interface/Service.h"
00005 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00006 
00007 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00008 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
00009 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
00010 
00011 #include "GeneratorInterface/CosmicMuonGenerator/interface/CosMuoGenProducer.h"
00012 
00013 edm::CosMuoGenProducer::CosMuoGenProducer( const ParameterSet & pset ) :
00014   //RanS(pset.getParameter<int>("RanSeed", 123456)), //get seed now from Framework
00015   MinP(pset.getParameter<double>("MinP")),
00016   MinP_CMS(pset.getParameter<double>("MinP_CMS")),
00017   MaxP(pset.getParameter<double>("MaxP")),
00018   MinT(pset.getParameter<double>("MinTheta")),
00019   MaxT(pset.getParameter<double>("MaxTheta")),
00020   MinPh(pset.getParameter<double>("MinPhi")),
00021   MaxPh(pset.getParameter<double>("MaxPhi")),
00022   MinS(pset.getParameter<double>("MinT0")),
00023   MaxS(pset.getParameter<double>("MaxT0")),
00024   ELSF(pset.getParameter<double>("ElossScaleFactor")),
00025   RTarget(pset.getParameter<double>("RadiusOfTarget")),
00026   ZTarget(pset.getParameter<double>("ZDistOfTarget")),
00027   ZCTarget(pset.getParameter<double>("ZCentrOfTarget")),
00028   TrackerOnly(pset.getParameter<bool>("TrackerOnly")),
00029   MultiMuon(pset.getParameter<bool>("MultiMuon")),
00030   MultiMuonFileName(pset.getParameter<std::string>("MultiMuonFileName")),
00031   MultiMuonFileFirstEvent(pset.getParameter<int>("MultiMuonFileFirstEvent")),
00032   MultiMuonNmin(pset.getParameter<int>("MultiMuonNmin")),
00033   TIFOnly_constant(pset.getParameter<bool>("TIFOnly_constant")),
00034   TIFOnly_linear(pset.getParameter<bool>("TIFOnly_linear")),
00035   MTCCHalf(pset.getParameter<bool>("MTCCHalf")),
00036   PlugVtx(pset.getParameter<double>("PlugVx")),
00037   PlugVtz(pset.getParameter<double>("PlugVz")),
00038   VarRhoAir(pset.getParameter<double>("RhoAir")),
00039   VarRhoWall(pset.getParameter<double>("RhoWall")),
00040   VarRhoRock(pset.getParameter<double>("RhoRock")),
00041   VarRhoClay(pset.getParameter<double>("RhoClay")),
00042   VarRhoPlug(pset.getParameter<double>("RhoPlug")),
00043   ClayLayerWidth(pset.getParameter<double>("ClayWidth")),
00044   MinEn(pset.getParameter<double>("MinEnu")),
00045   MaxEn(pset.getParameter<double>("MaxEnu")),
00046   NuPrdAlt(pset.getParameter<double>("NuProdAlt")),
00047   AllMu(pset.getParameter<bool>("AcptAllMu")),
00048   extCrossSect(pset.getUntrackedParameter<double>("crossSection", -1.)),
00049   extFilterEff(pset.getUntrackedParameter<double>("filterEfficiency", -1.)),
00050   cmVerbosity_(pset.getParameter<bool>("Verbosity"))
00051   {
00052     //if not specified (i.e. negative) then use MinP also for MinP_CMS
00053     if(MinP_CMS < 0) MinP_CMS = MinP;
00054 
00055     edm::Service<RandomNumberGenerator> rng;
00056     if (!rng.isAvailable())
00057       throw cms::Exception("Configuration")
00058          << "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
00059             "which appears to be absent.  Please add that service to your configuration\n"
00060             "or remove the modules that require it." << std::endl;
00061 
00062     // set up the generator
00063     CosMuoGen = new CosmicMuonGenerator();
00064 // Begin JMM change
00065 //  CosMuoGen->setNumberOfEvents(numberEventsInRun());
00066     CosMuoGen->setNumberOfEvents(999999999);
00067 // End of JMM change
00068     CosMuoGen->setRanSeed(RanS);
00069     CosMuoGen->setMinP(MinP);
00070     CosMuoGen->setMinP_CMS(MinP_CMS);
00071     CosMuoGen->setMaxP(MaxP);
00072     CosMuoGen->setMinTheta(MinT);
00073     CosMuoGen->setMaxTheta(MaxT);
00074     CosMuoGen->setMinPhi(MinPh);
00075     CosMuoGen->setMaxPhi(MaxPh);
00076     CosMuoGen->setMinT0(MinS);
00077     CosMuoGen->setMaxT0(MaxS);
00078     CosMuoGen->setElossScaleFactor(ELSF);
00079     CosMuoGen->setRadiusOfTarget(RTarget);
00080     CosMuoGen->setZDistOfTarget(ZTarget);
00081     CosMuoGen->setZCentrOfTarget(ZCTarget);
00082     CosMuoGen->setTrackerOnly(TrackerOnly);
00083     CosMuoGen->setMultiMuon(MultiMuon);
00084     CosMuoGen->setMultiMuonFileName(MultiMuonFileName);
00085     CosMuoGen->setMultiMuonFileFirstEvent(MultiMuonFileFirstEvent);
00086     CosMuoGen->setMultiMuonNmin(MultiMuonNmin);
00087     CosMuoGen->setTIFOnly_constant(TIFOnly_constant);
00088     CosMuoGen->setTIFOnly_linear(TIFOnly_linear);
00089     CosMuoGen->setMTCCHalf(MTCCHalf);
00090     CosMuoGen->setPlugVx(PlugVtx);
00091     CosMuoGen->setPlugVz(PlugVtz);    
00092     CosMuoGen->setRhoAir(VarRhoAir);
00093     CosMuoGen->setRhoWall(VarRhoWall);
00094     CosMuoGen->setRhoRock(VarRhoRock);
00095     CosMuoGen->setRhoClay(VarRhoClay);
00096     CosMuoGen->setRhoPlug(VarRhoPlug);
00097     CosMuoGen->setClayWidth(ClayLayerWidth);
00098     CosMuoGen->setMinEnu(MinEn);
00099     CosMuoGen->setMaxEnu(MaxEn);    
00100     CosMuoGen->setNuProdAlt(NuPrdAlt);
00101     CosMuoGen->setAcptAllMu(AllMu);
00102     CosMuoGen->initialize(&rng->getEngine());
00103     produces<HepMCProduct>();
00104     produces<GenEventInfoProduct>();
00105     produces<GenRunInfoProduct, edm::InRun>();
00106   }
00107 
00108 edm::CosMuoGenProducer::~CosMuoGenProducer(){
00109   //CosMuoGen->terminate();
00110   delete CosMuoGen;
00111   //  delete fEvt;
00112   clear();
00113 }
00114 
00115 void edm::CosMuoGenProducer::endRun( Run &run, const EventSetup& es )
00116 {
00117   std::auto_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());
00118 
00119   double cs = CosMuoGen->getRate(); // flux in Hz, not s^-1m^-2
00120   if (MultiMuon) genRunInfo->setInternalXSec(0.);
00121   else genRunInfo->setInternalXSec(cs);
00122   genRunInfo->setExternalXSecLO(extCrossSect);
00123   genRunInfo->setFilterEfficiency(extFilterEff);
00124 
00125   run.put(genRunInfo);
00126 
00127   CosMuoGen->terminate();
00128 }
00129 
00130 void edm::CosMuoGenProducer::clear(){}
00131 
00132 void edm::CosMuoGenProducer::produce(Event &e, const edm::EventSetup &es)
00133 {  
00134   // generate event
00135   if (!MultiMuon) {
00136     CosMuoGen->nextEvent();
00137   }
00138   else {
00139     bool success = CosMuoGen->nextMultiEvent();
00140     if (!success) std::cout << "CosMuoGenProducer.cc: CosMuoGen->nextMultiEvent() failed!" 
00141                             << std::endl;
00142   }
00143 
00144   if (Debug) {
00145     std::cout << "CosMuoGenSource.cc: CosMuoGen->EventWeight=" << CosMuoGen->EventWeight 
00146               << "  CosMuoGen: Nmuons=" << CosMuoGen->Id_sf.size() << std::endl; 
00147     std::cout << "CosMuoGen->Id_at=" << CosMuoGen->Id_at
00148               << "  CosMuoGen->Vx_at=" << CosMuoGen->Vx_at 
00149               << "  CosMuoGen->Vy_at=" << CosMuoGen->Vy_at
00150               << "  CosMuoGen->Vz_at=" << CosMuoGen->Vz_at 
00151               << "  CosMuoGen->T0_at=" << CosMuoGen->T0_at << std::endl;
00152     std::cout << "  Px=" << CosMuoGen->Px_at
00153               << "  Py=" << CosMuoGen->Py_at
00154               << "  Pz=" << CosMuoGen->Pz_at << std::endl;
00155     for (unsigned int i=0; i<CosMuoGen->Id_sf.size(); ++i) {
00156       std::cout << "Id_sf[" << i << "]=" << CosMuoGen->Id_sf[i]
00157                 << "  Vx_sf[" << i << "]=" << CosMuoGen->Vx_sf[i]
00158                 << "  Vy_sf=" << CosMuoGen->Vy_sf[i]
00159                 << "  Vz_sf=" << CosMuoGen->Vz_sf[i]
00160                 << "  T0_sf=" << CosMuoGen->T0_sf[i]
00161                 << "  Px_sf=" << CosMuoGen->Px_sf[i]
00162                 << "  Py_sf=" << CosMuoGen->Py_sf[i]
00163                 << "  Pz_sf=" << CosMuoGen->Pz_sf[i] << std::endl;
00164       std::cout << "phi_sf=" << atan2(CosMuoGen->Px_sf[i],CosMuoGen->Pz_sf[i]) << std::endl;
00165       std::cout << "Id_ug[" << i << "]=" << CosMuoGen->Id_ug[i] 
00166                 << "  Vx_ug[" << i << "]=" << CosMuoGen->Vx_ug[i] 
00167                 << "  Vy_ug=" << CosMuoGen->Vy_ug[i]
00168                 << "  Vz_ug=" << CosMuoGen->Vz_ug[i]
00169                 << "  T0_ug=" << CosMuoGen->T0_ug[i]
00170                 << "  Px_ug=" << CosMuoGen->Px_ug[i]
00171                 << "  Py_ug=" << CosMuoGen->Py_ug[i]
00172                 << "  Pz_ug=" << CosMuoGen->Pz_ug[i] << std::endl;
00173       std::cout << "phi_ug=" << atan2(CosMuoGen->Px_ug[i],CosMuoGen->Pz_ug[i]) << std::endl;;
00174     }
00175   }
00176 
00177 
00178   fEvt = new HepMC::GenEvent();
00179   
00180   HepMC::GenVertex* Vtx_at = new  HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_at, //[mm]
00181                                                              CosMuoGen->Vy_at, //[mm]
00182                                                              CosMuoGen->Vz_at, //[mm]
00183                                                              CosMuoGen->T0_at)); //[mm]
00184   //cout << "CosMuoGenSource.cc: Vy_at=" << CosMuoGen->Vy_at << endl;
00185   HepMC::FourVector p_at(CosMuoGen->Px_at,CosMuoGen->Py_at,CosMuoGen->Pz_at,CosMuoGen->E_at);
00186   HepMC::GenParticle* Part_at =
00187     new HepMC::GenParticle(p_at,CosMuoGen->Id_at, 3);//Comment mother particle in
00188   Vtx_at->add_particle_in(Part_at);
00189 
00190 
00191   //loop here in case of multi muon events (else just one iteration)
00192   for (unsigned int i=0; i<CosMuoGen->Id_sf.size(); ++i) {
00193 
00194     HepMC::FourVector p_sf(CosMuoGen->Px_sf[i],CosMuoGen->Py_sf[i],CosMuoGen->Pz_sf[i],CosMuoGen->E_sf[i]);
00195     HepMC::GenParticle* Part_sf_in =
00196       new HepMC::GenParticle(p_sf,CosMuoGen->Id_sf[i], 3); //Comment daughter particle
00197     Vtx_at->add_particle_out(Part_sf_in);
00198     
00199     HepMC::GenVertex* Vtx_sf = new HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_sf[i],                             CosMuoGen->Vy_sf[i], CosMuoGen->Vz_sf[i], CosMuoGen->T0_sf[i])); //[mm]
00200     HepMC::GenParticle* Part_sf_out =
00201       new HepMC::GenParticle(p_sf,CosMuoGen->Id_sf[i], 3); //Comment daughter particle
00202     
00203     Vtx_sf->add_particle_in(Part_sf_in);
00204     Vtx_sf->add_particle_out(Part_sf_out);
00205     
00206     fEvt->add_vertex(Vtx_sf); //one per muon
00207 
00208     HepMC::GenVertex* Vtx_ug = new HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_ug[i],                             CosMuoGen->Vy_ug[i], CosMuoGen->Vz_ug[i], CosMuoGen->T0_ug[i])); //[mm]
00209     
00210     HepMC::FourVector p_ug(CosMuoGen->Px_ug[i],CosMuoGen->Py_ug[i],CosMuoGen->Pz_ug[i],CosMuoGen->E_ug[i]);
00211     HepMC::GenParticle* Part_ug =
00212       new HepMC::GenParticle(p_ug,CosMuoGen->Id_ug[i], 1);//Final state daughter particle
00213 
00214     Vtx_ug->add_particle_in(Part_sf_out);
00215     Vtx_ug->add_particle_out(Part_ug);
00216 
00217     fEvt->add_vertex(Vtx_ug); //one per muon
00218 
00219   }
00220 
00221   fEvt->add_vertex(Vtx_at);
00222   fEvt->set_signal_process_vertex(Vtx_at);
00223 
00224   fEvt->set_event_number(e.id().event());
00225   fEvt->set_signal_process_id(13);
00226 
00227   fEvt->weights().push_back( CosMuoGen->EventWeight ); // just one event weight 
00228   fEvt->weights().push_back( CosMuoGen->Trials ); // int Trials number (unweighted) 
00229 
00230 
00231   if (cmVerbosity_) fEvt->print();
00232 
00233   std::auto_ptr<HepMCProduct> CMProduct(new HepMCProduct());
00234   CMProduct->addHepMCData( fEvt );
00235   e.put(CMProduct);
00236 
00237   std::auto_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct( fEvt ));
00238   e.put(genEventInfo);
00239 
00240 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/GeneratorInterface/CosmicMuonGenerator/src/CosMuoGenProducer.cc
