Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:26:40

0001 /**
0002  * file ProtonTaggerFilter.cc
0003  *
0004  * Selection of very forward protons, generated and from pileup,
0005  * in clockwise and anti-clockwise beam directions.
0006  * Access to near-beam detector acceptances.
0007  *
0008  * Author: Dmitry Zaborov
0009  */
0010 
0011 // Version: $Id: ProtonTaggerFilter.cc,v 1.3 2009/03/03 14:02:39 abdullin Exp $
0012 
0013 #include "FastSimulation/ForwardDetectors/plugins/ProtonTaggerFilter.h"
0014 
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/FileInPath.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/MakerMacros.h"
0019 #include "FWCore/Utilities/interface/Exception.h"
0020 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0021 #include "DataFormats/Common/interface/Handle.h"
0022 
0023 //#include "CLHEP/Random/RandGaussQ.h"
0024 
0025 #include <iostream>
0026 #include <list>
0027 
0028 /** read (and verify) parameters */
0029 
0030 ProtonTaggerFilter::ProtonTaggerFilter(edm::ParameterSet const& p) {
0031   std::cout << "ProtonTaggerFilter: Initializing ..." << std::endl;
0032 
0033   // ... get parameters
0034 
0035   beam1mode = p.getParameter<unsigned int>("beam1mode");
0036   beam2mode = p.getParameter<unsigned int>("beam2mode");
0037 
0038   beamCombiningMode = p.getParameter<unsigned int>("beamCombiningMode");
0039 
0040   switch (beam1mode) {
0041     case 0:
0042       std::cout << "Option chosen for beam 1: ingnore" << std::endl;
0043       break;
0044     case 1:
0045       std::cout << "Option chosen for beam 1: 420" << std::endl;
0046       break;
0047     case 2:
0048       std::cout << "Option chosen for beam 1: 220" << std::endl;
0049       break;
0050     case 3:
0051       std::cout << "Option chosen for beam 1: 420 and 220" << std::endl;
0052       break;
0053     case 4:
0054       std::cout << "Option chosen for beam 1: 420 or 220" << std::endl;
0055       break;
0056     default:
0057       throw cms::Exception("FastSimulation/ProtonTaggers") << "Error: beam1mode cannot be " << beam1mode;
0058   }
0059 
0060   switch (beam2mode) {
0061     case 0:
0062       std::cout << "Option chosen for beam 2: ingnore" << std::endl;
0063       break;
0064     case 1:
0065       std::cout << "Option chosen for beam 2: 420" << std::endl;
0066       break;
0067     case 2:
0068       std::cout << "Option chosen for beam 2: 220" << std::endl;
0069       break;
0070     case 3:
0071       std::cout << "Option chosen for beam 2: 420 and 220" << std::endl;
0072       break;
0073     case 4:
0074       std::cout << "Option chosen for beam 2: 420 or 220" << std::endl;
0075       break;
0076     default:
0077       throw cms::Exception("FastSimulation/ProtonTaggers") << "Error: beam2mode cannot be " << beam2mode;
0078   }
0079 
0080   switch (beamCombiningMode) {
0081     case 1:
0082       std::cout << "Option chosen: one proton is sufficient" << std::endl;
0083       break;
0084     case 2:
0085       std::cout << "Option chosen: two protons should be tagged" << std::endl;
0086       break;
0087     case 3:
0088       std::cout << "Option chosen: two protons should be tagged as 220+220 or 420+420" << std::endl;
0089       break;
0090     case 4:
0091       std::cout << "Option chosen: two protons should be tagged as 220+420 or 420+220" << std::endl;
0092       break;
0093     default:
0094       throw cms::Exception("FastSimulation/ProtonTaggers")
0095           << "Error: beamCombiningMode cannot be " << beamCombiningMode;
0096   }
0097 
0098   if (((beam1mode != 4) || (beam2mode != 4)) && (beamCombiningMode > 2)) {
0099     std::cerr << "Warning: beamCombiningMode = " << beamCombiningMode
0100               << " only makes sence with beam1mode = beam2mode = 4" << std::endl;
0101   }
0102 
0103   if (((beam1mode == 0) || (beam2mode == 0)) && (beamCombiningMode > 1)) {
0104     std::cerr << "Warning: You ask for 2 protons while one of the beams is set to ignore" << std::endl;
0105   }
0106 
0107   if ((beam1mode == 0) && (beam2mode == 0)) {
0108     std::cerr << "Warning: Both beams are set to ignore." << std::endl;
0109   }
0110 
0111   std::cout << "ProtonTaggerFilter: Initialized" << std::endl;
0112 }
0113 
0114 /** just empty */
0115 
0116 ProtonTaggerFilter::~ProtonTaggerFilter() { ; }
0117 
0118 /** initialize detector acceptance table */
0119 
0120 void ProtonTaggerFilter::beginJob() {
0121   std::cout << "ProtonTaggerFilter: Getting ready ..." << std::endl;
0122 
0123   edm::FileInPath myDataFile("FastSimulation/ForwardDetectors/data/acceptance_420_220.root");
0124   std::string fullPath = myDataFile.fullPath();
0125 
0126   std::cout << "Opening " << fullPath << std::endl;
0127   TFile f(fullPath.c_str());
0128 
0129   if (f.Get("description") != nullptr)
0130     std::cout << "Description found: " << f.Get("description")->GetTitle() << std::endl;
0131 
0132   std::cout << "Reading acceptance tables @#@#%@$%@$#%@%" << std::endl;
0133 
0134   helper420beam1.Init(f, "a420");
0135   helper420beam2.Init(f, "a420_b2");
0136 
0137   helper220beam1.Init(f, "a220");
0138   helper220beam2.Init(f, "a220_b2");
0139 
0140   helper420a220beam1.Init(f, "a420a220");
0141   helper420a220beam2.Init(f, "a420a220_b2");
0142 
0143   f.Close();
0144 
0145   std::cout << "ProtonTaggerFilter: Ready" << std::endl;
0146 }
0147 
0148 /** nothing to be done here */
0149 
0150 void ProtonTaggerFilter::endJob() { ; }
0151 
0152 /** Compute the detector acceptances and decide whether to filter the event */
0153 
0154 bool ProtonTaggerFilter::filter(edm::Event& iEvent, const edm::EventSetup& es) {
0155   // ... get generated event
0156 
0157   edm::Handle<edm::HepMCProduct> evtSource;
0158   iEvent.getByLabel("generatorSmeared", evtSource);
0159   const HepMC::GenEvent* genEvent = evtSource->GetEvent();
0160 
0161   //std::cout << "event contains " << genEvent->particles_size() << " particles " << std::endl;
0162   if (genEvent->particles_empty()) {
0163     std::cout << "empty source event" << std::endl;
0164     return false;
0165   }
0166 
0167   // ... get pileup event
0168 
0169   edm::Handle<edm::HepMCProduct> pileUpSource;
0170   const HepMC::GenEvent* pileUpEvent = nullptr;
0171   bool isPileUp = true;
0172 
0173   bool isProduct = iEvent.getByLabel("famosPileUp", "PileUpEvents", pileUpSource);
0174 
0175   if (isProduct) {
0176     pileUpEvent = pileUpSource->GetEvent();
0177     //std::cout << "got pileup" << std::endl;
0178   } else {
0179     isPileUp = false;
0180     //std::cout << "no pileup in the event" << std::endl;
0181   }
0182 
0183   //if (isPileUp) std::cout << "event contains " << pileUpEvent->particles_size() << " pileup particles " << std::endl;
0184 
0185   // ... some constants
0186 
0187   const double mp = 0.938272029;      // just a proton mass
0188   const double Ebeam = 7000.0;        // beam energy - would be better to read from parameters
0189   const float pzCut = 2500;           // ignore particles with less than |Pz| < pzCut
0190   const float acceptThreshold = 0.5;  // proton will be "accepted" if the value of acceptance is above this threshold
0191 
0192   // ... loop over particles, find the most energetic proton in either direction
0193 
0194   std::list<HepMC::GenParticle*> veryForwardParicles;
0195 
0196   for (HepMC::GenEvent::particle_const_iterator piter = genEvent->particles_begin(); piter != genEvent->particles_end();
0197        ++piter) {
0198     HepMC::GenParticle* p = *piter;
0199 
0200     float pz = p->momentum().pz();
0201     if (((pz > pzCut) || (pz < -pzCut)) && ((p->status() == 0) || (p->status() == 1))) {
0202       veryForwardParicles.push_back(p);
0203       //std::cout << "pdgid: " << p->pdg_id() << " status: " << p->status() << std::endl;
0204     }
0205   }
0206 
0207   //std::cout << "# generated forward particles  : " << veryForwardParicles.size() << std::endl;
0208 
0209   if (isPileUp) {
0210     //std::cout << "Adding pileup " << std::endl;
0211 
0212     for (HepMC::GenEvent::particle_const_iterator piter = pileUpEvent->particles_begin();
0213          piter != pileUpEvent->particles_end();
0214          ++piter) {
0215       HepMC::GenParticle* p = *piter;
0216 
0217       float pz = p->momentum().pz();
0218       if (((pz > pzCut) || (pz < -pzCut)) && ((p->status() == 0) || (p->status() == 1))) {
0219         veryForwardParicles.push_back(p);
0220         //std::cout << "pdgid: " << p->pdg_id() << " status: " << p->status() << std::endl;
0221       }
0222     }
0223   }
0224 
0225   //std::cout << "# forward particles to be tried: " << veryForwardParicles.size() << std::endl;
0226 
0227   // ... return false if no forward protons found
0228 
0229   if (veryForwardParicles.empty())
0230     return false;
0231 
0232   // ... set all acceptances to zero
0233 
0234   float acc420b1, acc220b1, acc420and220b1, acc420or220b1;  // beam 1 (clockwise)
0235   float acc420b2, acc220b2, acc420and220b2, acc420or220b2;  // beam 2 (anti-clockwise)
0236 
0237   acc420b1 = acc220b1 = acc420and220b1 = acc420or220b1 = 0;
0238   acc420b2 = acc220b2 = acc420and220b2 = acc420or220b2 = 0;
0239 
0240   int nP1at220m = 0;
0241   int nP1at420m = 0;
0242 
0243   int nP2at220m = 0;
0244   int nP2at420m = 0;
0245 
0246   // ... loop over (pre-selected) forward particles
0247 
0248   for (std::list<HepMC::GenParticle*>::const_iterator part = veryForwardParicles.begin();
0249        part != veryForwardParicles.end();
0250        part++) {
0251     HepMC::GenParticle* p = *part;
0252 
0253     float pz = p->momentum().pz();
0254     float pt = p->momentum().perp();
0255     float phi = p->momentum().phi();
0256     ;
0257 
0258     if ((pz > Ebeam) || (pz < -Ebeam))
0259       continue;
0260 
0261     // ... compute kinimatical variable
0262 
0263     float xi = 1.0;  // fractional momentum loss
0264     if (pz > 0)
0265       xi -= pz / Ebeam;
0266     else
0267       xi += pz / Ebeam;
0268 
0269     double t = (-pt * pt - mp * mp * xi * xi) / (1 - xi);  // "t"
0270 
0271     //std::cout << " pdg_id: "  << p->pdg_id() << " eta: " << p->momentum().eta() << " e: "
0272     //          <<  p->momentum().e() << std::endl;
0273     //std::cout << "pz: " << pz << " pt: " <<  pt << " xi: " << xi
0274     //          << " t: " << t << " phi: " << phi << std::endl;
0275 
0276     if (xi < 0.0)
0277       xi = -10.0;
0278     if (xi > 1.0)
0279       xi = 10.0;
0280 
0281     //float rnd1 = RandGauss::shoot(0.,1.);
0282     //float rnd2 = RandGauss::shoot(0.,1.);
0283 
0284     // ... get acceptance from tables: beam 1 (if it is not ignored)
0285 
0286     if ((pz > 0) && (beam1mode != 0)) {
0287       acc420b1 = helper420beam1.GetAcceptance(t, xi, phi);
0288       acc220b1 = helper220beam1.GetAcceptance(t, xi, phi);
0289       acc420and220b1 = helper420a220beam1.GetAcceptance(t, xi, phi);
0290 
0291       acc420or220b1 = acc420b1 + acc220b1 - acc420and220b1;
0292 
0293       //std::cout << "+acc420b1: " << acc420b1 << " acc220b1: " << acc220b1 << " acc420and220b1: " << acc420and220b1 << " acc420or220b1: " << acc420or220b1  << std::endl;
0294 
0295       bool res420and220 = (acc420and220b1 > acceptThreshold);
0296       bool res420 = (acc420b1 > acceptThreshold);
0297       bool res220 = (acc220b1 > acceptThreshold);
0298 
0299       if (res420and220) {
0300         nP1at220m++;
0301         nP1at420m++;
0302       } else if (res420)
0303         nP1at420m++;
0304       else if (res220)
0305         nP1at220m++;
0306 
0307       if ((p->pdg_id() != 2212) && (res220 || res420 || res420and220)) {
0308         std::cout << " !!! P got proton 1 at 420 m: pz = " << pz << std::endl;
0309         if (res220)
0310           std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m along beam 1, pz = " << pz << std::endl;
0311         if (res420)
0312           std::cout << "got a particle with pid" << p->pdg_id() << " at 420 m along beam 1, pz = " << pz << std::endl;
0313         if (res420and220)
0314           std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m & 420 m  along beam 1, pz = " << pz
0315                     << std::endl;
0316       }
0317     }
0318 
0319     // ... the same for beam 2
0320 
0321     if ((pz < 0) && (beam2mode != 0)) {
0322       acc420b2 = helper420beam2.GetAcceptance(t, xi, phi);
0323       acc220b2 = helper220beam2.GetAcceptance(t, xi, phi);
0324       acc420and220b2 = helper420a220beam2.GetAcceptance(t, xi, phi);
0325 
0326       acc420or220b2 = acc420b2 + acc220b2 - acc420and220b2;
0327 
0328       //std::cout << "+acc420b2: " << acc420b2 << " acc220b2: " << acc220b2 << " acc420and220b2: " << acc420and220b2 << " acc420or220b2: " << acc420or220b2 << std::endl;
0329 
0330       bool res420and220 = (acc420and220b2 > acceptThreshold);
0331       bool res420 = (acc420b2 > acceptThreshold);
0332       bool res220 = (acc220b2 > acceptThreshold);
0333 
0334       if (res420and220) {
0335         nP2at220m++;
0336         nP2at420m++;
0337       } else if (res420)
0338         nP2at420m++;
0339       else if (res220)
0340         nP2at220m++;
0341 
0342       if ((p->pdg_id() != 2212) && (res220 || res420 || res420and220)) {
0343         std::cout << " !!! P got proton 1 at 420 m: pz = " << pz << std::endl;
0344         if (res220)
0345           std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m along beam 2, pz = " << pz << std::endl;
0346         if (res420)
0347           std::cout << "got a particle with pid" << p->pdg_id() << " at 420 m along beam 2, pz = " << pz << std::endl;
0348         if (res420and220)
0349           std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m & 420 m along beam 2, pz = " << pz
0350                     << std::endl;
0351       }
0352     }
0353   }
0354 
0355   // ... boolean result for each detector
0356 
0357   bool p1at220m = (nP1at220m > 0) ? true : false;
0358   bool p1at420m = (nP1at420m > 0) ? true : false;
0359   bool p2at220m = (nP2at220m > 0) ? true : false;
0360   bool p2at420m = (nP2at420m > 0) ? true : false;
0361 
0362   if ((nP1at220m > 1) && (beam1mode != 1))
0363     std::cout << " !!! " << nP1at220m << " proton(s) from beam 1 at 220 m" << std::endl;
0364   if ((nP1at420m > 1) && (beam1mode != 2))
0365     std::cout << " !!! " << nP1at420m << " proton(s) from beam 1 at 420 m" << std::endl;
0366   if ((nP2at220m > 1) && (beam2mode != 1))
0367     std::cout << " !!! " << nP2at220m << " proton(s) from beam 2 at 220 m" << std::endl;
0368   if ((nP2at420m > 1) && (beam2mode != 2))
0369     std::cout << " !!! " << nP2at420m << " proton(s) from beam 2 at 420 m" << std::endl;
0370 
0371   // ... make a decision based on requested filter configuration
0372 
0373   bool p1accepted = false;
0374   bool p2accepted = false;
0375 
0376   if ((beam1mode == 1) && p1at420m)
0377     p1accepted = true;
0378   if ((beam1mode == 2) && p1at220m)
0379     p1accepted = true;
0380   if ((beam1mode == 3) && p1at220m && p1at420m)
0381     p1accepted = true;
0382   if ((beam1mode == 4) && (p1at220m || p1at420m))
0383     p1accepted = true;
0384 
0385   if ((beam2mode == 1) && p2at420m)
0386     p2accepted = true;
0387   if ((beam2mode == 2) && p2at220m)
0388     p2accepted = true;
0389   if ((beam2mode == 3) && p2at220m && p2at420m)
0390     p2accepted = true;
0391   if ((beam2mode == 4) && (p2at220m || p2at420m))
0392     p2accepted = true;
0393 
0394   //if (p1accepted) std::cout << "proton 1 accepted" << std::endl;
0395   //if (p2accepted) std::cout << "proton 2 accepted" << std::endl;
0396 
0397   switch (beamCombiningMode) {
0398     case 1:  // ... either of two protons
0399       if (p1accepted || p2accepted)
0400         return true;
0401       else
0402         return false;
0403 
0404     case 2:  // ... both protons
0405       if (p1accepted && p2accepted)
0406         return true;
0407       else
0408         return false;
0409 
0410     case 3:  // ... 220+220 or 420+420
0411       if ((p1at220m && p2at220m) || (p1at420m && p2at420m))
0412         return true;
0413       else
0414         return false;
0415 
0416     case 4:  // ... 220+420 or 420+220
0417       if ((p1at220m && p2at420m) || (p1at420m && p2at220m))
0418         return true;
0419       else
0420         return false;
0421   }
0422 
0423   return false;
0424 }
0425 
0426 // ... define CMSSW module
0427 
0428 DEFINE_FWK_MODULE(ProtonTaggerFilter);