Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:13:28 

0001 
0002 #include "FWCore/Framework/interface/Run.h"
0003 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0004 
0005 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0006 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0007 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0008 
0009 #include "GeneratorInterface/CosmicMuonGenerator/interface/CosMuoGenProducer.h"
0010 
0011 edm::CosMuoGenProducer::CosMuoGenProducer(const ParameterSet& pset)
0012     :  //RanS(pset.getParameter<int>("RanSeed", 123456)), //get seed now from Framework
0013       MinP(pset.getParameter<double>("MinP")),
0014       MinP_CMS(pset.getParameter<double>("MinP_CMS")),
0015       MaxP(pset.getParameter<double>("MaxP")),
0016       MinT(pset.getParameter<double>("MinTheta")),
0017       MaxT(pset.getParameter<double>("MaxTheta")),
0018       MinPh(pset.getParameter<double>("MinPhi")),
0019       MaxPh(pset.getParameter<double>("MaxPhi")),
0020       MinS(pset.getParameter<double>("MinT0")),
0021       MaxS(pset.getParameter<double>("MaxT0")),
0022       ELSF(pset.getParameter<double>("ElossScaleFactor")),
0023       RTarget(pset.getParameter<double>("RadiusOfTarget")),
0024       ZTarget(pset.getParameter<double>("ZDistOfTarget")),
0025       ZCTarget(pset.getParameter<double>("ZCentrOfTarget")),
0026       TrackerOnly(pset.getParameter<bool>("TrackerOnly")),
0027       MultiMuon(pset.getParameter<bool>("MultiMuon")),
0028       MultiMuonFileName(pset.getParameter<std::string>("MultiMuonFileName")),
0029       MultiMuonFileFirstEvent(pset.getParameter<int>("MultiMuonFileFirstEvent")),
0030       MultiMuonNmin(pset.getParameter<int>("MultiMuonNmin")),
0031       TIFOnly_constant(pset.getParameter<bool>("TIFOnly_constant")),
0032       TIFOnly_linear(pset.getParameter<bool>("TIFOnly_linear")),
0033       MTCCHalf(pset.getParameter<bool>("MTCCHalf")),
0034       PlugVtx(pset.getParameter<double>("PlugVx")),
0035       PlugVtz(pset.getParameter<double>("PlugVz")),
0036       VarRhoAir(pset.getParameter<double>("RhoAir")),
0037       VarRhoWall(pset.getParameter<double>("RhoWall")),
0038       VarRhoRock(pset.getParameter<double>("RhoRock")),
0039       VarRhoClay(pset.getParameter<double>("RhoClay")),
0040       VarRhoPlug(pset.getParameter<double>("RhoPlug")),
0041       ClayLayerWidth(pset.getParameter<double>("ClayWidth")),
0042       MinEn(pset.getParameter<double>("MinEnu")),
0043       MaxEn(pset.getParameter<double>("MaxEnu")),
0044       NuPrdAlt(pset.getParameter<double>("NuProdAlt")),
0045       AllMu(pset.getParameter<bool>("AcptAllMu")),
0046       extCrossSect(pset.getUntrackedParameter<double>("crossSection", -1.)),
0047       extFilterEff(pset.getUntrackedParameter<double>("filterEfficiency", -1.)),
0048       cmVerbosity_(pset.getParameter<bool>("Verbosity")),
0049       isInitialized_(false) {
0050   //if not specified (i.e. negative) then use MinP also for MinP_CMS
0051   if (MinP_CMS < 0)
0052     MinP_CMS = MinP;
0053 
0054   // set up the generator
0055   CosMuoGen = std::make_unique<CosmicMuonGenerator>();
0056   // Begin JMM change
0057   //  CosMuoGen->setNumberOfEvents(numberEventsInRun());
0058   CosMuoGen->setNumberOfEvents(999999999);
0059   // End of JMM change
0060   CosMuoGen->setRanSeed(RanS);
0061   CosMuoGen->setMinP(MinP);
0062   CosMuoGen->setMinP_CMS(MinP_CMS);
0063   CosMuoGen->setMaxP(MaxP);
0064   CosMuoGen->setMinTheta(MinT);
0065   CosMuoGen->setMaxTheta(MaxT);
0066   CosMuoGen->setMinPhi(MinPh);
0067   CosMuoGen->setMaxPhi(MaxPh);
0068   CosMuoGen->setMinT0(MinS);
0069   CosMuoGen->setMaxT0(MaxS);
0070   CosMuoGen->setElossScaleFactor(ELSF);
0071   CosMuoGen->setRadiusOfTarget(RTarget);
0072   CosMuoGen->setZDistOfTarget(ZTarget);
0073   CosMuoGen->setZCentrOfTarget(ZCTarget);
0074   CosMuoGen->setTrackerOnly(TrackerOnly);
0075   CosMuoGen->setMultiMuon(MultiMuon);
0076   CosMuoGen->setMultiMuonFileName(MultiMuonFileName);
0077   CosMuoGen->setMultiMuonFileFirstEvent(MultiMuonFileFirstEvent);
0078   CosMuoGen->setMultiMuonNmin(MultiMuonNmin);
0079   CosMuoGen->setTIFOnly_constant(TIFOnly_constant);
0080   CosMuoGen->setTIFOnly_linear(TIFOnly_linear);
0081   CosMuoGen->setMTCCHalf(MTCCHalf);
0082   CosMuoGen->setPlugVx(PlugVtx);
0083   CosMuoGen->setPlugVz(PlugVtz);
0084   CosMuoGen->setRhoAir(VarRhoAir);
0085   CosMuoGen->setRhoWall(VarRhoWall);
0086   CosMuoGen->setRhoRock(VarRhoRock);
0087   CosMuoGen->setRhoClay(VarRhoClay);
0088   CosMuoGen->setRhoPlug(VarRhoPlug);
0089   CosMuoGen->setClayWidth(ClayLayerWidth);
0090   CosMuoGen->setMinEnu(MinEn);
0091   CosMuoGen->setMaxEnu(MaxEn);
0092   CosMuoGen->setNuProdAlt(NuPrdAlt);
0093   CosMuoGen->setAcptAllMu(AllMu);
0094   produces<HepMCProduct>("unsmeared");
0095   produces<GenEventInfoProduct>();
0096   produces<GenRunInfoProduct, edm::Transition::EndRun>();
0097 }
0098 
0099 edm::CosMuoGenProducer::~CosMuoGenProducer() { clear(); }
0100 
0101 void edm::CosMuoGenProducer::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {
0102   if (!isInitialized_) {
0103     isInitialized_ = true;
0104     RandomEngineSentry<CosmicMuonGenerator> randomEngineSentry(CosMuoGen.get(), lumi.index());
0105     CosMuoGen->initialize(randomEngineSentry.randomEngine());
0106   }
0107 }
0108 
0109 void edm::CosMuoGenProducer::endRunProduce(Run& run, const EventSetup& es) {
0110   std::unique_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());
0111 
0112   double cs = CosMuoGen->getRate();  // flux in Hz, not s^-1m^-2
0113   if (MultiMuon)
0114     genRunInfo->setInternalXSec(0.);
0115   else
0116     genRunInfo->setInternalXSec(cs);
0117   genRunInfo->setExternalXSecLO(extCrossSect);
0118   genRunInfo->setFilterEfficiency(extFilterEff);
0119 
0120   run.put(std::move(genRunInfo));
0121 
0122   CosMuoGen->terminate();
0123 }
0124 
0125 void edm::CosMuoGenProducer::clear() {}
0126 
0127 void edm::CosMuoGenProducer::produce(Event& e, const edm::EventSetup& es) {
0128   RandomEngineSentry<CosmicMuonGenerator> randomEngineSentry(CosMuoGen.get(), e.streamID());
0129 
0130   // generate event
0131   if (!MultiMuon) {
0132     CosMuoGen->nextEvent();
0133   } else {
0134     bool success = CosMuoGen->nextMultiEvent();
0135     if (!success)
0136       std::cout << "CosMuoGenProducer.cc: CosMuoGen->nextMultiEvent() failed!" << std::endl;
0137   }
0138 
0139   if (Debug) {
0140     std::cout << "CosMuoGenProducer.cc: CosMuoGen->EventWeight=" << CosMuoGen->EventWeight
0141               << "  CosMuoGen: Nmuons=" << CosMuoGen->Id_sf.size() << std::endl;
0142     std::cout << "CosMuoGen->Id_at=" << CosMuoGen->Id_at << "  CosMuoGen->Vx_at=" << CosMuoGen->Vx_at
0143               << "  CosMuoGen->Vy_at=" << CosMuoGen->Vy_at << "  CosMuoGen->Vz_at=" << CosMuoGen->Vz_at
0144               << "  CosMuoGen->T0_at=" << CosMuoGen->T0_at << std::endl;
0145     std::cout << "  Px=" << CosMuoGen->Px_at << "  Py=" << CosMuoGen->Py_at << "  Pz=" << CosMuoGen->Pz_at << std::endl;
0146     for (unsigned int i = 0; i < CosMuoGen->Id_sf.size(); ++i) {
0147       std::cout << "Id_sf[" << i << "]=" << CosMuoGen->Id_sf[i] << "  Vx_sf[" << i << "]=" << CosMuoGen->Vx_sf[i]
0148                 << "  Vy_sf=" << CosMuoGen->Vy_sf[i] << "  Vz_sf=" << CosMuoGen->Vz_sf[i]
0149                 << "  T0_sf=" << CosMuoGen->T0_sf[i] << "  Px_sf=" << CosMuoGen->Px_sf[i]
0150                 << "  Py_sf=" << CosMuoGen->Py_sf[i] << "  Pz_sf=" << CosMuoGen->Pz_sf[i] << std::endl;
0151       std::cout << "phi_sf=" << atan2(CosMuoGen->Px_sf[i], CosMuoGen->Pz_sf[i]) << std::endl;
0152       std::cout << "Id_ug[" << i << "]=" << CosMuoGen->Id_ug[i] << "  Vx_ug[" << i << "]=" << CosMuoGen->Vx_ug[i]
0153                 << "  Vy_ug=" << CosMuoGen->Vy_ug[i] << "  Vz_ug=" << CosMuoGen->Vz_ug[i]
0154                 << "  T0_ug=" << CosMuoGen->T0_ug[i] << "  Px_ug=" << CosMuoGen->Px_ug[i]
0155                 << "  Py_ug=" << CosMuoGen->Py_ug[i] << "  Pz_ug=" << CosMuoGen->Pz_ug[i] << std::endl;
0156       std::cout << "phi_ug=" << atan2(CosMuoGen->Px_ug[i], CosMuoGen->Pz_ug[i]) << std::endl;
0157       ;
0158     }
0159   }
0160 
0161   auto fEvt = std::make_unique<HepMC::GenEvent>();
0162 
0163   HepMC::GenVertex* Vtx_at = new HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_at,    //[mm]
0164                                                                     CosMuoGen->Vy_at,    //[mm]
0165                                                                     CosMuoGen->Vz_at,    //[mm]
0166                                                                     CosMuoGen->T0_at));  //[mm]
0167   //cout << "CosMuoGenProducer.cc: Vy_at=" << CosMuoGen->Vy_at << endl;
0168   HepMC::FourVector p_at(CosMuoGen->Px_at, CosMuoGen->Py_at, CosMuoGen->Pz_at, CosMuoGen->E_at);
0169   HepMC::GenParticle* Part_at = new HepMC::GenParticle(p_at, CosMuoGen->Id_at, 3);  //Comment mother particle in
0170   Vtx_at->add_particle_in(Part_at);
0171 
0172   //loop here in case of multi muon events (else just one iteration)
0173   for (unsigned int i = 0; i < CosMuoGen->Id_sf.size(); ++i) {
0174     HepMC::FourVector p_sf(CosMuoGen->Px_sf[i], CosMuoGen->Py_sf[i], CosMuoGen->Pz_sf[i], CosMuoGen->E_sf[i]);
0175     HepMC::GenParticle* Part_sf_in = new HepMC::GenParticle(p_sf, CosMuoGen->Id_sf[i], 3);  //Comment daughter particle
0176     Vtx_at->add_particle_out(Part_sf_in);
0177 
0178     HepMC::GenVertex* Vtx_sf = new HepMC::GenVertex(
0179         HepMC::FourVector(CosMuoGen->Vx_sf[i], CosMuoGen->Vy_sf[i], CosMuoGen->Vz_sf[i], CosMuoGen->T0_sf[i]));  //[mm]
0180     HepMC::GenParticle* Part_sf_out = new HepMC::GenParticle(p_sf, CosMuoGen->Id_sf[i], 3);  //Comment daughter particle
0181 
0182     Vtx_sf->add_particle_in(Part_sf_in);
0183     Vtx_sf->add_particle_out(Part_sf_out);
0184 
0185     fEvt->add_vertex(Vtx_sf);  //one per muon
0186 
0187     HepMC::GenVertex* Vtx_ug = new HepMC::GenVertex(
0188         HepMC::FourVector(CosMuoGen->Vx_ug[i], CosMuoGen->Vy_ug[i], CosMuoGen->Vz_ug[i], CosMuoGen->T0_ug[i]));  //[mm]
0189 
0190     HepMC::FourVector p_ug(CosMuoGen->Px_ug[i], CosMuoGen->Py_ug[i], CosMuoGen->Pz_ug[i], CosMuoGen->E_ug[i]);
0191     HepMC::GenParticle* Part_ug = new HepMC::GenParticle(p_ug, CosMuoGen->Id_ug[i], 1);  //Final state daughter particle
0192 
0193     Vtx_ug->add_particle_in(Part_sf_out);
0194     Vtx_ug->add_particle_out(Part_ug);
0195 
0196     fEvt->add_vertex(Vtx_ug);  //one per muon
0197   }
0198 
0199   fEvt->add_vertex(Vtx_at);
0200   fEvt->set_signal_process_vertex(Vtx_at);
0201 
0202   fEvt->set_event_number(e.id().event());
0203   fEvt->set_signal_process_id(13);
0204 
0205   fEvt->weights().push_back(CosMuoGen->EventWeight);  // just one event weight
0206   fEvt->weights().push_back(CosMuoGen->Trials);       // int Trials number (unweighted)
0207 
0208   if (cmVerbosity_)
0209     fEvt->print();
0210 
0211   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt.get()));
0212   e.put(std::move(genEventInfo));
0213 
0214   //This causes fEvt to be deleted
0215   std::unique_ptr<HepMCProduct> CMProduct(new HepMCProduct());
0216   CMProduct->addHepMCData(fEvt.release());
0217   e.put(std::move(CMProduct), "unsmeared");
0218 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!