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
#	Content
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