HiGenCommon/src/EcalGenEvtSelectorFrag.cc

0001 #include <iostream>
0002 #include "GeneratorInterface/HiGenCommon/interface/EcalGenEvtSelectorFrag.h"
0003 #include "FWCore/Utilities/interface/EDMException.h"
0004 #include "HepMC/GenVertex.h"
0005 using namespace std;
0006
0007 EcalGenEvtSelectorFrag::EcalGenEvtSelectorFrag(const edm::ParameterSet& pset) : BaseHiGenEvtSelector(pset) {
0008   partonId_ = pset.getParameter<vector<int> >("partons");
0009   partonStatus_ = pset.getParameter<vector<int> >("partonStatus");
0010   partonPt_ = pset.getParameter<vector<double> >("partonPt");
0011
0012   particleId_ = pset.getParameter<vector<int> >("particles");
0013   particleStatus_ = pset.getParameter<vector<int> >("particleStatus");
0014   particlePt_ = pset.getParameter<vector<double> >("particlePt");
0015
0016   etaMax_ = pset.getParameter<double>("etaMax");
0017
0018   int id = partonId_.size();
0019   int st = partonStatus_.size();
0020   int pt = partonPt_.size();
0021
0022   if (partonId_.size() != partonStatus_.size() || partonId_.size() != partonPt_.size()) {
0023     throw edm::Exception(edm::errors::LogicError) << id << st << pt << endl;
0024   }
0025
0026   id = particleId_.size();
0027   st = particleStatus_.size();
0028   pt = particlePt_.size();
0029
0030   if (particleId_.size() != particleStatus_.size() || particleId_.size() != particlePt_.size()) {
0031     throw edm::Exception(edm::errors::LogicError) << id << st << pt << endl;
0032   }
0033 }
0034
0035 bool EcalGenEvtSelectorFrag::filter(HepMC::GenEvent* evt) {
0036   HepMC::GenEvent::particle_const_iterator begin = evt->particles_begin();
0037   HepMC::GenEvent::particle_const_iterator end = evt->particles_end();
0038
0039   bool foundParticle = false;
0040   bool foundParton = false;
0041
0042   HepMC::GenEvent::particle_const_iterator it = begin;
0043   while ((!foundParton || !foundParticle) && it != end) {
0044     bool isFrag = false;
0045     bool isThisPhoton = false;
0046
0047     foundParton = true;
0048     /*for(unsigned i = 0; i < partonId_.size(); ++i){
0049        if(selectParticle(*it, partonStatus_[i], partonId_[i], partonPt_[i], etaMax_)) foundParton = true;
0050        }*/
0051
0052     for (unsigned i = 0; i < particleId_.size(); ++i) {
0053       if (selectParticle(*it, particleStatus_[i], particleId_[i], particlePt_[i], etaMax_))
0054         isThisPhoton = true;
0055     }
0056
0057     // Now you have to requre the partcile is "prompt", meaning its mom is parton
0058
0059     if (!((*it)->production_vertex())) {
0060       isFrag = false;
0061     } else {
0062       const HepMC::GenVertex* productionVertex = (*it)->production_vertex();
0063
0064       size_t numberOfMothers = productionVertex->particles_in_size();
0065       if (numberOfMothers <= 0) {
0066         isFrag = false;
0067         //   cout << "number of mothers = " << numberOfMothers << endl;
0068       } else {
0069         //   cout << "number of mothers = " << numberOfMothers << endl;
0070         HepMC::GenVertex::particles_in_const_iterator motherIt = productionVertex->particles_in_const_begin();
0071         for (; motherIt != productionVertex->particles_in_const_end(); motherIt++) {
0072           if (fabs((*motherIt)->pdg_id()) <= 22) {
0073             isFrag = true;
0074           }
0075         }
0076       }
0077     }
0078
0079     if ((isFrag == true) && (isThisPhoton == true)) {
0080       //cout << "Found fragmentation photon!!" << endl ;
0081       foundParticle = true;
0082     }
0083
0084     ++it;
0085   }
0086
0087   foundParton = true;  // We don't care of the parton
0088   return (foundParton && foundParticle);
0089 }