SFrameAnalysis/src/Cleaner.cxx

#include "../include/Cleaner.h"


Cleaner::Cleaner( BaseCycleContainer* input_){

  bcc = input_;

}

void Cleaner::JetEnergyResolutionShifter(int syst_shift){

  double met = 0;
  if(bcc->met) met=bcc->met->pt();

  for(unsigned int i=0; i<bcc->jets->size(); ++i){
    Jet jet = bcc->jets->at(i);
    float gen_pt = jet.genjet_pt();
    //ignore unmatched jets (which have zero vector) or jets with very low pt:
    if(gen_pt < 15.0) continue;

    met += jet.pt()*jet.JEC_factor_raw();
    
    float recopt = jet.pt();
    double factor = -10;
    double abseta = fabs(jet.eta());

    //numbers taken from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
    if(syst_shift==0){
      if(abseta < 0.5)
        factor = 0.052;
      else if(abseta >= 0.5 && abseta <1.1)
        factor = 0.057;
      else if(abseta >= 1.1 && abseta <1.7)
        factor = 0.096;
      else if(abseta >= 1.7 && abseta <2.3)
        factor = 0.134;
      else if(abseta >= 2.3)
        factor = 0.288;
    }
    else if(syst_shift>0){
      if(abseta < 0.5)
        factor = 0.11;
      else if(abseta >= 0.5 && abseta <1.1)
        factor = 0.12;
      else if(abseta >= 1.1 && abseta <1.7)
        factor = 0.16;
      else if(abseta >= 1.7 && abseta <2.3)
        factor = 0.23;
      else if(abseta >= 2.3)
        factor = 0.49;
    }
    else{
      if(abseta < 0.5)
        factor = -0.01;
      else if(abseta >= 0.5 && abseta <1.1)
        factor = 0.00;
      else if(abseta >= 1.1 && abseta <1.7)
        factor = 0.04;
      else if(abseta >= 1.7 && abseta <2.3)
        factor = 0.03;
      else if(abseta >= 2.3)
        factor = 0.09;
    }

    float deltapt = (recopt - gen_pt) * factor;
    float ptscale = std::max(0.0f, (recopt + deltapt) / recopt);
    jet.set_pt(jet.pt() * ptscale);
    bcc->jets->at(i).set_pt( jet.pt());

    //propagate JER shifts to MET
    met -= jet.pt()*jet.JEC_factor_raw();
  }
  
  //store changed MET, flip phi if new MET is negative
  if(bcc->met){
    if(met>=0){
      bcc->met->set_pt( met);
    }
    else{
      bcc->met->set_pt( -1*met);
      if(bcc->met->phi()<=0)
        bcc->met->set_phi (bcc->met->phi()+PI);
      else
        bcc->met->set_phi ( bcc->met->phi()-PI);
    }
  }

  sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
}


//tight ele ID from https://twiki.cern.ch/twiki/bin/view/CMS/EgammaCutBasedIdentification
bool Cleaner::eleID(Electron ele){

  bool pass=false;
  
  if(fabs(ele.supercluster_eta())<1.4442){
    if(ele.dEtaIn()<0.004 && ele.dPhiIn()<0.03 && ele.sigmaIEtaIEta()<0.01 && ele.HoverE()<0.12 && fabs(1./ele.EcalEnergy()-1./ele.v4().P())<0.05) pass=true; 
  }
  else if( fabs(ele.supercluster_eta())>1.5660){
    if(ele.dEtaIn()<0.005 && ele.dPhiIn()<0.02 && ele.sigmaIEtaIEta()<0.03 && ele.HoverE()<0.10 && fabs(1./ele.EcalEnergy()-1./ele.v4().P())<0.05) pass=true; 
  }

  return pass;

}

//pf ID has already been applied when using goodPatJets
bool Cleaner::pfID(Jet jet){

  if(jet.numberOfDaughters()>1 
     && jet.neutralHadronEnergyFraction()<0.99
     && jet.neutralEmEnergyFraction()<0.99){
    
    if(fabs(jet.eta())>=2.4)
      return true;
    
    if(jet.chargedEmEnergyFraction()<0.99
       && jet.chargedHadronEnergyFraction()>0
       && jet.chargedMultiplicity()>0)
      return true;

  }
  return false;
     
}


void Cleaner::ElectronCleaner(double ptmin, double etamax){

  std::vector<Electron> good_eles;
  for(unsigned int i=0; i<bcc->electrons->size(); ++i){
    Electron ele = bcc->electrons->at(i);
    if(ele.pt()>ptmin){
      if(fabs(ele.eta())<etamax){
        if(fabs(ele.supercluster_eta())<1.4442 || fabs(ele.supercluster_eta())>1.5660){
          if(bcc->pvs->size()>0){
            if(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y())<0.02){
              if(ele.passconversionveto()){
                //if(ele.mvaTrigV0>0.0){
                if(eleID(ele)){
                  if(ele.relIso()<0.1){
                    good_eles.push_back(ele);
                  }
                }
                //}
              }
            }
          }
        }
      }
    }
  }

  bcc->electrons->clear();

  for(unsigned int i=0; i<good_eles.size(); ++i){
    bcc->electrons->push_back(good_eles[i]);
  }
  sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
}


void Cleaner::MuonCleaner(double ptmin, double etamax){

  std::vector<Muon> good_mus;
  for(unsigned int i=0; i<bcc->muons->size(); ++i){
    Muon mu = bcc->muons->at(i);
    if(mu.pt()>ptmin){
      if(fabs(mu.eta())<etamax){
        if(mu.isGlobalMuon()){
          if(mu.globalTrack_chi2()/mu.globalTrack_ndof()<10){
            if(mu.innerTrack_trackerLayersWithMeasurement()>8){
              if(mu.dB()<0.02){
                if(mu.relIso()<0.125){
                  if(fabs(mu.vertex_z()-bcc->pvs->at(0).z())<1){
                    if(mu.innerTrack_numberOfValidPixelHits()>0){
                      if(mu.numberOfMatchedStations()>1){
                        good_mus.push_back(mu);
                      }
                    }
                  }
                }
              }
            }
          }
        }
      }
    }
  }
  bcc->muons->clear();

  for(unsigned int i=0; i<good_mus.size(); ++i){
    bcc->muons->push_back(good_mus[i]);
  }
  sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
}

void Cleaner::TauCleaner(double ptmin, double etamax){

  std::vector<Tau> good_taus;
  for(unsigned int i=0; i<bcc->taus->size(); ++i){
    Tau tau = bcc->taus->at(i);
    if(tau.pt()>ptmin){
      if(fabs(tau.eta())<etamax){
        if(bcc->taus->at(i).decayModeFinding()){
          if(bcc->taus->at(i).byMediumCombinedIsolationDeltaBetaCorr()){
            if(bcc->taus->at(i).againstElectronTight()){
              if(bcc->taus->at(i).againstMuonTight()){
                good_taus.push_back(tau);
              }
            }
          }
        }
      }
    }
  }
  
  bcc->taus->clear();

  for(unsigned int i=0; i<good_taus.size(); ++i){
    bcc->taus->push_back(good_taus[i]);
  }
  sort(bcc->taus->begin(), bcc->taus->end(), HigherPt());
}


void Cleaner::JetCleaner(double ptmin, double etamax, bool doPFID){

  std::vector<Jet> good_jets;
  for(unsigned int i=0; i<bcc->jets->size(); ++i){
    Jet jet = bcc->jets->at(i);
    if(jet.pt()>ptmin){
      if(fabs(jet.eta())<etamax){
        if(!doPFID || pfID(jet)){
          good_jets.push_back(jet);
        }
      }
    }
  }

  bcc->jets->clear();

  for(unsigned int i=0; i<good_jets.size(); ++i){
    bcc->jets->push_back(good_jets[i]);
  }
  sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
}

void Cleaner::TopJetCleaner(double ptmin, double etamax, bool doPFID){

  std::vector<TopJet> good_topjets;
  for(unsigned int i=0; i<bcc->topjets->size(); ++i){
    TopJet topjet = bcc->topjets->at(i);
    if(topjet.pt()>ptmin){
      if(fabs(topjet.eta())<etamax){
        if(!doPFID || pfID(topjet)){
          good_topjets.push_back(topjet);
        }
      }
    }
  }
  bcc->topjets->clear();

  for(unsigned int i=0; i<good_topjets.size(); ++i){
    bcc->topjets->push_back(good_topjets[i]);
  }
  sort(bcc->topjets->begin(), bcc->topjets->end(), HigherPt());
}
Revision:	1.6
Committed:	Fri May 25 09:32:31 2012 UTC (12 years, 11 months ago) by peiffer
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+14 -10 lines
Log Message:	move files to SFrameTools
#	User	Rev	Content
1	peiffer	1.1	#include "../include/Cleaner.h"
2
3
4			Cleaner::Cleaner( BaseCycleContainer* input_){
5
6			bcc = input_;
7
8			}
9
10	peiffer	1.2	void Cleaner::JetEnergyResolutionShifter(int syst_shift){
11
12	peiffer	1.6	double met = 0;
13			if(bcc->met) met=bcc->met->pt();
14	peiffer	1.2
15			for(unsigned int i=0; i<bcc->jets->size(); ++i){
16			Jet jet = bcc->jets->at(i);
17	peiffer	1.5	float gen_pt = jet.genjet_pt();
18	peiffer	1.2	//ignore unmatched jets (which have zero vector) or jets with very low pt:
19			if(gen_pt < 15.0) continue;
20
21	peiffer	1.5	met += jet.pt()*jet.JEC_factor_raw();
22	peiffer	1.2
23	peiffer	1.5	float recopt = jet.pt();
24	peiffer	1.2	double factor = -10;
25	peiffer	1.5	double abseta = fabs(jet.eta());
26	peiffer	1.2
27			//numbers taken from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
28			if(syst_shift==0){
29			if(abseta < 0.5)
30			factor = 0.052;
31			else if(abseta >= 0.5 && abseta <1.1)
32			factor = 0.057;
33			else if(abseta >= 1.1 && abseta <1.7)
34			factor = 0.096;
35			else if(abseta >= 1.7 && abseta <2.3)
36			factor = 0.134;
37			else if(abseta >= 2.3)
38			factor = 0.288;
39			}
40			else if(syst_shift>0){
41			if(abseta < 0.5)
42			factor = 0.11;
43			else if(abseta >= 0.5 && abseta <1.1)
44			factor = 0.12;
45			else if(abseta >= 1.1 && abseta <1.7)
46			factor = 0.16;
47			else if(abseta >= 1.7 && abseta <2.3)
48			factor = 0.23;
49			else if(abseta >= 2.3)
50			factor = 0.49;
51			}
52			else{
53			if(abseta < 0.5)
54			factor = -0.01;
55			else if(abseta >= 0.5 && abseta <1.1)
56			factor = 0.00;
57			else if(abseta >= 1.1 && abseta <1.7)
58			factor = 0.04;
59			else if(abseta >= 1.7 && abseta <2.3)
60			factor = 0.03;
61			else if(abseta >= 2.3)
62			factor = 0.09;
63			}
64
65			float deltapt = (recopt - gen_pt) * factor;
66			float ptscale = std::max(0.0f, (recopt + deltapt) / recopt);
67	peiffer	1.5	jet.set_pt(jet.pt() * ptscale);
68			bcc->jets->at(i).set_pt( jet.pt());
69	peiffer	1.2
70			//propagate JER shifts to MET
71	peiffer	1.5	met -= jet.pt()*jet.JEC_factor_raw();
72	peiffer	1.2	}
73
74			//store changed MET, flip phi if new MET is negative
75	peiffer	1.6	if(bcc->met){
76			if(met>=0){
77			bcc->met->set_pt( met);
78			}
79			else{
80			bcc->met->set_pt( -1*met);
81			if(bcc->met->phi()<=0)
82			bcc->met->set_phi (bcc->met->phi()+PI);
83			else
84			bcc->met->set_phi ( bcc->met->phi()-PI);
85			}
86	peiffer	1.2	}
87	peiffer	1.6
88	peiffer	1.5	sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
89	peiffer	1.2	}
90
91	peiffer	1.1
92			//tight ele ID from https://twiki.cern.ch/twiki/bin/view/CMS/EgammaCutBasedIdentification
93			bool Cleaner::eleID(Electron ele){
94
95			bool pass=false;
96
97	peiffer	1.5	if(fabs(ele.supercluster_eta())<1.4442){
98			if(ele.dEtaIn()<0.004 && ele.dPhiIn()<0.03 && ele.sigmaIEtaIEta()<0.01 && ele.HoverE()<0.12 && fabs(1./ele.EcalEnergy()-1./ele.v4().P())<0.05) pass=true;
99	peiffer	1.1	}
100	peiffer	1.5	else if( fabs(ele.supercluster_eta())>1.5660){
101			if(ele.dEtaIn()<0.005 && ele.dPhiIn()<0.02 && ele.sigmaIEtaIEta()<0.03 && ele.HoverE()<0.10 && fabs(1./ele.EcalEnergy()-1./ele.v4().P())<0.05) pass=true;
102	peiffer	1.1	}
103
104			return pass;
105
106			}
107
108	peiffer	1.2	//pf ID has already been applied when using goodPatJets
109			bool Cleaner::pfID(Jet jet){
110
111	peiffer	1.5	if(jet.numberOfDaughters()>1
112			&& jet.neutralHadronEnergyFraction()<0.99
113			&& jet.neutralEmEnergyFraction()<0.99){
114	peiffer	1.2
115	peiffer	1.5	if(fabs(jet.eta())>=2.4)
116	peiffer	1.2	return true;
117
118	peiffer	1.5	if(jet.chargedEmEnergyFraction()<0.99
119			&& jet.chargedHadronEnergyFraction()>0
120			&& jet.chargedMultiplicity()>0)
121	peiffer	1.2	return true;
122
123			}
124			return false;
125
126			}
127
128	peiffer	1.1
129	peiffer	1.3	void Cleaner::ElectronCleaner(double ptmin, double etamax){
130	peiffer	1.1
131	peiffer	1.3	std::vector<Electron> good_eles;
132	peiffer	1.1	for(unsigned int i=0; i<bcc->electrons->size(); ++i){
133			Electron ele = bcc->electrons->at(i);
134	peiffer	1.5	if(ele.pt()>ptmin){
135			if(fabs(ele.eta())<etamax){
136			if(fabs(ele.supercluster_eta())<1.4442 \|\| fabs(ele.supercluster_eta())>1.5660){
137	peiffer	1.1	if(bcc->pvs->size()>0){
138	peiffer	1.5	if(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y())<0.02){
139			if(ele.passconversionveto()){
140	peiffer	1.1	//if(ele.mvaTrigV0>0.0){
141			if(eleID(ele)){
142			if(ele.relIso()<0.1){
143	peiffer	1.3	good_eles.push_back(ele);
144	peiffer	1.1	}
145			}
146			//}
147			}
148			}
149			}
150			}
151			}
152			}
153			}
154
155	peiffer	1.3	bcc->electrons->clear();
156
157			for(unsigned int i=0; i<good_eles.size(); ++i){
158			bcc->electrons->push_back(good_eles[i]);
159			}
160	peiffer	1.5	sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
161	peiffer	1.1	}
162
163
164	peiffer	1.3	void Cleaner::MuonCleaner(double ptmin, double etamax){
165
166			std::vector<Muon> good_mus;
167	peiffer	1.1	for(unsigned int i=0; i<bcc->muons->size(); ++i){
168			Muon mu = bcc->muons->at(i);
169	peiffer	1.5	if(mu.pt()>ptmin){
170			if(fabs(mu.eta())<etamax){
171			if(mu.isGlobalMuon()){
172			if(mu.globalTrack_chi2()/mu.globalTrack_ndof()<10){
173			if(mu.innerTrack_trackerLayersWithMeasurement()>8){
174			if(mu.dB()<0.02){
175	peiffer	1.4	if(mu.relIso()<0.125){
176	peiffer	1.5	if(fabs(mu.vertex_z()-bcc->pvs->at(0).z())<1){
177			if(mu.innerTrack_numberOfValidPixelHits()>0){
178			if(mu.numberOfMatchedStations()>1){
179	peiffer	1.3	good_mus.push_back(mu);
180	peiffer	1.1	}
181			}
182			}
183			}
184			}
185			}
186			}
187			}
188			}
189			}
190			}
191	peiffer	1.3	bcc->muons->clear();
192
193			for(unsigned int i=0; i<good_mus.size(); ++i){
194			bcc->muons->push_back(good_mus[i]);
195			}
196	peiffer	1.5	sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
197	peiffer	1.1	}
198	peiffer	1.2
199	peiffer	1.4	void Cleaner::TauCleaner(double ptmin, double etamax){
200
201			std::vector<Tau> good_taus;
202			for(unsigned int i=0; i<bcc->taus->size(); ++i){
203			Tau tau = bcc->taus->at(i);
204	peiffer	1.5	if(tau.pt()>ptmin){
205			if(fabs(tau.eta())<etamax){
206			if(bcc->taus->at(i).decayModeFinding()){
207			if(bcc->taus->at(i).byMediumCombinedIsolationDeltaBetaCorr()){
208			if(bcc->taus->at(i).againstElectronTight()){
209			if(bcc->taus->at(i).againstMuonTight()){
210	peiffer	1.4	good_taus.push_back(tau);
211			}
212			}
213			}
214			}
215			}
216			}
217			}
218
219			bcc->taus->clear();
220
221			for(unsigned int i=0; i<good_taus.size(); ++i){
222			bcc->taus->push_back(good_taus[i]);
223			}
224	peiffer	1.5	sort(bcc->taus->begin(), bcc->taus->end(), HigherPt());
225	peiffer	1.4	}
226
227	peiffer	1.2
228	peiffer	1.3	void Cleaner::JetCleaner(double ptmin, double etamax, bool doPFID){
229
230			std::vector<Jet> good_jets;
231	peiffer	1.2	for(unsigned int i=0; i<bcc->jets->size(); ++i){
232			Jet jet = bcc->jets->at(i);
233	peiffer	1.5	if(jet.pt()>ptmin){
234			if(fabs(jet.eta())<etamax){
235	peiffer	1.2	if(!doPFID \|\| pfID(jet)){
236	peiffer	1.3	good_jets.push_back(jet);
237	peiffer	1.2	}
238			}
239			}
240			}
241
242	peiffer	1.3	bcc->jets->clear();
243
244			for(unsigned int i=0; i<good_jets.size(); ++i){
245			bcc->jets->push_back(good_jets[i]);
246			}
247	peiffer	1.5	sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
248	peiffer	1.2	}
249
250	peiffer	1.3	void Cleaner::TopJetCleaner(double ptmin, double etamax, bool doPFID){
251	peiffer	1.2
252	peiffer	1.3	std::vector<TopJet> good_topjets;
253	peiffer	1.2	for(unsigned int i=0; i<bcc->topjets->size(); ++i){
254			TopJet topjet = bcc->topjets->at(i);
255	peiffer	1.5	if(topjet.pt()>ptmin){
256			if(fabs(topjet.eta())<etamax){
257	peiffer	1.2	if(!doPFID \|\| pfID(topjet)){
258	peiffer	1.3	good_topjets.push_back(topjet);
259	peiffer	1.2	}
260			}
261			}
262			}
263	peiffer	1.3	bcc->topjets->clear();
264	peiffer	1.2
265	peiffer	1.3	for(unsigned int i=0; i<good_topjets.size(); ++i){
266			bcc->topjets->push_back(good_topjets[i]);
267			}
268	peiffer	1.5	sort(bcc->topjets->begin(), bcc->topjets->end(), HigherPt());
269	peiffer	1.2	}