SFrameAnalysis/src/Cleaner.cxx

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


Cleaner::Cleaner( BaseCycleContainer* input){

  bcc = input;

}

Cleaner::Cleaner(){

  ObjectHandler* objs = ObjectHandler::Instance();
  bcc = objs->GetBaseCycleContainer();

}

void Cleaner::resetEventCalc(){

  EventCalc* evc = EventCalc::Instance();
  evc->Reset();

}

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());
  resetEventCalc();
}


void Cleaner::JetLeptonSubtractor(FactorizedJetCorrector *corrector){

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

  for(unsigned int i=0; i<bcc->jets->size(); ++i){

    LorentzVector jet_v4_raw = bcc->jets->at(i).v4()*bcc->jets->at(i).JEC_factor_raw();

    met-=bcc->jets->at(i).pt();
  
    //subtract lepton momenta from raw jet momentum
    if(bcc->electrons){
      for(unsigned int j=0; j<bcc->electrons->size(); ++j){
        if(bcc->jets->at(i).deltaR(bcc->electrons->at(j))<0.5){
          jet_v4_raw -= bcc->electrons->at(j).v4();
        }
      }
    }
    if(bcc->muons){
      for(unsigned int j=0; j<bcc->muons->size(); ++j){
        jet_v4_raw -= bcc->muons->at(j).v4();
      }
    }

    //apply jet energy corrections to modified raw momentum
    corrector->setJetPt(jet_v4_raw.Pt());
    corrector->setJetEta(jet_v4_raw.Eta());
    corrector->setJetE(jet_v4_raw.E());  
    corrector->setJetA(bcc->jets->at(i).jetArea());
    corrector->setRho(bcc->rho);

    float correctionfactor = corrector->getCorrection();
    LorentzVector jet_v4_corrected = jet_v4_raw *correctionfactor;

    bcc->jets->at(i).set_v4(jet_v4_corrected);
    bcc->jets->at(i).set_JEC_factor_raw(1./correctionfactor);

    met+=bcc->jets->at(i).pt();
    
  }

  //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());
  resetEventCalc();
}

//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_noIso(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)){
                    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::ElectronCleaner(double ptmin, double etamax, double relisomax){

  ElectronCleaner_noIso(ptmin, etamax);
  std::vector<Electron> good_eles;
  for(unsigned int i=0; i<bcc->electrons->size(); ++i){
    Electron ele = bcc->electrons->at(i);
    if(ele.relIso()<relisomax){
      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());
  resetEventCalc();
}


void Cleaner::MuonCleaner_noIso(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(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());
  resetEventCalc();
}

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

  MuonCleaner_noIso(ptmin, etamax);

  std::vector<Muon> good_mus;
  for(unsigned int i=0; i<bcc->muons->size(); ++i){
    Muon mu = bcc->muons->at(i);
    if(mu.relIso()<relisomax){
      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());
  resetEventCalc();
}


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());
  resetEventCalc();
}

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());
  resetEventCalc();
}


void Cleaner::PrimaryVertexCleaner(int ndofmax, double zmax, double rhomax){

  std::vector<PrimaryVertex> good_pvs;
  for(unsigned int i=0; i<bcc->pvs->size(); ++i){
    PrimaryVertex pv = bcc->pvs->at(i);
    if(pv.ndof()>=ndofmax && fabs(pv.z())<zmax && pv.rho()<=rhomax){
      good_pvs.push_back(pv);
    } 
  }
  bcc->pvs->clear();

  for(unsigned int i=0; i<good_pvs.size(); ++i){
    bcc->pvs->push_back(good_pvs[i]);
  }
  resetEventCalc();
}
Revision:	1.9
Committed:	Wed Jun 6 15:50:55 2012 UTC (12 years, 10 months ago) by peiffer
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+23 -4 lines
Log Message:	use EventCalc and ObjectHandler
#	User	Rev	Content
1	peiffer	1.1	#include "../include/Cleaner.h"
2
3
4	peiffer	1.7	Cleaner::Cleaner( BaseCycleContainer* input){
5	peiffer	1.1
6	peiffer	1.7	bcc = input;
7	peiffer	1.1
8			}
9
10	peiffer	1.9	Cleaner::Cleaner(){
11
12			ObjectHandler* objs = ObjectHandler::Instance();
13			bcc = objs->GetBaseCycleContainer();
14
15			}
16
17			void Cleaner::resetEventCalc(){
18
19			EventCalc* evc = EventCalc::Instance();
20			evc->Reset();
21
22			}
23
24	peiffer	1.2	void Cleaner::JetEnergyResolutionShifter(int syst_shift){
25
26	peiffer	1.6	double met = 0;
27			if(bcc->met) met=bcc->met->pt();
28	peiffer	1.2
29			for(unsigned int i=0; i<bcc->jets->size(); ++i){
30			Jet jet = bcc->jets->at(i);
31	peiffer	1.5	float gen_pt = jet.genjet_pt();
32	peiffer	1.2	//ignore unmatched jets (which have zero vector) or jets with very low pt:
33			if(gen_pt < 15.0) continue;
34
35	peiffer	1.5	met += jet.pt()*jet.JEC_factor_raw();
36	peiffer	1.2
37	peiffer	1.5	float recopt = jet.pt();
38	peiffer	1.2	double factor = -10;
39	peiffer	1.5	double abseta = fabs(jet.eta());
40	peiffer	1.2
41			//numbers taken from https://twiki.cern.ch/twiki/bin/view/CMS/JetResolution
42			if(syst_shift==0){
43			if(abseta < 0.5)
44			factor = 0.052;
45			else if(abseta >= 0.5 && abseta <1.1)
46			factor = 0.057;
47			else if(abseta >= 1.1 && abseta <1.7)
48			factor = 0.096;
49			else if(abseta >= 1.7 && abseta <2.3)
50			factor = 0.134;
51			else if(abseta >= 2.3)
52			factor = 0.288;
53			}
54			else if(syst_shift>0){
55			if(abseta < 0.5)
56			factor = 0.11;
57			else if(abseta >= 0.5 && abseta <1.1)
58			factor = 0.12;
59			else if(abseta >= 1.1 && abseta <1.7)
60			factor = 0.16;
61			else if(abseta >= 1.7 && abseta <2.3)
62			factor = 0.23;
63			else if(abseta >= 2.3)
64			factor = 0.49;
65			}
66			else{
67			if(abseta < 0.5)
68			factor = -0.01;
69			else if(abseta >= 0.5 && abseta <1.1)
70			factor = 0.00;
71			else if(abseta >= 1.1 && abseta <1.7)
72			factor = 0.04;
73			else if(abseta >= 1.7 && abseta <2.3)
74			factor = 0.03;
75			else if(abseta >= 2.3)
76			factor = 0.09;
77			}
78
79			float deltapt = (recopt - gen_pt) * factor;
80			float ptscale = std::max(0.0f, (recopt + deltapt) / recopt);
81	peiffer	1.5	jet.set_pt(jet.pt() * ptscale);
82			bcc->jets->at(i).set_pt( jet.pt());
83	peiffer	1.2
84			//propagate JER shifts to MET
85	peiffer	1.5	met -= jet.pt()*jet.JEC_factor_raw();
86	peiffer	1.2	}
87
88			//store changed MET, flip phi if new MET is negative
89	peiffer	1.6	if(bcc->met){
90			if(met>=0){
91			bcc->met->set_pt( met);
92			}
93			else{
94			bcc->met->set_pt( -1*met);
95			if(bcc->met->phi()<=0)
96			bcc->met->set_phi (bcc->met->phi()+PI);
97			else
98	peiffer	1.7	bcc->met->set_phi (bcc->met->phi()-PI);
99	peiffer	1.6	}
100	peiffer	1.2	}
101	peiffer	1.6
102	peiffer	1.5	sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
103	peiffer	1.9	resetEventCalc();
104	peiffer	1.2	}
105
106	peiffer	1.1
107	peiffer	1.7	void Cleaner::JetLeptonSubtractor(FactorizedJetCorrector *corrector){
108
109			double met = 0;
110			if(bcc->met) met=bcc->met->pt();
111
112			for(unsigned int i=0; i<bcc->jets->size(); ++i){
113
114			LorentzVector jet_v4_raw = bcc->jets->at(i).v4()*bcc->jets->at(i).JEC_factor_raw();
115
116			met-=bcc->jets->at(i).pt();
117	peiffer	1.9
118	peiffer	1.7	//subtract lepton momenta from raw jet momentum
119			if(bcc->electrons){
120			for(unsigned int j=0; j<bcc->electrons->size(); ++j){
121			if(bcc->jets->at(i).deltaR(bcc->electrons->at(j))<0.5){
122			jet_v4_raw -= bcc->electrons->at(j).v4();
123			}
124			}
125			}
126			if(bcc->muons){
127			for(unsigned int j=0; j<bcc->muons->size(); ++j){
128			jet_v4_raw -= bcc->muons->at(j).v4();
129			}
130			}
131
132			//apply jet energy corrections to modified raw momentum
133			corrector->setJetPt(jet_v4_raw.Pt());
134			corrector->setJetEta(jet_v4_raw.Eta());
135			corrector->setJetE(jet_v4_raw.E());
136			corrector->setJetA(bcc->jets->at(i).jetArea());
137			corrector->setRho(bcc->rho);
138
139			float correctionfactor = corrector->getCorrection();
140			LorentzVector jet_v4_corrected = jet_v4_raw *correctionfactor;
141
142			bcc->jets->at(i).set_v4(jet_v4_corrected);
143			bcc->jets->at(i).set_JEC_factor_raw(1./correctionfactor);
144
145			met+=bcc->jets->at(i).pt();
146
147			}
148
149			//store changed MET, flip phi if new MET is negative
150			if(bcc->met){
151			if(met>=0){
152			bcc->met->set_pt( met);
153			}
154			else{
155			bcc->met->set_pt( -1*met);
156			if(bcc->met->phi()<=0)
157			bcc->met->set_phi (bcc->met->phi()+PI);
158			else
159			bcc->met->set_phi (bcc->met->phi()-PI);
160			}
161			}
162
163	peiffer	1.8	sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
164	peiffer	1.9	resetEventCalc();
165	peiffer	1.7	}
166
167	peiffer	1.1	//tight ele ID from https://twiki.cern.ch/twiki/bin/view/CMS/EgammaCutBasedIdentification
168			bool Cleaner::eleID(Electron ele){
169
170			bool pass=false;
171
172	peiffer	1.5	if(fabs(ele.supercluster_eta())<1.4442){
173			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;
174	peiffer	1.1	}
175	peiffer	1.5	else if( fabs(ele.supercluster_eta())>1.5660){
176			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;
177	peiffer	1.1	}
178
179			return pass;
180
181			}
182
183	peiffer	1.2	//pf ID has already been applied when using goodPatJets
184			bool Cleaner::pfID(Jet jet){
185
186	peiffer	1.5	if(jet.numberOfDaughters()>1
187			&& jet.neutralHadronEnergyFraction()<0.99
188			&& jet.neutralEmEnergyFraction()<0.99){
189	peiffer	1.2
190	peiffer	1.5	if(fabs(jet.eta())>=2.4)
191	peiffer	1.2	return true;
192
193	peiffer	1.5	if(jet.chargedEmEnergyFraction()<0.99
194			&& jet.chargedHadronEnergyFraction()>0
195			&& jet.chargedMultiplicity()>0)
196	peiffer	1.2	return true;
197
198			}
199			return false;
200
201			}
202
203	peiffer	1.1
204	peiffer	1.7	void Cleaner::ElectronCleaner_noIso(double ptmin, double etamax){
205	peiffer	1.1
206	peiffer	1.3	std::vector<Electron> good_eles;
207	peiffer	1.1	for(unsigned int i=0; i<bcc->electrons->size(); ++i){
208			Electron ele = bcc->electrons->at(i);
209	peiffer	1.5	if(ele.pt()>ptmin){
210			if(fabs(ele.eta())<etamax){
211			if(fabs(ele.supercluster_eta())<1.4442 \|\| fabs(ele.supercluster_eta())>1.5660){
212	peiffer	1.1	if(bcc->pvs->size()>0){
213	peiffer	1.5	if(ele.gsfTrack_dxy_vertex(bcc->pvs->at(0).x(), bcc->pvs->at(0).y())<0.02){
214			if(ele.passconversionveto()){
215	peiffer	1.1	//if(ele.mvaTrigV0>0.0){
216			if(eleID(ele)){
217	peiffer	1.3	good_eles.push_back(ele);
218	peiffer	1.1	}
219			//}
220			}
221			}
222			}
223			}
224			}
225			}
226			}
227
228	peiffer	1.3	bcc->electrons->clear();
229
230			for(unsigned int i=0; i<good_eles.size(); ++i){
231			bcc->electrons->push_back(good_eles[i]);
232			}
233	peiffer	1.5	sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
234	peiffer	1.1	}
235
236	peiffer	1.7	void Cleaner::ElectronCleaner(double ptmin, double etamax, double relisomax){
237	peiffer	1.1
238	peiffer	1.7	ElectronCleaner_noIso(ptmin, etamax);
239			std::vector<Electron> good_eles;
240			for(unsigned int i=0; i<bcc->electrons->size(); ++i){
241			Electron ele = bcc->electrons->at(i);
242			if(ele.relIso()<relisomax){
243			good_eles.push_back(ele);
244			}
245			}
246			bcc->electrons->clear();
247
248			for(unsigned int i=0; i<good_eles.size(); ++i){
249			bcc->electrons->push_back(good_eles[i]);
250			}
251			sort(bcc->electrons->begin(), bcc->electrons->end(), HigherPt());
252	peiffer	1.9	resetEventCalc();
253	peiffer	1.7	}
254
255
256			void Cleaner::MuonCleaner_noIso(double ptmin, double etamax){
257	peiffer	1.3
258			std::vector<Muon> good_mus;
259	peiffer	1.1	for(unsigned int i=0; i<bcc->muons->size(); ++i){
260			Muon mu = bcc->muons->at(i);
261	peiffer	1.5	if(mu.pt()>ptmin){
262			if(fabs(mu.eta())<etamax){
263			if(mu.isGlobalMuon()){
264			if(mu.globalTrack_chi2()/mu.globalTrack_ndof()<10){
265			if(mu.innerTrack_trackerLayersWithMeasurement()>8){
266			if(mu.dB()<0.02){
267	peiffer	1.7	if(fabs(mu.vertex_z()-bcc->pvs->at(0).z())<1){
268			if(mu.innerTrack_numberOfValidPixelHits()>0){
269			if(mu.numberOfMatchedStations()>1){
270			good_mus.push_back(mu);
271	peiffer	1.1	}
272			}
273			}
274			}
275			}
276			}
277			}
278			}
279			}
280			}
281	peiffer	1.3	bcc->muons->clear();
282
283			for(unsigned int i=0; i<good_mus.size(); ++i){
284			bcc->muons->push_back(good_mus[i]);
285			}
286	peiffer	1.5	sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
287	peiffer	1.9	resetEventCalc();
288	peiffer	1.1	}
289	peiffer	1.2
290	peiffer	1.7	void Cleaner::MuonCleaner(double ptmin, double etamax, double relisomax){
291
292			MuonCleaner_noIso(ptmin, etamax);
293
294			std::vector<Muon> good_mus;
295			for(unsigned int i=0; i<bcc->muons->size(); ++i){
296			Muon mu = bcc->muons->at(i);
297			if(mu.relIso()<relisomax){
298			good_mus.push_back(mu);
299			}
300			}
301			bcc->muons->clear();
302
303			for(unsigned int i=0; i<good_mus.size(); ++i){
304			bcc->muons->push_back(good_mus[i]);
305			}
306			sort(bcc->muons->begin(), bcc->muons->end(), HigherPt());
307			}
308
309	peiffer	1.4	void Cleaner::TauCleaner(double ptmin, double etamax){
310
311			std::vector<Tau> good_taus;
312			for(unsigned int i=0; i<bcc->taus->size(); ++i){
313			Tau tau = bcc->taus->at(i);
314	peiffer	1.5	if(tau.pt()>ptmin){
315			if(fabs(tau.eta())<etamax){
316			if(bcc->taus->at(i).decayModeFinding()){
317			if(bcc->taus->at(i).byMediumCombinedIsolationDeltaBetaCorr()){
318			if(bcc->taus->at(i).againstElectronTight()){
319			if(bcc->taus->at(i).againstMuonTight()){
320	peiffer	1.4	good_taus.push_back(tau);
321			}
322			}
323			}
324			}
325			}
326			}
327			}
328
329			bcc->taus->clear();
330
331			for(unsigned int i=0; i<good_taus.size(); ++i){
332			bcc->taus->push_back(good_taus[i]);
333			}
334	peiffer	1.5	sort(bcc->taus->begin(), bcc->taus->end(), HigherPt());
335	peiffer	1.9	resetEventCalc();
336	peiffer	1.4	}
337
338	peiffer	1.2
339	peiffer	1.3	void Cleaner::JetCleaner(double ptmin, double etamax, bool doPFID){
340
341			std::vector<Jet> good_jets;
342	peiffer	1.2	for(unsigned int i=0; i<bcc->jets->size(); ++i){
343			Jet jet = bcc->jets->at(i);
344	peiffer	1.5	if(jet.pt()>ptmin){
345			if(fabs(jet.eta())<etamax){
346	peiffer	1.2	if(!doPFID \|\| pfID(jet)){
347	peiffer	1.3	good_jets.push_back(jet);
348	peiffer	1.2	}
349			}
350			}
351			}
352
353	peiffer	1.3	bcc->jets->clear();
354
355			for(unsigned int i=0; i<good_jets.size(); ++i){
356			bcc->jets->push_back(good_jets[i]);
357			}
358	peiffer	1.5	sort(bcc->jets->begin(), bcc->jets->end(), HigherPt());
359	peiffer	1.9	resetEventCalc();
360	peiffer	1.2	}
361
362	peiffer	1.3	void Cleaner::TopJetCleaner(double ptmin, double etamax, bool doPFID){
363	peiffer	1.2
364	peiffer	1.3	std::vector<TopJet> good_topjets;
365	peiffer	1.2	for(unsigned int i=0; i<bcc->topjets->size(); ++i){
366			TopJet topjet = bcc->topjets->at(i);
367	peiffer	1.5	if(topjet.pt()>ptmin){
368			if(fabs(topjet.eta())<etamax){
369	peiffer	1.2	if(!doPFID \|\| pfID(topjet)){
370	peiffer	1.3	good_topjets.push_back(topjet);
371	peiffer	1.2	}
372			}
373			}
374			}
375	peiffer	1.3	bcc->topjets->clear();
376	peiffer	1.2
377	peiffer	1.3	for(unsigned int i=0; i<good_topjets.size(); ++i){
378			bcc->topjets->push_back(good_topjets[i]);
379			}
380	peiffer	1.5	sort(bcc->topjets->begin(), bcc->topjets->end(), HigherPt());
381	peiffer	1.9	resetEventCalc();
382	peiffer	1.2	}
383	peiffer	1.7
384
385			void Cleaner::PrimaryVertexCleaner(int ndofmax, double zmax, double rhomax){
386
387			std::vector<PrimaryVertex> good_pvs;
388			for(unsigned int i=0; i<bcc->pvs->size(); ++i){
389			PrimaryVertex pv = bcc->pvs->at(i);
390			if(pv.ndof()>=ndofmax && fabs(pv.z())<zmax && pv.rho()<=rhomax){
391			good_pvs.push_back(pv);
392			}
393			}
394			bcc->pvs->clear();
395
396			for(unsigned int i=0; i<good_pvs.size(); ++i){
397			bcc->pvs->push_back(good_pvs[i]);
398			}
399	peiffer	1.9	resetEventCalc();
400	peiffer	1.7	}