ParticleGuns/src/MultiParticleInConeGunProducer.cc

0001 /*
0002  *  \author Jean-Roch Vlimant
0003  */
0004
0005 #include <ostream>
0006
0007 #include "IOMC/ParticleGuns/interface/MultiParticleInConeGunProducer.h"
0008
0009 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0010 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/ServiceRegistry/interface/Service.h"
0015 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0016
0017 #include "CLHEP/Random/RandFlat.h"
0018
0019 using namespace edm;
0020 using namespace std;
0021
0022 MultiParticleInConeGunProducer::MultiParticleInConeGunProducer(const ParameterSet& pset) : BaseFlatGunProducer(pset) {
0023   ParameterSet defpset;
0024   ParameterSet pgun_params = pset.getParameter<ParameterSet>("PGunParameters");
0025
0026   fMinPt = pgun_params.getParameter<double>("MinPt");
0027   fMaxPt = pgun_params.getParameter<double>("MaxPt");
0028
0029   fInConeIds = pgun_params.getParameter<vector<int> >("InConeID");
0030   fMinDeltaR = pgun_params.getParameter<double>("MinDeltaR");
0031   fMaxDeltaR = pgun_params.getParameter<double>("MaxDeltaR");
0032   fMinMomRatio = pgun_params.getParameter<double>("MinMomRatio");
0033   fMaxMomRatio = pgun_params.getParameter<double>("MaxMomRatio");
0034
0035   fInConeMinEta = pgun_params.getParameter<double>("InConeMinEta");
0036   fInConeMaxEta = pgun_params.getParameter<double>("InConeMaxEta");
0037   fInConeMinPhi = pgun_params.getParameter<double>("InConeMinPhi");
0038   fInConeMaxPhi = pgun_params.getParameter<double>("InConeMaxPhi");
0039   fInConeMaxTry = pgun_params.getParameter<unsigned int>("InConeMaxTry");
0040
0041   produces<HepMCProduct>("unsmeared");
0042   produces<GenEventInfoProduct>();
0043 }
0044
0045 MultiParticleInConeGunProducer::~MultiParticleInConeGunProducer() {
0046   // no need to cleanup GenEvent memory - done in HepMCProduct
0047 }
0048
0049 void MultiParticleInConeGunProducer::produce(Event& e, const EventSetup& es) {
0050   edm::Service<edm::RandomNumberGenerator> rng;
0051   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0052
0053   if (fVerbosity > 0) {
0054     cout << " MultiParticleInConeGunProducer : Begin New Event Generation" << endl;
0055   }
0056   // event loop (well, another step in it...)
0057
0058   // no need to clean up GenEvent memory - done in HepMCProduct
0059   //
0060
0061   // here re-create fEvt (memory)
0062   //
0063   fEvt = new HepMC::GenEvent();
0064
0065   // now actualy, cook up the event from PDGTable and gun parameters
0066   //
0067   // 1st, primary vertex
0068   //
0069   //HepMC::GenVertex* Vtx = new HepMC::GenVertex(CLHEP::HepLorentzVector(0.,0.,0.));
0070   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0071
0072   // loop over particles
0073   //
0074   int barcode = 1;
0075   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0076     double pt = CLHEP::RandFlat::shoot(engine, fMinPt, fMaxPt);
0077     double eta = CLHEP::RandFlat::shoot(engine, fMinEta, fMaxEta);
0078     double phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0079     int PartID = fPartIDs[ip];
0080     const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0081     double mass = PData->mass().value();
0082     double theta = 2. * atan(exp(-eta));
0083     double mom = pt / sin(theta);
0084     double px = pt * cos(phi);
0085     double py = pt * sin(phi);
0086     double pz = mom * cos(theta);
0087     double energy2 = mom * mom + mass * mass;
0088     double energy = sqrt(energy2);
0089
0090     HepMC::FourVector p(px, py, pz, energy);
0091     HepMC::GenParticle* Part = new HepMC::GenParticle(p, PartID, 1);
0092     Part->suggest_barcode(barcode);
0093     barcode++;
0094     Vtx->add_particle_out(Part);
0095
0096     if (fAddAntiParticle) {
0097     }
0098
0099     // now add the particles in the cone
0100     for (unsigned iPic = 0; iPic != fInConeIds.size(); iPic++) {
0101       unsigned int nTry = 0;
0102       while (true) {
0103         //shoot flat Deltar
0104         double dR = CLHEP::RandFlat::shoot(engine, fMinDeltaR, fMaxDeltaR);
0105         //shoot flat eta/phi mixing
0106         double alpha = CLHEP::RandFlat::shoot(engine, -3.14159265358979323846, 3.14159265358979323846);
0107         double dEta = dR * cos(alpha);
0108         double dPhi = dR * sin(alpha);
0109
0110         /*
0111        //shoot Energy of associated particle
0112        double energyIc = CLHEP::RandFlat::shoot(engine, fMinEInCone, fMaxEInCone);
0113        if (mom2Ic>0){ momIC = sqrt(mom2Ic);}
0114        */
0115         //   get kinematics
0116         double etaIc = eta + dEta;
0117         double phiIc = phi + dPhi;
0118         //put it back in -Pi:Pi if necessary. multiple time might be necessary if dR > 3
0119         const unsigned int maxL = 100;
0120         unsigned int iL = 0;
0121         while (iL++ < maxL) {
0122           if (phiIc > 3.14159265358979323846)
0123             phiIc -= 2 * 3.14159265358979323846;
0124           else if (phiIc < -3.14159265358979323846)
0125             phiIc += 2 * 3.14159265358979323846;
0126
0127           if (abs(phiIc) < 3.14159265358979323846)
0128             break;
0129         }
0130
0131         //allow to skip it if you did not run out of possible drawing
0132         if (nTry++ <= fInConeMaxTry) {
0133           //draw another one if this one is not in acceptance
0134           if (etaIc < fInConeMinEta || etaIc > fInConeMaxEta)
0135             continue;
0136           if (phiIc < fInConeMinPhi || phiIc > fInConeMaxPhi)
0137             continue;
0138         } else {
0139           if (fVerbosity > 0) {
0140             cout << " MultiParticleInConeGunProducer : could not produce a particle " << fInConeIds[iPic] << " in cone "
0141                  << fMinDeltaR << " to " << fMaxDeltaR << " within the " << fInConeMaxTry << " allowed tries." << endl;
0142           }
0143           /*         edm::LogError("MultiParticleInConeGunProducer")<< " MultiParticleInConeGunProducer : could not produce a particle "<<
0144            fInConeIds[iPic]<<" in cone "<<
0145            fMinDeltaR<<" to "<<fMaxDeltaR<<" within the "<<fInConeMaxTry<<" allowed tries.";*/
0146         }
0147         int PartIDIc = fInConeIds[iPic];
0148         const HepPDT::ParticleData* PDataIc = fPDGTable->particle(HepPDT::ParticleID(abs(PartIDIc)));
0149
0150         //shoot momentum ratio
0151         double momR = CLHEP::RandFlat::shoot(engine, fMinMomRatio, fMaxMomRatio);
0152         double massIc = PDataIc->mass().value();
0153         double momIc = momR * mom;
0154         double energyIc = sqrt(momIc * momIc + massIc * massIc);
0155
0156         double thetaIc = 2. * atan(exp(-etaIc));
0157         double pxIc = momIc * sin(thetaIc) * cos(phiIc);
0158         double pyIc = momIc * sin(thetaIc) * sin(phiIc);
0159         double pzIc = momIc * cos(thetaIc);
0160
0161         HepMC::FourVector pIc(pxIc, pyIc, pzIc, energyIc);
0162         HepMC::GenParticle* PartIc = new HepMC::GenParticle(pIc, PartIDIc, 1);
0163         PartIc->suggest_barcode(barcode);
0164         barcode++;
0165         Vtx->add_particle_out(PartIc);
0166         break;
0167       }  //try many times while not in acceptance
0168     }    //loop over the particle Ids in the cone
0169   }
0170
0171   fEvt->add_vertex(Vtx);
0172   fEvt->set_event_number(e.id().event());
0173   fEvt->set_signal_process_id(20);
0174
0175   if (fVerbosity > 0) {
0176     fEvt->print();
0177   }
0178
0179   unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0180   BProduct->addHepMCData(fEvt);
0181   e.put(std::move(BProduct), "unsmeared");
0182
0183   unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0184   e.put(std::move(genEventInfo));
0185
0186   if (fVerbosity > 0) {
0187     cout << " MultiParticleInConeGunProducer : Event Generation Done " << endl;
0188   }
0189 }