CosMuoGenProducer.cc

Go to the documentation of this file.

#include "FWCore/Framework/interface/Run.h"
#include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"

#include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"

#include "GeneratorInterface/CosmicMuonGenerator/interface/CosMuoGenProducer.h"

edm::CosMuoGenProducer::CosMuoGenProducer(const ParameterSet& pset)
    :  //RanS(pset.getParameter<int>("RanSeed", 123456)), //get seed now from Framework
      MinP(pset.getParameter<double>("MinP")),
      MinP_CMS(pset.getParameter<double>("MinP_CMS")),
      MaxP(pset.getParameter<double>("MaxP")),
      MinT(pset.getParameter<double>("MinTheta")),
      MaxT(pset.getParameter<double>("MaxTheta")),
      MinPh(pset.getParameter<double>("MinPhi")),
      MaxPh(pset.getParameter<double>("MaxPhi")),
      MinS(pset.getParameter<double>("MinT0")),
      MaxS(pset.getParameter<double>("MaxT0")),
      ELSF(pset.getParameter<double>("ElossScaleFactor")),
      RTarget(pset.getParameter<double>("RadiusOfTarget")),
      ZTarget(pset.getParameter<double>("ZDistOfTarget")),
      ZCTarget(pset.getParameter<double>("ZCentrOfTarget")),
      TrackerOnly(pset.getParameter<bool>("TrackerOnly")),
      MultiMuon(pset.getParameter<bool>("MultiMuon")),
      MultiMuonFileName(pset.getParameter<std::string>("MultiMuonFileName")),
      MultiMuonFileFirstEvent(pset.getParameter<int>("MultiMuonFileFirstEvent")),
      MultiMuonNmin(pset.getParameter<int>("MultiMuonNmin")),
      TIFOnly_constant(pset.getParameter<bool>("TIFOnly_constant")),
      TIFOnly_linear(pset.getParameter<bool>("TIFOnly_linear")),
      MTCCHalf(pset.getParameter<bool>("MTCCHalf")),
      PlugVtx(pset.getParameter<double>("PlugVx")),
      PlugVtz(pset.getParameter<double>("PlugVz")),
      VarRhoAir(pset.getParameter<double>("RhoAir")),
      VarRhoWall(pset.getParameter<double>("RhoWall")),
      VarRhoRock(pset.getParameter<double>("RhoRock")),
      VarRhoClay(pset.getParameter<double>("RhoClay")),
      VarRhoPlug(pset.getParameter<double>("RhoPlug")),
      ClayLayerWidth(pset.getParameter<double>("ClayWidth")),
      MinEn(pset.getParameter<double>("MinEnu")),
      MaxEn(pset.getParameter<double>("MaxEnu")),
      NuPrdAlt(pset.getParameter<double>("NuProdAlt")),
      AllMu(pset.getParameter<bool>("AcptAllMu")),
      extCrossSect(pset.getUntrackedParameter<double>("crossSection", -1.)),
      extFilterEff(pset.getUntrackedParameter<double>("filterEfficiency", -1.)),
      cmVerbosity_(pset.getParameter<bool>("Verbosity")),
      isInitialized_(false) {
  //if not specified (i.e. negative) then use MinP also for MinP_CMS
  if (MinP_CMS < 0)
    MinP_CMS = MinP;

  // set up the generator
  CosMuoGen = std::make_unique<CosmicMuonGenerator>();
  // Begin JMM change
  //  CosMuoGen->setNumberOfEvents(numberEventsInRun());
  CosMuoGen->setNumberOfEvents(999999999);
  // End of JMM change
  CosMuoGen->setRanSeed(RanS);
  CosMuoGen->setMinP(MinP);
  CosMuoGen->setMinP_CMS(MinP_CMS);
  CosMuoGen->setMaxP(MaxP);
  CosMuoGen->setMinTheta(MinT);
  CosMuoGen->setMaxTheta(MaxT);
  CosMuoGen->setMinPhi(MinPh);
  CosMuoGen->setMaxPhi(MaxPh);
  CosMuoGen->setMinT0(MinS);
  CosMuoGen->setMaxT0(MaxS);
  CosMuoGen->setElossScaleFactor(ELSF);
  CosMuoGen->setRadiusOfTarget(RTarget);
  CosMuoGen->setZDistOfTarget(ZTarget);
  CosMuoGen->setZCentrOfTarget(ZCTarget);
  CosMuoGen->setTrackerOnly(TrackerOnly);
  CosMuoGen->setMultiMuon(MultiMuon);
  CosMuoGen->setMultiMuonFileName(MultiMuonFileName);
  CosMuoGen->setMultiMuonFileFirstEvent(MultiMuonFileFirstEvent);
  CosMuoGen->setMultiMuonNmin(MultiMuonNmin);
  CosMuoGen->setTIFOnly_constant(TIFOnly_constant);
  CosMuoGen->setTIFOnly_linear(TIFOnly_linear);
  CosMuoGen->setMTCCHalf(MTCCHalf);
  CosMuoGen->setPlugVx(PlugVtx);
  CosMuoGen->setPlugVz(PlugVtz);
  CosMuoGen->setRhoAir(VarRhoAir);
  CosMuoGen->setRhoWall(VarRhoWall);
  CosMuoGen->setRhoRock(VarRhoRock);
  CosMuoGen->setRhoClay(VarRhoClay);
  CosMuoGen->setRhoPlug(VarRhoPlug);
  CosMuoGen->setClayWidth(ClayLayerWidth);
  CosMuoGen->setMinEnu(MinEn);
  CosMuoGen->setMaxEnu(MaxEn);
  CosMuoGen->setNuProdAlt(NuPrdAlt);
  CosMuoGen->setAcptAllMu(AllMu);
  produces<HepMCProduct>("unsmeared");
  produces<GenEventInfoProduct>();
  produces<GenRunInfoProduct, edm::Transition::EndRun>();
}

edm::CosMuoGenProducer::~CosMuoGenProducer() { clear(); }

void edm::CosMuoGenProducer::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const&) {
  if (!isInitialized_) {
    isInitialized_ = true;
    RandomEngineSentry<CosmicMuonGenerator> randomEngineSentry(CosMuoGen.get(), lumi.index());
    CosMuoGen->initialize(randomEngineSentry.randomEngine());
  }
}

void edm::CosMuoGenProducer::endRunProduce(Run& run, const EventSetup& es) {
  std::unique_ptr<GenRunInfoProduct> genRunInfo(new GenRunInfoProduct());

  double cs = CosMuoGen->getRate();  // flux in Hz, not s^-1m^-2
  if (MultiMuon)
    genRunInfo->setInternalXSec(0.);
  else
    genRunInfo->setInternalXSec(cs);
  genRunInfo->setExternalXSecLO(extCrossSect);
  genRunInfo->setFilterEfficiency(extFilterEff);

  run.put(std::move(genRunInfo));

  CosMuoGen->terminate();
}

void edm::CosMuoGenProducer::clear() {}

void edm::CosMuoGenProducer::produce(Event& e, const edm::EventSetup& es) {
  RandomEngineSentry<CosmicMuonGenerator> randomEngineSentry(CosMuoGen.get(), e.streamID());

  // generate event
  if (!MultiMuon) {
    CosMuoGen->nextEvent();
  } else {
    bool success = CosMuoGen->nextMultiEvent();
    if (!success)
      std::cout << "CosMuoGenProducer.cc: CosMuoGen->nextMultiEvent() failed!" << std::endl;
  }

  if (Debug) {
    std::cout << "CosMuoGenProducer.cc: CosMuoGen->EventWeight=" << CosMuoGen->EventWeight
              << "  CosMuoGen: Nmuons=" << CosMuoGen->Id_sf.size() << std::endl;
    std::cout << "CosMuoGen->Id_at=" << CosMuoGen->Id_at << "  CosMuoGen->Vx_at=" << CosMuoGen->Vx_at
              << "  CosMuoGen->Vy_at=" << CosMuoGen->Vy_at << "  CosMuoGen->Vz_at=" << CosMuoGen->Vz_at
              << "  CosMuoGen->T0_at=" << CosMuoGen->T0_at << std::endl;
    std::cout << "  Px=" << CosMuoGen->Px_at << "  Py=" << CosMuoGen->Py_at << "  Pz=" << CosMuoGen->Pz_at << std::endl;
    for (unsigned int i = 0; i < CosMuoGen->Id_sf.size(); ++i) {
      std::cout << "Id_sf[" << i << "]=" << CosMuoGen->Id_sf[i] << "  Vx_sf[" << i << "]=" << CosMuoGen->Vx_sf[i]
                << "  Vy_sf=" << CosMuoGen->Vy_sf[i] << "  Vz_sf=" << CosMuoGen->Vz_sf[i]
                << "  T0_sf=" << CosMuoGen->T0_sf[i] << "  Px_sf=" << CosMuoGen->Px_sf[i]
                << "  Py_sf=" << CosMuoGen->Py_sf[i] << "  Pz_sf=" << CosMuoGen->Pz_sf[i] << std::endl;
      std::cout << "phi_sf=" << atan2(CosMuoGen->Px_sf[i], CosMuoGen->Pz_sf[i]) << std::endl;
      std::cout << "Id_ug[" << i << "]=" << CosMuoGen->Id_ug[i] << "  Vx_ug[" << i << "]=" << CosMuoGen->Vx_ug[i]
                << "  Vy_ug=" << CosMuoGen->Vy_ug[i] << "  Vz_ug=" << CosMuoGen->Vz_ug[i]
                << "  T0_ug=" << CosMuoGen->T0_ug[i] << "  Px_ug=" << CosMuoGen->Px_ug[i]
                << "  Py_ug=" << CosMuoGen->Py_ug[i] << "  Pz_ug=" << CosMuoGen->Pz_ug[i] << std::endl;
      std::cout << "phi_ug=" << atan2(CosMuoGen->Px_ug[i], CosMuoGen->Pz_ug[i]) << std::endl;
      ;
    }
  }

  auto fEvt = std::make_unique<HepMC::GenEvent>();

  HepMC::GenVertex* Vtx_at = new HepMC::GenVertex(HepMC::FourVector(CosMuoGen->Vx_at,    //[mm]
                                                                    CosMuoGen->Vy_at,    //[mm]
                                                                    CosMuoGen->Vz_at,    //[mm]
                                                                    CosMuoGen->T0_at));  //[mm]
  //cout << "CosMuoGenProducer.cc: Vy_at=" << CosMuoGen->Vy_at << endl;
  HepMC::FourVector p_at(CosMuoGen->Px_at, CosMuoGen->Py_at, CosMuoGen->Pz_at, CosMuoGen->E_at);
  HepMC::GenParticle* Part_at = new HepMC::GenParticle(p_at, CosMuoGen->Id_at, 3);  //Comment mother particle in
  Vtx_at->add_particle_in(Part_at);

  //loop here in case of multi muon events (else just one iteration)
  for (unsigned int i = 0; i < CosMuoGen->Id_sf.size(); ++i) {
    HepMC::FourVector p_sf(CosMuoGen->Px_sf[i], CosMuoGen->Py_sf[i], CosMuoGen->Pz_sf[i], CosMuoGen->E_sf[i]);
    HepMC::GenParticle* Part_sf_in = new HepMC::GenParticle(p_sf, CosMuoGen->Id_sf[i], 3);  //Comment daughter particle
    Vtx_at->add_particle_out(Part_sf_in);

    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]
    HepMC::GenParticle* Part_sf_out = new HepMC::GenParticle(p_sf, CosMuoGen->Id_sf[i], 3);  //Comment daughter particle

    Vtx_sf->add_particle_in(Part_sf_in);
    Vtx_sf->add_particle_out(Part_sf_out);

    fEvt->add_vertex(Vtx_sf);  //one per muon

    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]

    HepMC::FourVector p_ug(CosMuoGen->Px_ug[i], CosMuoGen->Py_ug[i], CosMuoGen->Pz_ug[i], CosMuoGen->E_ug[i]);
    HepMC::GenParticle* Part_ug = new HepMC::GenParticle(p_ug, CosMuoGen->Id_ug[i], 1);  //Final state daughter particle

    Vtx_ug->add_particle_in(Part_sf_out);
    Vtx_ug->add_particle_out(Part_ug);

    fEvt->add_vertex(Vtx_ug);  //one per muon
  }

  fEvt->add_vertex(Vtx_at);
  fEvt->set_signal_process_vertex(Vtx_at);

  fEvt->set_event_number(e.id().event());
  fEvt->set_signal_process_id(13);

  fEvt->weights().push_back(CosMuoGen->EventWeight);  // just one event weight
  fEvt->weights().push_back(CosMuoGen->Trials);       // int Trials number (unweighted)

  if (cmVerbosity_)
    fEvt->print();

  std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt.get()));
  e.put(std::move(genEventInfo));

  //This causes fEvt to be deleted
  std::unique_ptr<HepMCProduct> CMProduct(new HepMCProduct());
  CMProduct->addHepMCData(fEvt.release());
  e.put(std::move(CMProduct), "unsmeared");
}