Reference Manual

Go to the documentation of this file.
00001 #include "FWCore/Framework/interface/Run.h"
00002 #include "FWCore/ServiceRegistry/interface/Service.h"
00003 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00004 
00005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00006 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
00007 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
00008 
00009 #include "GeneratorInterface/CosmicMuonGenerator/interface/CosMuoGenSource.h"
00010 
00011 
00012 edm::CosMuoGenSource::CosMuoGenSource( const ParameterSet & pset, InputSourceDescription const& desc ) :
00013   GeneratedInputSource(pset, desc ) ,  
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 
00053     //if not specified (i.e. negative) then use MinP also for MinP_CMS
00054     if(MinP_CMS < 0) MinP_CMS = MinP;
00055 
00056     //get seed now from Framework
00057     edm::Service<edm::RandomNumberGenerator> rng;
00058     RanS = rng->mySeed();
00059     // set up the generator
00060     CosMuoGen = new CosmicMuonGenerator();
00061     CosMuoGen->setNumberOfEvents(numberEventsInRun());
00062     CosMuoGen->setRanSeed(RanS);
00063     CosMuoGen->setMinP(MinP);
00064     CosMuoGen->setMinP_CMS(MinP_CMS);
00065     CosMuoGen->setMaxP(MaxP);
00066     CosMuoGen->setMinTheta(MinT);
00067     CosMuoGen->setMaxTheta(MaxT);
00068     CosMuoGen->setMinPhi(MinPh);
00069     CosMuoGen->setMaxPhi(MaxPh);
00070     CosMuoGen->setMinT0(MinS);
00071     CosMuoGen->setMaxT0(MaxS);
00072     CosMuoGen->setElossScaleFactor(ELSF);
00073     CosMuoGen->setRadiusOfTarget(RTarget);
00074     CosMuoGen->setZDistOfTarget(ZTarget);
00075     CosMuoGen->setZCentrOfTarget(ZCTarget);
00076     CosMuoGen->setTrackerOnly(TrackerOnly);
00077     CosMuoGen->setMultiMuon(MultiMuon);
00078     CosMuoGen->setMultiMuonFileName(MultiMuonFileName);
00079     CosMuoGen->setMultiMuonFileFirstEvent(MultiMuonFileFirstEvent);
00080     CosMuoGen->setMultiMuonNmin(MultiMuonNmin);
00081     CosMuoGen->setTIFOnly_constant(TIFOnly_constant);
00082     CosMuoGen->setTIFOnly_linear(TIFOnly_linear);
00083     CosMuoGen->setMTCCHalf(MTCCHalf);
00084     CosMuoGen->setPlugVx(PlugVtx);
00085     CosMuoGen->setPlugVz(PlugVtz);    
00086     CosMuoGen->setRhoAir(VarRhoAir);
00087     CosMuoGen->setRhoWall(VarRhoWall);
00088     CosMuoGen->setRhoRock(VarRhoRock);
00089     CosMuoGen->setRhoClay(VarRhoClay);
00090     CosMuoGen->setRhoPlug(VarRhoPlug);
00091     CosMuoGen->setClayWidth(ClayLayerWidth);
00092     CosMuoGen->setMinEnu(MinEn);
00093     CosMuoGen->setMaxEnu(MaxEn);    
00094     CosMuoGen->setNuProdAlt(NuPrdAlt);
00095     CosMuoGen->setAcptAllMu(AllMu);
00096     CosMuoGen->initialize();
00097     produces<HepMCProduct>();
00098     //  fEvt = new HepMC::GenEvent();
00099     produces<GenRunInfoProduct, edm::InRun>();
00100   }
00101 
00102 edm::CosMuoGenSource::~CosMuoGenSource(){
00103   //CosMuoGen->terminate();
00104   delete CosMuoGen;
00105   //  delete fEvt;
00106   clear();
00107 }
00108 
00109 
00110 void edm::CosMuoGenSource::endRun(Run & r) {
00111 
00112   std::auto_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());
00113 
00114   double cs = CosMuoGen->getRate(); // flux in Hz, not s^-1m^-2
00115   genRunInfo->setInternalXSec(cs);
00116   genRunInfo->setExternalXSecLO(extCrossSect);
00117   genRunInfo->setFilterEfficiency(extFilterEff);
00118 
00119   r.put(genRunInfo);
00120   
00121   CosMuoGen->terminate();
00122 }
00123 
00124 
00125 void edm::CosMuoGenSource::clear(){}
00126 
00127 bool edm::CosMuoGenSource::produce(Event &e)
00128 {  
00129   // generate event
00130   if (!MultiMuon) {
00131     CosMuoGen->nextEvent();
00132   }
00133   else {
00134     bool success = CosMuoGen->nextMultiEvent();
00135     if (!success) return false;
00136   }
00137 
00138   if (Debug) {
00139     std::cout << "CosMuoGenSource.cc: CosMuoGen->EventWeight=" << CosMuoGen->EventWeight 
00140               << "  CosMuoGen: Nmuons=" << CosMuoGen->Id_sf.size() << std::endl; 
00141     std::cout << "CosMuoGen->Id_at=" << CosMuoGen->Id_at
00142               << "  CosMuoGen->Vx_at=" << CosMuoGen->Vx_at 
00143               << "  CosMuoGen->Vy_at=" << CosMuoGen->Vy_at
00144               << "  CosMuoGen->Vz_at=" << CosMuoGen->Vz_at 
00145               << "  CosMuoGen->T0_at=" << CosMuoGen->T0_at << std::endl;
00146     std::cout << "  Px=" << CosMuoGen->Px_at
00147               << "  Py=" << CosMuoGen->Py_at
00148               << "  Pz=" << CosMuoGen->Pz_at << std::endl;
00149     for (unsigned int i=0; i<CosMuoGen->Id_sf.size(); ++i) {
00150       std::cout << "Id_sf[" << i << "]=" << CosMuoGen->Id_sf[i]
00151                 << "  Vx_sf[" << i << "]=" << CosMuoGen->Vx_sf[i]
00152                 << "  Vy_sf=" << CosMuoGen->Vy_sf[i]
00153                 << "  Vz_sf=" << CosMuoGen->Vz_sf[i]
00154                 << "  T0_sf=" << CosMuoGen->T0_sf[i]
00155                 << "  Px_sf=" << CosMuoGen->Px_sf[i]
00156                 << "  Py_sf=" << CosMuoGen->Py_sf[i]
00157                 << "  Pz_sf=" << CosMuoGen->Pz_sf[i] << std::endl;
00158       std::cout << "phi_sf=" << atan2(CosMuoGen->Px_sf[i],CosMuoGen->Pz_sf[i]) << std::endl;
00159       std::cout << "Id_ug[" << i << "]=" << CosMuoGen->Id_ug[i] 
00160                 << "  Vx_ug[" << i << "]=" << CosMuoGen->Vx_ug[i] 
00161                 << "  Vy_ug=" << CosMuoGen->Vy_ug[i]
00162                 << "  Vz_ug=" << CosMuoGen->Vz_ug[i]
00163                 << "  T0_ug=" << CosMuoGen->T0_ug[i]
00164                 << "  Px_ug=" << CosMuoGen->Px_ug[i]
00165                 << "  Py_ug=" << CosMuoGen->Py_ug[i]
00166                 << "  Pz_ug=" << CosMuoGen->Pz_ug[i] << std::endl;
00167       std::cout << "phi_ug=" << atan2(CosMuoGen->Px_ug[i],CosMuoGen->Pz_ug[i]) << std::endl;;
00168     }
00169   }
00170 
00171 
00172   fEvt = new HepMC::GenEvent();
00173   
00174   HepMC::GenVertex* Vtx_at = new  HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_at, //[mm]
00175                                                              CosMuoGen->Vy_at, //[mm]
00176                                                              CosMuoGen->Vz_at, //[mm]
00177                                                              CosMuoGen->T0_at)); //[mm]
00178   //cout << "CosMuoGenSource.cc: Vy_at=" << CosMuoGen->Vy_at << endl;
00179   HepMC::FourVector p_at(CosMuoGen->Px_at,CosMuoGen->Py_at,CosMuoGen->Pz_at,CosMuoGen->E_at);
00180   HepMC::GenParticle* Part_at =
00181     new HepMC::GenParticle(p_at,CosMuoGen->Id_at, 3);//Comment mother particle in
00182   Vtx_at->add_particle_in(Part_at);
00183 
00184 
00185   //loop here in case of multi muon events (else just one iteration)
00186   for (unsigned int i=0; i<CosMuoGen->Id_sf.size(); ++i) {
00187 
00188     HepMC::FourVector p_sf(CosMuoGen->Px_sf[i],CosMuoGen->Py_sf[i],CosMuoGen->Pz_sf[i],CosMuoGen->E_sf[i]);
00189     HepMC::GenParticle* Part_sf_in =
00190       new HepMC::GenParticle(p_sf,CosMuoGen->Id_sf[i], 3); //Comment daughter particle
00191     Vtx_at->add_particle_out(Part_sf_in);
00192     
00193     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]
00194     HepMC::GenParticle* Part_sf_out =
00195       new HepMC::GenParticle(p_sf,CosMuoGen->Id_sf[i], 3); //Comment daughter particle
00196     
00197     Vtx_sf->add_particle_in(Part_sf_in);
00198     Vtx_sf->add_particle_out(Part_sf_out);
00199     
00200     fEvt->add_vertex(Vtx_sf); //one per muon
00201 
00202     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]
00203     
00204     HepMC::FourVector p_ug(CosMuoGen->Px_ug[i],CosMuoGen->Py_ug[i],CosMuoGen->Pz_ug[i],CosMuoGen->E_ug[i]);
00205     HepMC::GenParticle* Part_ug =
00206       new HepMC::GenParticle(p_ug,CosMuoGen->Id_ug[i], 1);//Final state daughter particle
00207 
00208     Vtx_ug->add_particle_in(Part_sf_out);
00209     Vtx_ug->add_particle_out(Part_ug);
00210 
00211     fEvt->add_vertex(Vtx_ug); //one per muon
00212 
00213   }
00214 
00215   fEvt->add_vertex(Vtx_at);
00216   fEvt->set_signal_process_vertex(Vtx_at);
00217 
00218   fEvt->set_event_number(event());
00219   fEvt->set_signal_process_id(13);
00220 
00221   fEvt->weights().push_back( CosMuoGen->EventWeight ); // just one event weight 
00222   fEvt->weights().push_back( CosMuoGen->Trials ); // int Trials number (unweighted) 
00223 
00224 
00225   if (cmVerbosity_) fEvt->print();
00226   
00227   std::auto_ptr<HepMCProduct> CMProduct(new HepMCProduct());
00228   CMProduct->addHepMCData( fEvt );
00229   e.put(CMProduct);
00230 
00231   return true;
00232 }
00233
CMSSW_4_4_3_patch1/src/GeneratorInterface/CosmicMuonGenerator/src/CosMuoGenSource.cc
