Reference Manual

Go to the documentation of this file.
00001 
00011 // Version: $Id: ProtonTaggerFilter.cc,v 1.4 2009/04/28 12:01:09 dzaborov Exp $
00012 
00013 #include "FastSimulation/ForwardDetectors/plugins/ProtonTaggerFilter.h"
00014 
00015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00016 #include "FWCore/ParameterSet/interface/FileInPath.h"
00017 #include "FWCore/Framework/interface/Event.h"
00018 #include "FWCore/Framework/interface/MakerMacros.h"
00019 #include "FWCore/Utilities/interface/Exception.h"
00020 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00021 #include "DataFormats/Common/interface/Handle.h"
00022 
00023 
00024 //#include "CLHEP/Random/RandGaussQ.h"
00025 
00026 
00027 #include <iostream>
00028 #include <list>
00029 
00030 
00033 ProtonTaggerFilter::ProtonTaggerFilter(edm::ParameterSet const & p)      
00034 {
00035   std::cout << "ProtonTaggerFilter: Initializing ..." << std::endl;
00036   
00037   // ... get parameters
00038 
00039   beam1mode = p.getParameter<unsigned int>("beam1mode");
00040   beam2mode = p.getParameter<unsigned int>("beam2mode");
00041 
00042   beamCombiningMode = p.getParameter<unsigned int>("beamCombiningMode");
00043   
00044   switch (beam1mode)
00045   {
00046     case 0: std::cout << "Option chosen for beam 1: ingnore" << std::endl; break;
00047     case 1: std::cout << "Option chosen for beam 1: 420" << std::endl; break;
00048     case 2: std::cout << "Option chosen for beam 1: 220" << std::endl; break;
00049     case 3: std::cout << "Option chosen for beam 1: 420 and 220" << std::endl; break;
00050     case 4: std::cout << "Option chosen for beam 1: 420 or 220" << std::endl; break;
00051     default: throw cms::Exception("FastSimulation/ProtonTaggers") 
00052                            << "Error: beam1mode cannot be " << beam1mode;
00053   }
00054   
00055   switch (beam2mode)
00056   {
00057     case 0: std::cout << "Option chosen for beam 2: ingnore" << std::endl; break;
00058     case 1: std::cout << "Option chosen for beam 2: 420" << std::endl; break;
00059     case 2: std::cout << "Option chosen for beam 2: 220" << std::endl; break;
00060     case 3: std::cout << "Option chosen for beam 2: 420 and 220" << std::endl; break;
00061     case 4: std::cout << "Option chosen for beam 2: 420 or 220" << std::endl; break;
00062     default: throw cms::Exception("FastSimulation/ProtonTaggers")
00063                            << "Error: beam2mode cannot be " << beam2mode;
00064   }
00065 
00066   switch (beamCombiningMode)
00067   {
00068     case 1: std::cout << "Option chosen: one proton is sufficient" << std::endl; break;
00069     case 2: std::cout << "Option chosen: two protons should be tagged" << std::endl; break;
00070     case 3: std::cout << "Option chosen: two protons should be tagged as 220+220 or 420+420" << std::endl; break;
00071     case 4: std::cout << "Option chosen: two protons should be tagged as 220+420 or 420+220" << std::endl; break;
00072     default: throw cms::Exception("FastSimulation/ProtonTaggers") 
00073                              << "Error: beamCombiningMode cannot be " << beamCombiningMode;
00074   }
00075 
00076   if (((beam1mode != 4) || (beam2mode != 4)) && (beamCombiningMode > 2)) {
00077     std::cerr << "Warning: beamCombiningMode = " << beamCombiningMode 
00078               << " only makes sence with beam1mode = beam2mode = 4" << std::endl;
00079   }
00080 
00081   if (((beam1mode == 0) || (beam2mode == 0)) && (beamCombiningMode > 1)) {
00082     std::cerr << "Warning: You ask for 2 protons while one of the beams is set to ignore" << std::endl;
00083   }
00084 
00085   if ((beam1mode == 0) && (beam2mode == 0)) {
00086     std::cerr << "Warning: Both beams are set to ignore." << std::endl;
00087   }
00088 
00089   std::cout << "ProtonTaggerFilter: Initialized" << std::endl;
00090 }
00091 
00092 
00095 ProtonTaggerFilter::~ProtonTaggerFilter() {;}
00096 
00097 
00100 void ProtonTaggerFilter::beginJob()
00101 {
00102   std::cout << "ProtonTaggerFilter: Getting ready ..." << std::endl;
00103 
00104   edm::FileInPath myDataFile("FastSimulation/ForwardDetectors/data/acceptance_420_220.root");
00105   std::string fullPath = myDataFile.fullPath();
00106 
00107   std::cout << "Opening " << fullPath << std::endl;
00108   TFile f(fullPath.c_str());
00109 
00110   if (f.Get("description") != NULL)
00111       std::cout << "Description found: " << f.Get("description")->GetTitle() << std::endl;
00112   
00113   std::cout << "Reading acceptance tables @#@#%@$%@$#%@%" << std::endl;
00114   
00115   helper420beam1.Init(f, "a420");
00116   helper420beam2.Init(f, "a420_b2");
00117 
00118   helper220beam1.Init(f, "a220");
00119   helper220beam2.Init(f, "a220_b2");
00120 
00121   helper420a220beam1.Init(f, "a420a220");
00122   helper420a220beam2.Init(f, "a420a220_b2");
00123 
00124   f.Close();
00125 
00126   std::cout << "ProtonTaggerFilter: Ready" << std::endl;
00127 }
00128 
00129 
00132 void ProtonTaggerFilter::endJob() {;}
00133 
00134 
00137 bool ProtonTaggerFilter::filter(edm::Event & iEvent, const edm::EventSetup & es)
00138 {
00139   // ... get generated event
00140 
00141   edm::Handle<edm::HepMCProduct> evtSource;
00142   iEvent.getByLabel("generator",evtSource);
00143   const HepMC::GenEvent* genEvent = evtSource->GetEvent();
00144 
00145   //std::cout << "event contains " << genEvent->particles_size() << " particles " << std::endl;
00146   if (genEvent->particles_empty()) {
00147     std::cout << "empty source event" << std::endl;
00148     return false;
00149   }
00150 
00151   // ... get pileup event
00152 
00153   edm::Handle<edm::HepMCProduct> pileUpSource;
00154   const HepMC::GenEvent* pileUpEvent = 0;
00155   bool isPileUp = true;
00156   
00157   bool isProduct = iEvent.getByLabel("famosPileUp", "PileUpEvents", pileUpSource);
00158 
00159   if (isProduct) {
00160     pileUpEvent = pileUpSource->GetEvent();
00161     //std::cout << "got pileup" << std::endl;
00162   } else {
00163     isPileUp = false;
00164     //std::cout << "no pileup in the event" << std::endl;
00165   }
00166 
00167   //if (isPileUp) std::cout << "event contains " << pileUpEvent->particles_size() << " pileup particles " << std::endl;  
00168 
00169   // ... some constants
00170 
00171   const double mp = 0.938272029; // just a proton mass
00172   const double Ebeam = 7000.0;   // beam energy - would be better to read from parameters
00173   const float  pzCut = 2500;     // ignore particles with less than |Pz| < pzCut  
00174   const float  acceptThreshold = 0.5; // proton will be "accepted" if the value of acceptance is above this threshold
00175 
00176   // ... loop over particles, find the most energetic proton in either direction
00177 
00178   std::list<HepMC::GenParticle*> veryForwardParicles;
00179   
00180   for ( HepMC::GenEvent::particle_const_iterator piter = genEvent->particles_begin();
00181           piter != genEvent->particles_end();
00182           ++piter )
00183   {
00184     
00185     HepMC::GenParticle* p = *piter;
00186 
00187     float pz = p->momentum().pz();
00188     if (((pz > pzCut) || (pz < -pzCut)) && ((p->status() == 0) || (p->status() == 1))) {
00189       veryForwardParicles.push_back(p);
00190       //std::cout << "pdgid: " << p->pdg_id() << " status: " << p->status() << std::endl;
00191     }
00192   }
00193 
00194   //std::cout << "# generated forward particles  : " << veryForwardParicles.size() << std::endl;
00195 
00196   if ( isPileUp ) {
00197 
00198     //std::cout << "Adding pileup " << std::endl;
00199 
00200     for ( HepMC::GenEvent::particle_const_iterator piter = pileUpEvent->particles_begin();
00201             piter != pileUpEvent->particles_end();
00202             ++piter )
00203     {
00204       
00205       HepMC::GenParticle* p = *piter;
00206 
00207       float pz = p->momentum().pz();
00208       if (((pz > pzCut) || (pz < -pzCut)) && ((p->status() == 0) || (p->status() == 1))) {
00209         veryForwardParicles.push_back(p);
00210         //std::cout << "pdgid: " << p->pdg_id() << " status: " << p->status() << std::endl;
00211       }
00212     }
00213   }
00214 
00215   //std::cout << "# forward particles to be tried: " << veryForwardParicles.size() << std::endl;
00216 
00217 
00218   // ... return false if no forward protons found
00219 
00220   if (veryForwardParicles.empty()) return false;
00221 
00222 
00223   // ... set all acceptances to zero
00224 
00225   float acc420b1, acc220b1, acc420and220b1, acc420or220b1; // beam 1 (clockwise)
00226   float acc420b2, acc220b2, acc420and220b2, acc420or220b2; // beam 2 (anti-clockwise)
00227 
00228   acc420b1 = acc220b1 = acc420and220b1 = acc420or220b1 = 0;
00229   acc420b2 = acc220b2 = acc420and220b2 = acc420or220b2 = 0;
00230 
00231   int nP1at220m = 0;
00232   int nP1at420m = 0;
00233 
00234   int nP2at220m = 0;
00235   int nP2at420m = 0;
00236 
00237   // ... loop over (pre-selected) forward particles
00238   
00239   for (std::list<HepMC::GenParticle*>::const_iterator part = veryForwardParicles.begin();
00240        part != veryForwardParicles.end();
00241        part++)
00242   {
00243 
00244     HepMC::GenParticle* p = *part;
00245 
00246     float pz  = p->momentum().pz();
00247     float pt  = p->momentum().perp();
00248     float phi = p->momentum().phi();;
00249 
00250     if ((pz > Ebeam) || (pz < -Ebeam)) continue;
00251     
00252     // ... compute kinimatical variable
00253 
00254     float xi  = 1.0;    // fractional momentum loss
00255     if (pz>0) xi -= pz/Ebeam;
00256     else xi += pz/Ebeam;
00257 
00258     double t   = (-pt*pt - mp*mp*xi*xi) / (1-xi); // "t"
00259 
00260     //std::cout << " pdg_id: "  << p->pdg_id() << " eta: " << p->momentum().eta() << " e: " 
00261     //          <<  p->momentum().e() << std::endl;
00262     //std::cout << "pz: " << pz << " pt: " <<  pt << " xi: " << xi
00263     //          << " t: " << t << " phi: " << phi << std::endl;
00264 
00265     if (xi<0.0) xi=-10.0;
00266     if (xi>1.0) xi=10.0;
00267 
00268     //float rnd1 = RandGauss::shoot(0.,1.);
00269     //float rnd2 = RandGauss::shoot(0.,1.);
00270 
00271     // ... get acceptance from tables: beam 1 (if it is not ignored)
00272 
00273     if ((pz > 0) && (beam1mode != 0)) {    
00274 
00275       acc420b1       = helper420beam1.GetAcceptance(t, xi, phi);
00276       acc220b1       = helper220beam1.GetAcceptance(t, xi, phi);
00277       acc420and220b1 = helper420a220beam1.GetAcceptance(t, xi, phi);
00278 
00279       acc420or220b1  = acc420b1 + acc220b1 - acc420and220b1;
00280 
00281       //std::cout << "+acc420b1: " << acc420b1 << " acc220b1: " << acc220b1 << " acc420and220b1: " << acc420and220b1 << " acc420or220b1: " << acc420or220b1  << std::endl;
00282 
00283       bool res420and220 = (acc420and220b1 > acceptThreshold);
00284       bool res420       = (acc420b1       > acceptThreshold);
00285       bool res220       = (acc220b1       > acceptThreshold);
00286 
00287       if (res420and220) {nP1at220m++; nP1at420m++;}
00288       else if (res420) nP1at420m++;
00289       else if (res220) nP1at220m++;
00290   
00291       if ((p->pdg_id() != 2212) && (res220 || res420 || res420and220)) {
00292         std::cout << " !!! P got proton 1 at 420 m: pz = " << pz << std::endl;
00293         if (res220) std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m along beam 1, pz = " << pz << std::endl;
00294         if (res420) std::cout << "got a particle with pid" << p->pdg_id() << " at 420 m along beam 1, pz = " << pz << std::endl;
00295         if (res420and220) std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m & 420 m  along beam 1, pz = " << pz << std::endl;
00296       }
00297     }
00298     
00299     // ... the same for beam 2
00300 
00301     if ((pz < 0) && (beam2mode != 0)) {
00302 
00303       acc420b2       = helper420beam2.GetAcceptance(t, xi, phi);
00304       acc220b2       = helper220beam2.GetAcceptance(t, xi, phi);
00305       acc420and220b2 = helper420a220beam2.GetAcceptance(t, xi, phi);
00306     
00307       acc420or220b2  = acc420b2 + acc220b2 - acc420and220b2;
00308 
00309       //std::cout << "+acc420b2: " << acc420b2 << " acc220b2: " << acc220b2 << " acc420and220b2: " << acc420and220b2 << " acc420or220b2: " << acc420or220b2 << std::endl;
00310 
00311       bool res420and220 = (acc420and220b2 > acceptThreshold);
00312       bool res420       = (acc420b2       > acceptThreshold);
00313       bool res220       = (acc220b2       > acceptThreshold);
00314 
00315       if (res420and220) {nP2at220m++; nP2at420m++;}
00316       else if (res420) nP2at420m++;
00317       else if (res220) nP2at220m++;
00318   
00319       if ((p->pdg_id() != 2212) && (res220 || res420 || res420and220)) {
00320         std::cout << " !!! P got proton 1 at 420 m: pz = " << pz << std::endl;
00321         if (res220) std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m along beam 2, pz = " << pz << std::endl;
00322         if (res420) std::cout << "got a particle with pid" << p->pdg_id() << " at 420 m along beam 2, pz = " << pz << std::endl;
00323         if (res420and220) std::cout << "got a particle with pid" << p->pdg_id() << " at 220 m & 420 m along beam 2, pz = " << pz << std::endl;
00324       }
00325     }
00326   }
00327 
00328   // ... boolean result for each detector
00329   
00330   bool p1at220m = (nP1at220m>0) ? true : false;
00331   bool p1at420m = (nP1at420m>0) ? true : false;
00332   bool p2at220m = (nP2at220m>0) ? true : false;
00333   bool p2at420m = (nP2at420m>0) ? true : false;
00334 
00335   if ((nP1at220m > 1) && (beam1mode!=1)) std::cout << " !!! " << nP1at220m << " proton(s) from beam 1 at 220 m" << std::endl;
00336   if ((nP1at420m > 1) && (beam1mode!=2)) std::cout << " !!! " << nP1at420m << " proton(s) from beam 1 at 420 m" << std::endl;
00337   if ((nP2at220m > 1) && (beam2mode!=1)) std::cout << " !!! " << nP2at220m << " proton(s) from beam 2 at 220 m" << std::endl;
00338   if ((nP2at420m > 1) && (beam2mode!=2)) std::cout << " !!! " << nP2at420m << " proton(s) from beam 2 at 420 m" << std::endl;
00339 
00340 
00341   // ... make a decision based on requested filter configuration
00342     
00343   bool p1accepted = false;
00344   bool p2accepted = false;
00345   
00346   if ((beam1mode==1) &&  p1at420m ) p1accepted = true;
00347   if ((beam1mode==2) &&  p1at220m ) p1accepted = true;
00348   if ((beam1mode==3) &&  p1at220m && p1at420m) p1accepted = true;
00349   if ((beam1mode==4) && (p1at220m || p1at420m)) p1accepted = true;
00350   
00351   if ((beam2mode==1) &&  p2at420m ) p2accepted = true;
00352   if ((beam2mode==2) &&  p2at220m ) p2accepted = true;
00353   if ((beam2mode==3) &&  p2at220m && p2at420m ) p2accepted = true;
00354   if ((beam2mode==4) && (p2at220m || p2at420m)) p2accepted = true;
00355 
00356   //if (p1accepted) std::cout << "proton 1 accepted" << std::endl;
00357   //if (p2accepted) std::cout << "proton 2 accepted" << std::endl;
00358 
00359   switch (beamCombiningMode) {
00360 
00361     case 1:  // ... either of two protons
00362      if (p1accepted || p2accepted) return true; else return false;
00363 
00364     case 2:  // ... both protons
00365      if (p1accepted && p2accepted) return true; else return false;
00366 
00367     case 3:  // ... 220+220 or 420+420
00368      if ((p1at220m && p2at220m) || (p1at420m && p2at420m)) return true; else return false;
00369 
00370     case 4:  // ... 220+420 or 420+220
00371      if ((p1at220m && p2at420m) || (p1at420m && p2at220m)) return true; else return false;
00372 
00373   }
00374   
00375   return false;
00376 }
00377 
00378 // ... define CMSSW module
00379 
00380 DEFINE_FWK_MODULE(ProtonTaggerFilter);
/data/doxygen/doxygen-1.7.3/gen/CMSSW_4_2_8/src/FastSimulation/ForwardDetectors/plugins/ProtonTaggerFilter.cc
