ParticleGuns/src/RandomtXiGunProducer.cc

0001 /*
0002  *  $Date: 2011/12/19 23:10:40 $
0003  *  $Revision: 1.4 $
0004  *  \author Luiz Mundim
0005  */
0006
0007 #include <ostream>
0008
0009 #include "IOMC/ParticleGuns/interface/RandomtXiGunProducer.h"
0010
0011 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0012 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0013
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0019
0020 #include <CLHEP/Random/RandFlat.h>
0021
0022 using namespace edm;
0023 using namespace std;
0024
0025 RandomtXiGunProducer::RandomtXiGunProducer(const edm::ParameterSet& pset) : BaseRandomtXiGunProducer(pset) {
0026   ParameterSet defpset;
0027   edm::ParameterSet pgun_params = pset.getParameter<edm::ParameterSet>("PGunParameters");
0028
0029   fMint = pgun_params.getParameter<double>("Mint");
0030   fMaxt = pgun_params.getParameter<double>("Maxt");
0031   fMinXi = pgun_params.getParameter<double>("MinXi");
0032   fMaxXi = pgun_params.getParameter<double>("MaxXi");
0033   fLog_t = pgun_params.getUntrackedParameter<bool>("Log_t", false);
0034
0035   produces<HepMCProduct>("unsmeared");
0036   produces<GenEventInfoProduct>();
0037 }
0038
0039 RandomtXiGunProducer::~RandomtXiGunProducer() {
0040   // no need to cleanup GenEvent memory - done in HepMCProduct
0041 }
0042
0043 void RandomtXiGunProducer::produce(edm::Event& e, const edm::EventSetup& es) {
0044   edm::Service<edm::RandomNumberGenerator> rng;
0045   CLHEP::HepRandomEngine* engine = &rng->getEngine(e.streamID());
0046
0047   if (fVerbosity > 0) {
0048     edm::LogInfo("RandomtXiGunProducer") << "Begin New Event Generation\n";
0049   }
0050   // event loop (well, another step in it...)
0051
0052   // no need to clean up GenEvent memory - done in HepMCProduct
0053   //
0054
0055   // here re-create fEvt (memory)
0056   //
0057   fEvt = new HepMC::GenEvent();
0058
0059   // now actualy, cook up the event from PDGTable and gun parameters
0060   //
0061   // 1st, primary vertex
0062   //
0063   HepMC::GenVertex* Vtx = new HepMC::GenVertex(HepMC::FourVector(0., 0., 0.));
0064
0065   // loop over particles
0066   //
0067   int barcode = 1;
0068   for (unsigned int ip = 0; ip < fPartIDs.size(); ++ip) {
0069     int PartID = fPartIDs[ip];
0070     //  t = -2*P*P'*(1-cos(theta)) -> t/(2*P*P')+1=cos(theta)
0071     // xi = 1 - P'/P  --> P'= (1-xi)*P
0072     //
0073     PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID)));
0074     if (!PData)
0075       exit(1);
0076     double t = 0;
0077     double Xi = 0;
0078     double phi = 0;
0079     if (fFireForward) {
0080       while (true) {
0081         Xi = CLHEP::RandFlat::shoot(engine, fMinXi, fMaxXi);
0082         double min_t = std::max(fMint, Minimum_t(Xi));
0083         double max_t = fMaxt;
0084         if (min_t > max_t) {
0085           edm::LogInfo("RandomtXiGunProducer") << "WARNING: t limits redefined (unphysical values for given xi).\n";
0086           max_t = min_t;
0087         }
0088         t = (fLog_t) ? pow(CLHEP::RandFlat::shoot(engine, log10(min_t), log10(max_t)), 10)
0089                      : CLHEP::RandFlat::shoot(engine, min_t, max_t);
0090         break;
0091       }
0092       phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0093       HepMC::GenParticle* Part = new HepMC::GenParticle(make_particle(t, Xi, phi, PartID, 1), PartID, 1);
0094       Part->suggest_barcode(barcode);
0095       barcode++;
0096       Vtx->add_particle_out(Part);
0097     }
0098     if (fFireBackward) {
0099       while (true) {
0100         Xi = CLHEP::RandFlat::shoot(engine, fMinXi, fMaxXi);
0101         double min_t = std::max(fMint, Minimum_t(Xi));
0102         double max_t = fMaxt;
0103         if (min_t > max_t) {
0104           edm::LogInfo("RandomtXiGunProducer")
0105               << "WARNING: t limits redefined (unphysical values for given xi)." << endl;
0106           max_t = min_t;
0107         }
0108         t = (fLog_t) ? pow(CLHEP::RandFlat::shoot(engine, log10(min_t), log10(max_t)), 10)
0109                      : CLHEP::RandFlat::shoot(engine, min_t, max_t);
0110         break;
0111       }
0112       phi = CLHEP::RandFlat::shoot(engine, fMinPhi, fMaxPhi);
0113       HepMC::GenParticle* Part2 = new HepMC::GenParticle(make_particle(t, Xi, phi, PartID, -1), PartID, 1);
0114       Part2->suggest_barcode(barcode);
0115       barcode++;
0116       Vtx->add_particle_out(Part2);
0117     }
0118   }
0119
0120   fEvt->add_vertex(Vtx);
0121   fEvt->set_event_number(e.id().event());
0122   fEvt->set_signal_process_id(20);
0123
0124   if (fVerbosity > 0) {
0125     fEvt->print();
0126   }
0127
0128   std::unique_ptr<HepMCProduct> BProduct(new HepMCProduct());
0129   BProduct->addHepMCData(fEvt);
0130   e.put(std::move(BProduct), "unsmeared");
0131
0132   std::unique_ptr<GenEventInfoProduct> genEventInfo(new GenEventInfoProduct(fEvt));
0133   e.put(std::move(genEventInfo));
0134
0135   if (fVerbosity > 0) {
0136     // for testing purpose only
0137     // fEvt->print() ; // prints empty info after it's made into edm::Event
0138     edm::LogInfo("RandomtXiGunProducer") << " Event Generation Done \n";
0139   }
0140 }
0141 HepMC::FourVector RandomtXiGunProducer::make_particle(double t, double Xi, double phi, int PartID, int direction) {
0142   double mass = PData->mass().value();
0143   double fpMom = sqrt(fpEnergy * fpEnergy - mass * mass);  // momentum of beam proton
0144   double sEnergy = (1.0 - Xi) * fpEnergy;                  // energy of scattered particle
0145   double sMom = sqrt(sEnergy * sEnergy - mass * mass);     // momentum of scattered particle
0146   double min_t = -2. * (fpMom * sMom - fpEnergy * sEnergy + mass * mass);
0147   if (t < min_t)
0148     t = min_t;  // protect against kinemactically forbiden region
0149   long double theta = acos((-t / 2. - mass * mass + fpEnergy * sEnergy) / (sMom * fpMom));  // use t = -t
0150
0151   if (direction < 1)
0152     theta = acos(-1.) - theta;
0153
0154   double px = sMom * cos(phi) * sin(theta);
0155   double py = sMom * sin(phi) * sin(theta);
0156   double pz = sMom * cos(theta);  // the direction is already set by the theta angle
0157   if (fVerbosity > 0)
0158     edm::LogInfo("RandomXiGunProducer")
0159         << "-----------------------------------------------------------------------------------------------------\n"
0160         << "Produced a proton with  phi : " << phi << " theta: " << theta << " t: " << t << " Xi: " << Xi << "\n"
0161         << "                         Px : " << px << " Py : " << py << " Pz : " << pz << "\n"
0162         << "                   direction: " << direction << "\n"
0163         << "-----------------------------------------------------------------------------------------------------"
0164         << std::endl;
0165
0166   return HepMC::FourVector(px, py, pz, sEnergy);
0167 }
0168 //#include "FWCore/Framework/interface/MakerMacros.h"
0169 //DEFINE_FWK_MODULE(RandomtXiGunProducer);