Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:11:16 

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