Utilities/interface/FWLiteTools.h

#ifndef ForwardAnalysis_FWLiteTools_h
#define ForwardAnalysis_FWLiteTools_h

#if !defined(__CINT__) && !defined(__MAKECINT__)
#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"

#include <vector>
#include <algorithm>
#endif

namespace forwardAnalysis {

enum calo_region_t {Barrel,Endcap,Transition,Forward};

bool sortByEta( const math::XYZTLorentzVector& a, const math::XYZTLorentzVector& b){ 
   return a.eta() < b.eta();
}

void genRapidityGap(reco::GenParticleCollection const& genParticles, math::XYZTLorentzVector& genGapLow,
                                                                     math::XYZTLorentzVector& genGapHigh,
                                                                     double etaEdgeLow = -999.0,
                                                                     double etaEdgeHigh = 999.0){
   //double etaEdgeLow = -999.0;
   //double etaEdgeHigh = 999.0;
   // Copy and sort gen particles in eta
   std::vector<math::XYZTLorentzVector> genParticlesSort(0);
   reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
   reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();  
   for(; genpart != genpart_end; ++genpart){
      if( genpart->status() != 1 ) continue;

      if((genpart->eta() >= etaEdgeLow) && (genpart->eta() <= etaEdgeHigh))
         genParticlesSort.push_back( genpart->p4() );
   }
   std::stable_sort(genParticlesSort.begin(), genParticlesSort.end(), sortByEta);

   // Cases: 0, 1 or > 1 particles in selected range
   math::XYZTLorentzVector def_vec(0.,0.,0.,0.);
   if( genParticlesSort.size() == 0 ){
      genGapLow = def_vec; genGapHigh = def_vec;
      return;
   } else if( genParticlesSort.size() == 1 ){
      genGapLow = def_vec;
      genGapHigh = genParticlesSort[0];
      return;
   } else{
      //FIXME; There must be a STL algorithm for this
      double deltaEtaMax = 0.;
      std::vector<math::XYZTLorentzVector>::const_iterator genPartDeltaEtaMax = genParticlesSort.end();
      std::vector<math::XYZTLorentzVector>::const_iterator genpart = genParticlesSort.begin();
      std::vector<math::XYZTLorentzVector>::const_iterator genpart_end = genParticlesSort.end();
      for(; genpart != genpart_end; ++genpart){
         std::vector<math::XYZTLorentzVector>::const_iterator next = genpart + 1;
         double deltaEta = ( next != genpart_end ) ? ( next->eta() - genpart->eta() ) : 0.;
         if( deltaEta > (deltaEtaMax - 0.0001) ){
            deltaEtaMax = deltaEta;
            genPartDeltaEtaMax = genpart;
         } 
      }
      if( genPartDeltaEtaMax != genpart_end ){
         std::vector<math::XYZTLorentzVector>::const_iterator next = genPartDeltaEtaMax + 1;
         if( next != genpart_end ){
            genGapLow = (*genPartDeltaEtaMax);
            genGapHigh = (*next);
         } else{
            genGapLow = def_vec;
            genGapHigh = (*genPartDeltaEtaMax);
         }
      } else{
         genGapLow = def_vec; genGapHigh = def_vec;
         return;
      }
   } 
   
}

void setGenInfo(reco::GenParticleCollection const& genParticles, double Ebeam,
                                                                 math::XYZTLorentzVector& genAllParticles,
                                                                 math::XYZTLorentzVector& genAllParticlesInRange,
                                                                 math::XYZTLorentzVector& genAllParticlesHEPlus,
                                                                 math::XYZTLorentzVector& genAllParticlesHEMinus,
                                                                 math::XYZTLorentzVector& genAllParticlesHFPlus,
                                                                 math::XYZTLorentzVector& genAllParticlesHFMinus,
                                                                 math::XYZTLorentzVector& genEtaMax,
                                                                 math::XYZTLorentzVector& genEtaMin,
                                                                 math::XYZTLorentzVector& genProtonPlus,
                                                                 math::XYZTLorentzVector& genProtonMinus){

   math::XYZTLorentzVector allGenParticles(0.,0.,0.,0.);
   math::XYZTLorentzVector allGenParticlesInRange(0.,0.,0.,0.);
   math::XYZTLorentzVector allGenParticlesHEPlus(0.,0.,0.,0.);
   math::XYZTLorentzVector allGenParticlesHEMinus(0.,0.,0.,0.);
   math::XYZTLorentzVector allGenParticlesHFPlus(0.,0.,0.,0.);
   math::XYZTLorentzVector allGenParticlesHFMinus(0.,0.,0.,0.);
   
   reco::GenParticleCollection::const_iterator proton1 = genParticles.end();
   reco::GenParticleCollection::const_iterator proton2 = genParticles.end();
   for(reco::GenParticleCollection::const_iterator genpart = genParticles.begin(); genpart != genParticles.end();
 ++genpart){
      if( genpart->status() != 1 ) continue;
      if( genpart->pdgId() != 2212 ) continue;
 
      if( ( genpart->pz() > 0.50*Ebeam ) && ( ( proton1 == genParticles.end() ) ||
                                              ( genpart->pz() > proton1->pz() ) ) ) proton1 = genpart;
      if( ( genpart->pz() < -0.50*Ebeam ) && ( ( proton2 == genParticles.end() ) ||
                                               ( genpart->pz() < proton2->pz() ) ) ) proton2 = genpart;
   }

   reco::GenParticleCollection::const_iterator etaMaxParticle = genParticles.end();
   reco::GenParticleCollection::const_iterator etaMinParticle = genParticles.end(); 
   for(reco::GenParticleCollection::const_iterator genpart = genParticles.begin(); genpart != genParticles.end(); ++genpart){
      if(genpart->status() != 1) continue;

      allGenParticles += genpart->p4();
      if(fabs(genpart->eta()) < 5.0) allGenParticlesInRange += genpart->p4();
      if( (genpart->eta() >= 1.3) && (genpart->eta() < 3.0) ) allGenParticlesHEPlus += genpart->p4();
      if( (genpart->eta() > -3.0) && (genpart->eta() <= -1.3) ) allGenParticlesHEMinus += genpart->p4();
      if( (genpart->eta() >= 3.0) && (genpart->eta() < 5.0) ) allGenParticlesHFPlus += genpart->p4();
      if( (genpart->eta() > -5.0) && (genpart->eta() <= -3.0) ) allGenParticlesHFMinus += genpart->p4(); 

      if( (genpart != proton1) && (genpart != proton2) ){
         if( ( etaMaxParticle == genParticles.end() ) ||
             ( genpart->eta() > etaMaxParticle->eta() ) ) etaMaxParticle = genpart;
         if( ( etaMinParticle == genParticles.end() ) ||
             ( genpart->eta() < etaMinParticle->eta() ) ) etaMinParticle = genpart;
      }
   }

   // Commit
   if( proton1 != genParticles.end() ) allGenParticles -= proton1->p4();
   if( proton2 != genParticles.end() ) allGenParticles -= proton2->p4();

   genAllParticles = allGenParticles;
   genAllParticlesInRange = allGenParticlesInRange;
   genAllParticlesHEPlus = allGenParticlesHEPlus;
   genAllParticlesHEMinus = allGenParticlesHEMinus;
   genAllParticlesHFPlus = allGenParticlesHFPlus;
   genAllParticlesHFMinus = allGenParticlesHFMinus;
 
   if( proton1 != genParticles.end() ) genProtonPlus = proton1->p4();
   if( proton2 != genParticles.end() ) genProtonMinus = proton2->p4();

   if( etaMaxParticle != genParticles.end() ) genEtaMax = etaMaxParticle->p4();
   if( etaMinParticle != genParticles.end() ) genEtaMin = etaMinParticle->p4();
}

int pflowId(std::string const& name){
   // FIXME: The labels definition could go somewhere else
   std::vector<std::string> labels_X, labels_h, labels_e, labels_mu, labels_gamma, labels_h0, labels_h_HF, labels_egamma_HF;
   labels_X.push_back("X");
   labels_X.push_back("undefined");
   labels_h.push_back("h");
   labels_h.push_back("chargedHadron");
   labels_h.push_back("hadronCharged");
   labels_e.push_back("e");
   labels_e.push_back("electron");
   labels_mu.push_back("mu");
   labels_mu.push_back("muon");
   labels_gamma.push_back("gamma");
   labels_gamma.push_back("photon");
   labels_h0.push_back("h0");
   labels_h0.push_back("neutralHadron");
   labels_h0.push_back("hadronNeutral");
   labels_h_HF.push_back("h_HF");
   labels_h_HF.push_back("hadronHF");
   labels_egamma_HF.push_back("egamma_HF");
   labels_egamma_HF.push_back("emHF");
   // Find corresponding particle type   
   if( std::find(labels_X.begin(), labels_X.end(), name) != labels_X.end() )
      return reco::PFCandidate::X;
   else if( std::find(labels_h.begin(), labels_h.end(), name) != labels_h.end() )
      return reco::PFCandidate::h;
   else if( std::find(labels_e.begin(), labels_e.end(), name) != labels_e.end() )
      return reco::PFCandidate::e;
   else if( std::find(labels_mu.begin(), labels_mu.end(), name) != labels_mu.end() )
      return reco::PFCandidate::mu;
   else if( std::find(labels_gamma.begin(), labels_gamma.end(), name) != labels_gamma.end() ) 
      return reco::PFCandidate::gamma;
   else if( std::find(labels_h0.begin(), labels_h0.end(), name) != labels_h0.end() ) 
      return reco::PFCandidate::h0;
   else if( std::find(labels_h_HF.begin(), labels_h_HF.end(), name) != labels_h_HF.end() ) 
      return reco::PFCandidate::h_HF;
   else if( std::find(labels_egamma_HF.begin(), labels_egamma_HF.end(), name) != labels_egamma_HF.end() ) 
      return reco::PFCandidate::egamma_HF;
   else return -1;
}

bool pflowThreshold(reco::PFCandidate const& part, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){

   //FIXME
   // HF eta rings 29, 40, 41
   if( ( (fabs(part.eta()) >= 2.866) && (fabs(part.eta()) < 2.976) ) || 
         (fabs(part.eta()) >= 4.730) ) return false;

   bool accept = true;

   double eta = part.eta();
   int region = -1;
   if( (fabs(eta) >= 0.) && (fabs(eta) < 1.4) ) region = Barrel;
   else if( (fabs(eta) >= 1.4) && (fabs(eta) < 2.6) ) region = Endcap;
   else if( (fabs(eta) >= 2.6) && (fabs(eta) < 3.2) ) region = Transition;
   else if( (fabs(eta) >= 3.2) ) region = Forward;
   std::map<int,std::pair<double,double> > const& thresholds = thresholdMap.find(region)->second;
   
   double ptThreshold = -1.0;
   double eThreshold = -1.0;
   int partType = part.particleId();
   std::map<int,std::pair<double,double> >::const_iterator it_threshold = thresholds.find(partType);
   if(it_threshold != thresholds.end()) {
      ptThreshold = it_threshold->second.first;
      eThreshold = it_threshold->second.second;
   }

   if(part.pt() < ptThreshold) accept = false;
   if(part.energy() < eThreshold) accept = false;

   return accept;
}

/*bool caloTowerNoiseAccept(CaloTower const& tower, double emFracThreshold = 1.){
   bool accept = true;

   double emEnergy = tower.emEnergy();
   double hadEnergy = tower.hadEnergy();
   double emFrac = fabs(emEnergy/(emEnergy+hadEnergy));
   if(emFrac > emFracThreshold) accept = false;

   return accept;
}*/

double MassDissGen(reco::GenParticleCollection const& genParticles, double rangeEtaMin = -999.,
                                                                    double rangeEtaMax = 999.){
                                                                    
   math::XYZTLorentzVector allGenParticles(0.,0.,0.,0.);
   reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
   reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();
   for(; genpart != genpart_end; ++genpart){
      if( genpart->status() != 1 ) continue;

      if( ( genpart->eta() >= (rangeEtaMin - 0.0001) ) && 
          ( genpart->eta() <= (rangeEtaMax + 0.0001) ) ) allGenParticles += genpart->p4();
   }
   return allGenParticles.M();
}

// FIXME: Generalize for any collection with changeable threshold scheme
double MassColl(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
   math::XYZTLorentzVector allCands(0.,0.,0.,0.);
   reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
   reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
   for(; part != pfCands_end; ++part){
      if(pflowThreshold(*part,thresholdMap)) allCands += part->p4();
   }

   return allCands.M();
}

std::pair<double,double> xi(reco::PFCandidateCollection const& pflowCollection, double Ebeam, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){

   double xi_towers_plus = 0.;
   double xi_towers_minus = 0.;
   reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
   reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
   for(; part != pfCands_end; ++part){
     if(!pflowThreshold(*part,thresholdMap)) continue;

     xi_towers_plus += part->et()*TMath::Exp(part->eta());
     xi_towers_minus += part->et()*TMath::Exp(-part->eta());
   }

   xi_towers_plus /= 2*Ebeam;
   xi_towers_minus /= 2*Ebeam;
   
   return std::make_pair(xi_towers_plus,xi_towers_minus);
}

std::pair<double,double> EPlusPz(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
   double e_plus_pz = 0.;
   double e_minus_pz = 0.;
   reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
   reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
   for(; part != pfCands_end; ++part){
      if(!pflowThreshold(*part,thresholdMap)) continue;

      e_plus_pz += part->energy() + part->pz(); 
      e_minus_pz += part->energy() - part->pz();
   }

   return std::make_pair(e_plus_pz,e_minus_pz);
}

std::pair<double,double> etaMax(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
   std::vector<double> etaCands;
   reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
   reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
   for(; part != pfCands_end; ++part){                           
      if(!pflowThreshold(*part,thresholdMap)) continue;            
      etaCands.push_back( part->eta() );
   }                                                             
   double eta_max = etaCands.size() ? *( std::max_element(etaCands.begin(), etaCands.end()) ) : -999.;
   double eta_min = etaCands.size() ? *( std::min_element(etaCands.begin(), etaCands.end()) ) : -999.;

   return std::make_pair(eta_max,eta_min);
}

double castorEnergy(CastorRecHitCollection const& castorRecHitCollection, bool isRealData = true){
 
   double sumETotCastor = 0.,
          sumETotCastorNMod5 = 0.,
          sumETotCastorNMod4 = 0.,
          sumETotCastorNMod3 = 0.,
          sumETotCastorNMod2 = 0.;
  
   // Loop over rec hits
   CastorRecHitCollection::const_iterator castorRecHit = castorRecHitCollection.begin();
   CastorRecHitCollection::const_iterator castorRecHits_end = castorRecHitCollection.end();
   for(; castorRecHit != castorRecHits_end; ++castorRecHit) {
      const CastorRecHit& recHit = (*castorRecHit);

      int sectorId  = recHit.id().sector();
      int moduleId  = recHit.id().module();
      double energy = recHit.energy();
      double time   = recHit.time();

      if( !isRealData ) energy *= 62.5;

      if( moduleId > 5 ) continue;

      if( moduleId == 1 && sectorId == 5 ) continue;
      if( moduleId == 1 && sectorId == 6) continue;

      sumETotCastor += energy;

      if( moduleId <= 5 ) sumETotCastorNMod5 += energy;
      if( moduleId <= 4 ) sumETotCastorNMod4 += energy;
      if( moduleId <= 3 ) sumETotCastorNMod3 += energy;
      if( moduleId <= 2 ) sumETotCastorNMod2 += energy;
   }

   return sumETotCastor;
}

} // namespace
#endif
Revision:	1.1
Committed:	Wed Aug 17 15:13:01 2011 UTC (13 years, 8 months ago) by antoniov
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V01-01-01, V01-01-00, antoniov-forwardAnalysis-09Jul2012-v1, antoniov-forwardAnalysis-29Jun2012-v1, V01-00-00, antoniov-utilities-11Jun2012-v1, antoniov-forwardAnalysis-Oct072011-v1, sfonseca_10_04_2011, antoniov-forwardAnalysis-Sep182011-v1, antoniov-forwardAnalysis-Sep102011-v1, sfonseca_08_26_2011, forwardAnalysis-Aug232011-v1, forwardAnalysis-Aug172011-v1, HEAD
Log Message:	add FWLiteTools
#	User	Rev	Content
1	antoniov	1.1	#ifndef ForwardAnalysis_FWLiteTools_h
2			#define ForwardAnalysis_FWLiteTools_h
3
4			#if !defined(__CINT__) && !defined(__MAKECINT__)
5			#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
6			#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
7			#include "DataFormats/CaloTowers/interface/CaloTower.h"
8			#include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
9			#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
10
11			#include <vector>
12			#include <algorithm>
13			#endif
14
15			namespace forwardAnalysis {
16
17			enum calo_region_t {Barrel,Endcap,Transition,Forward};
18
19			bool sortByEta( const math::XYZTLorentzVector& a, const math::XYZTLorentzVector& b){
20			return a.eta() < b.eta();
21			}
22
23			void genRapidityGap(reco::GenParticleCollection const& genParticles, math::XYZTLorentzVector& genGapLow,
24			math::XYZTLorentzVector& genGapHigh,
25			double etaEdgeLow = -999.0,
26			double etaEdgeHigh = 999.0){
27			//double etaEdgeLow = -999.0;
28			//double etaEdgeHigh = 999.0;
29			// Copy and sort gen particles in eta
30			std::vector<math::XYZTLorentzVector> genParticlesSort(0);
31			reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
32			reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();
33			for(; genpart != genpart_end; ++genpart){
34			if( genpart->status() != 1 ) continue;
35
36			if((genpart->eta() >= etaEdgeLow) && (genpart->eta() <= etaEdgeHigh))
37			genParticlesSort.push_back( genpart->p4() );
38			}
39			std::stable_sort(genParticlesSort.begin(), genParticlesSort.end(), sortByEta);
40
41			// Cases: 0, 1 or > 1 particles in selected range
42			math::XYZTLorentzVector def_vec(0.,0.,0.,0.);
43			if( genParticlesSort.size() == 0 ){
44			genGapLow = def_vec; genGapHigh = def_vec;
45			return;
46			} else if( genParticlesSort.size() == 1 ){
47			genGapLow = def_vec;
48			genGapHigh = genParticlesSort[0];
49			return;
50			} else{
51			//FIXME; There must be a STL algorithm for this
52			double deltaEtaMax = 0.;
53			std::vector<math::XYZTLorentzVector>::const_iterator genPartDeltaEtaMax = genParticlesSort.end();
54			std::vector<math::XYZTLorentzVector>::const_iterator genpart = genParticlesSort.begin();
55			std::vector<math::XYZTLorentzVector>::const_iterator genpart_end = genParticlesSort.end();
56			for(; genpart != genpart_end; ++genpart){
57			std::vector<math::XYZTLorentzVector>::const_iterator next = genpart + 1;
58			double deltaEta = ( next != genpart_end ) ? ( next->eta() - genpart->eta() ) : 0.;
59			if( deltaEta > (deltaEtaMax - 0.0001) ){
60			deltaEtaMax = deltaEta;
61			genPartDeltaEtaMax = genpart;
62			}
63			}
64			if( genPartDeltaEtaMax != genpart_end ){
65			std::vector<math::XYZTLorentzVector>::const_iterator next = genPartDeltaEtaMax + 1;
66			if( next != genpart_end ){
67			genGapLow = (*genPartDeltaEtaMax);
68			genGapHigh = (*next);
69			} else{
70			genGapLow = def_vec;
71			genGapHigh = (*genPartDeltaEtaMax);
72			}
73			} else{
74			genGapLow = def_vec; genGapHigh = def_vec;
75			return;
76			}
77			}
78
79			}
80
81			void setGenInfo(reco::GenParticleCollection const& genParticles, double Ebeam,
82			math::XYZTLorentzVector& genAllParticles,
83			math::XYZTLorentzVector& genAllParticlesInRange,
84			math::XYZTLorentzVector& genAllParticlesHEPlus,
85			math::XYZTLorentzVector& genAllParticlesHEMinus,
86			math::XYZTLorentzVector& genAllParticlesHFPlus,
87			math::XYZTLorentzVector& genAllParticlesHFMinus,
88			math::XYZTLorentzVector& genEtaMax,
89			math::XYZTLorentzVector& genEtaMin,
90			math::XYZTLorentzVector& genProtonPlus,
91			math::XYZTLorentzVector& genProtonMinus){
92
93			math::XYZTLorentzVector allGenParticles(0.,0.,0.,0.);
94			math::XYZTLorentzVector allGenParticlesInRange(0.,0.,0.,0.);
95			math::XYZTLorentzVector allGenParticlesHEPlus(0.,0.,0.,0.);
96			math::XYZTLorentzVector allGenParticlesHEMinus(0.,0.,0.,0.);
97			math::XYZTLorentzVector allGenParticlesHFPlus(0.,0.,0.,0.);
98			math::XYZTLorentzVector allGenParticlesHFMinus(0.,0.,0.,0.);
99
100			reco::GenParticleCollection::const_iterator proton1 = genParticles.end();
101			reco::GenParticleCollection::const_iterator proton2 = genParticles.end();
102			for(reco::GenParticleCollection::const_iterator genpart = genParticles.begin(); genpart != genParticles.end();
103			++genpart){
104			if( genpart->status() != 1 ) continue;
105			if( genpart->pdgId() != 2212 ) continue;
106
107			if( ( genpart->pz() > 0.50*Ebeam ) && ( ( proton1 == genParticles.end() ) \|\|
108			( genpart->pz() > proton1->pz() ) ) ) proton1 = genpart;
109			if( ( genpart->pz() < -0.50*Ebeam ) && ( ( proton2 == genParticles.end() ) \|\|
110			( genpart->pz() < proton2->pz() ) ) ) proton2 = genpart;
111			}
112
113			reco::GenParticleCollection::const_iterator etaMaxParticle = genParticles.end();
114			reco::GenParticleCollection::const_iterator etaMinParticle = genParticles.end();
115			for(reco::GenParticleCollection::const_iterator genpart = genParticles.begin(); genpart != genParticles.end(); ++genpart){
116			if(genpart->status() != 1) continue;
117
118			allGenParticles += genpart->p4();
119			if(fabs(genpart->eta()) < 5.0) allGenParticlesInRange += genpart->p4();
120			if( (genpart->eta() >= 1.3) && (genpart->eta() < 3.0) ) allGenParticlesHEPlus += genpart->p4();
121			if( (genpart->eta() > -3.0) && (genpart->eta() <= -1.3) ) allGenParticlesHEMinus += genpart->p4();
122			if( (genpart->eta() >= 3.0) && (genpart->eta() < 5.0) ) allGenParticlesHFPlus += genpart->p4();
123			if( (genpart->eta() > -5.0) && (genpart->eta() <= -3.0) ) allGenParticlesHFMinus += genpart->p4();
124
125			if( (genpart != proton1) && (genpart != proton2) ){
126			if( ( etaMaxParticle == genParticles.end() ) \|\|
127			( genpart->eta() > etaMaxParticle->eta() ) ) etaMaxParticle = genpart;
128			if( ( etaMinParticle == genParticles.end() ) \|\|
129			( genpart->eta() < etaMinParticle->eta() ) ) etaMinParticle = genpart;
130			}
131			}
132
133			// Commit
134			if( proton1 != genParticles.end() ) allGenParticles -= proton1->p4();
135			if( proton2 != genParticles.end() ) allGenParticles -= proton2->p4();
136
137			genAllParticles = allGenParticles;
138			genAllParticlesInRange = allGenParticlesInRange;
139			genAllParticlesHEPlus = allGenParticlesHEPlus;
140			genAllParticlesHEMinus = allGenParticlesHEMinus;
141			genAllParticlesHFPlus = allGenParticlesHFPlus;
142			genAllParticlesHFMinus = allGenParticlesHFMinus;
143
144			if( proton1 != genParticles.end() ) genProtonPlus = proton1->p4();
145			if( proton2 != genParticles.end() ) genProtonMinus = proton2->p4();
146
147			if( etaMaxParticle != genParticles.end() ) genEtaMax = etaMaxParticle->p4();
148			if( etaMinParticle != genParticles.end() ) genEtaMin = etaMinParticle->p4();
149			}
150
151			int pflowId(std::string const& name){
152			// FIXME: The labels definition could go somewhere else
153			std::vector<std::string> labels_X, labels_h, labels_e, labels_mu, labels_gamma, labels_h0, labels_h_HF, labels_egamma_HF;
154			labels_X.push_back("X");
155			labels_X.push_back("undefined");
156			labels_h.push_back("h");
157			labels_h.push_back("chargedHadron");
158			labels_h.push_back("hadronCharged");
159			labels_e.push_back("e");
160			labels_e.push_back("electron");
161			labels_mu.push_back("mu");
162			labels_mu.push_back("muon");
163			labels_gamma.push_back("gamma");
164			labels_gamma.push_back("photon");
165			labels_h0.push_back("h0");
166			labels_h0.push_back("neutralHadron");
167			labels_h0.push_back("hadronNeutral");
168			labels_h_HF.push_back("h_HF");
169			labels_h_HF.push_back("hadronHF");
170			labels_egamma_HF.push_back("egamma_HF");
171			labels_egamma_HF.push_back("emHF");
172			// Find corresponding particle type
173			if( std::find(labels_X.begin(), labels_X.end(), name) != labels_X.end() )
174			return reco::PFCandidate::X;
175			else if( std::find(labels_h.begin(), labels_h.end(), name) != labels_h.end() )
176			return reco::PFCandidate::h;
177			else if( std::find(labels_e.begin(), labels_e.end(), name) != labels_e.end() )
178			return reco::PFCandidate::e;
179			else if( std::find(labels_mu.begin(), labels_mu.end(), name) != labels_mu.end() )
180			return reco::PFCandidate::mu;
181			else if( std::find(labels_gamma.begin(), labels_gamma.end(), name) != labels_gamma.end() )
182			return reco::PFCandidate::gamma;
183			else if( std::find(labels_h0.begin(), labels_h0.end(), name) != labels_h0.end() )
184			return reco::PFCandidate::h0;
185			else if( std::find(labels_h_HF.begin(), labels_h_HF.end(), name) != labels_h_HF.end() )
186			return reco::PFCandidate::h_HF;
187			else if( std::find(labels_egamma_HF.begin(), labels_egamma_HF.end(), name) != labels_egamma_HF.end() )
188			return reco::PFCandidate::egamma_HF;
189			else return -1;
190			}
191
192			bool pflowThreshold(reco::PFCandidate const& part, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
193
194			//FIXME
195			// HF eta rings 29, 40, 41
196			if( ( (fabs(part.eta()) >= 2.866) && (fabs(part.eta()) < 2.976) ) \|\|
197			(fabs(part.eta()) >= 4.730) ) return false;
198
199			bool accept = true;
200
201			double eta = part.eta();
202			int region = -1;
203			if( (fabs(eta) >= 0.) && (fabs(eta) < 1.4) ) region = Barrel;
204			else if( (fabs(eta) >= 1.4) && (fabs(eta) < 2.6) ) region = Endcap;
205			else if( (fabs(eta) >= 2.6) && (fabs(eta) < 3.2) ) region = Transition;
206			else if( (fabs(eta) >= 3.2) ) region = Forward;
207			std::map<int,std::pair<double,double> > const& thresholds = thresholdMap.find(region)->second;
208
209			double ptThreshold = -1.0;
210			double eThreshold = -1.0;
211			int partType = part.particleId();
212			std::map<int,std::pair<double,double> >::const_iterator it_threshold = thresholds.find(partType);
213			if(it_threshold != thresholds.end()) {
214			ptThreshold = it_threshold->second.first;
215			eThreshold = it_threshold->second.second;
216			}
217
218			if(part.pt() < ptThreshold) accept = false;
219			if(part.energy() < eThreshold) accept = false;
220
221			return accept;
222			}
223
224			/*bool caloTowerNoiseAccept(CaloTower const& tower, double emFracThreshold = 1.){
225			bool accept = true;
226
227			double emEnergy = tower.emEnergy();
228			double hadEnergy = tower.hadEnergy();
229			double emFrac = fabs(emEnergy/(emEnergy+hadEnergy));
230			if(emFrac > emFracThreshold) accept = false;
231
232			return accept;
233			}*/
234
235			double MassDissGen(reco::GenParticleCollection const& genParticles, double rangeEtaMin = -999.,
236			double rangeEtaMax = 999.){
237
238			math::XYZTLorentzVector allGenParticles(0.,0.,0.,0.);
239			reco::GenParticleCollection::const_iterator genpart = genParticles.begin();
240			reco::GenParticleCollection::const_iterator genpart_end = genParticles.end();
241			for(; genpart != genpart_end; ++genpart){
242			if( genpart->status() != 1 ) continue;
243
244			if( ( genpart->eta() >= (rangeEtaMin - 0.0001) ) &&
245			( genpart->eta() <= (rangeEtaMax + 0.0001) ) ) allGenParticles += genpart->p4();
246			}
247			return allGenParticles.M();
248			}
249
250			// FIXME: Generalize for any collection with changeable threshold scheme
251			double MassColl(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
252			math::XYZTLorentzVector allCands(0.,0.,0.,0.);
253			reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
254			reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
255			for(; part != pfCands_end; ++part){
256			if(pflowThreshold(*part,thresholdMap)) allCands += part->p4();
257			}
258
259			return allCands.M();
260			}
261
262			std::pair<double,double> xi(reco::PFCandidateCollection const& pflowCollection, double Ebeam, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
263
264			double xi_towers_plus = 0.;
265			double xi_towers_minus = 0.;
266			reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
267			reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
268			for(; part != pfCands_end; ++part){
269			if(!pflowThreshold(*part,thresholdMap)) continue;
270
271			xi_towers_plus += part->et()*TMath::Exp(part->eta());
272			xi_towers_minus += part->et()*TMath::Exp(-part->eta());
273			}
274
275			xi_towers_plus /= 2*Ebeam;
276			xi_towers_minus /= 2*Ebeam;
277
278			return std::make_pair(xi_towers_plus,xi_towers_minus);
279			}
280
281			std::pair<double,double> EPlusPz(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
282			double e_plus_pz = 0.;
283			double e_minus_pz = 0.;
284			reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
285			reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
286			for(; part != pfCands_end; ++part){
287			if(!pflowThreshold(*part,thresholdMap)) continue;
288
289			e_plus_pz += part->energy() + part->pz();
290			e_minus_pz += part->energy() - part->pz();
291			}
292
293			return std::make_pair(e_plus_pz,e_minus_pz);
294			}
295
296			std::pair<double,double> etaMax(reco::PFCandidateCollection const& pflowCollection, std::map<int, std::map<int,std::pair<double,double> > > const& thresholdMap){
297			std::vector<double> etaCands;
298			reco::PFCandidateCollection::const_iterator part = pflowCollection.begin();
299			reco::PFCandidateCollection::const_iterator pfCands_end = pflowCollection.end();
300			for(; part != pfCands_end; ++part){
301			if(!pflowThreshold(*part,thresholdMap)) continue;
302			etaCands.push_back( part->eta() );
303			}
304			double eta_max = etaCands.size() ? *( std::max_element(etaCands.begin(), etaCands.end()) ) : -999.;
305			double eta_min = etaCands.size() ? *( std::min_element(etaCands.begin(), etaCands.end()) ) : -999.;
306
307			return std::make_pair(eta_max,eta_min);
308			}
309
310			double castorEnergy(CastorRecHitCollection const& castorRecHitCollection, bool isRealData = true){
311
312			double sumETotCastor = 0.,
313			sumETotCastorNMod5 = 0.,
314			sumETotCastorNMod4 = 0.,
315			sumETotCastorNMod3 = 0.,
316			sumETotCastorNMod2 = 0.;
317
318			// Loop over rec hits
319			CastorRecHitCollection::const_iterator castorRecHit = castorRecHitCollection.begin();
320			CastorRecHitCollection::const_iterator castorRecHits_end = castorRecHitCollection.end();
321			for(; castorRecHit != castorRecHits_end; ++castorRecHit) {
322			const CastorRecHit& recHit = (*castorRecHit);
323
324			int sectorId = recHit.id().sector();
325			int moduleId = recHit.id().module();
326			double energy = recHit.energy();
327			double time = recHit.time();
328
329			if( !isRealData ) energy *= 62.5;
330
331			if( moduleId > 5 ) continue;
332
333			if( moduleId == 1 && sectorId == 5 ) continue;
334			if( moduleId == 1 && sectorId == 6) continue;
335
336			sumETotCastor += energy;
337
338			if( moduleId <= 5 ) sumETotCastorNMod5 += energy;
339			if( moduleId <= 4 ) sumETotCastorNMod4 += energy;
340			if( moduleId <= 3 ) sumETotCastorNMod3 += energy;
341			if( moduleId <= 2 ) sumETotCastorNMod2 += energy;
342			}
343
344			return sumETotCastor;
345			}
346
347			} // namespace
348			#endif