ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitHzz4l/Selection/src/SelectionFuncs.cc
Revision: 1.4
Committed: Wed Jun 20 17:40:22 2012 UTC (12 years, 10 months ago) by khahn
Content type: text/plain
Branch: MAIN
Changes since 1.3: +7 -5 lines
Log Message:
add 2012 offline SF and clean up

File Contents

# Content
1 #include "SelectionFuncs.h"
2
3 extern vector<SimpleLepton> failingLeptons;
4 extern vector<SimpleLepton> passingLeptons;
5
6 extern vector<unsigned> cutvec;
7 extern vector<vector<unsigned> > zcutvec;
8 extern vector<vector<unsigned> > zzcutvec;
9 extern map<unsigned,float> evtrhoMap;
10 extern bool passes_HLT_MC;
11
12 //----------------------------------------------------------------------------
13 bool setPV(ControlFlags ctrl,
14 const mithep::Array<mithep::Vertex> * vtxArr,
15 const mithep::Vertex* &vtx)
16 //----------------------------------------------------------------------------
17 {
18
19 const mithep::Vertex *bestPV = 0;
20 float best_sumpt=-1;
21
22 // good PV requirements
23 const UInt_t fMinNTracksFit = 0;
24 const Double_t fMinNdof = 4;
25 const Double_t fMaxAbsZ = 24;
26 const Double_t fMaxRho = 2;
27
28 for(int i=0; i<vtxArr->GetEntries(); ++i) {
29 const mithep::Vertex *pv = (mithep::Vertex*)(vtxArr->At(i));
30 if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl;
31
32 // Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used
33 if(!pv->IsValid()) continue;
34 if(pv->NTracksFit() < fMinNTracksFit) continue;
35 if(pv->Ndof() < fMinNdof) continue;
36 if(fabs(pv->Z()) > fMaxAbsZ) continue;
37 if(pv->Position().Rho() > fMaxRho) continue;
38
39 // take the first ...
40 bestPV = pv;
41 break;
42
43 // this never reached ...
44 float tmp_sumpt=0;
45 for( int t=0; t<pv->NTracks(); t++ )
46 tmp_sumpt += pv->Trk(t)->Pt();
47
48 if( tmp_sumpt > best_sumpt ) {
49 bestPV = pv;
50 best_sumpt = tmp_sumpt;
51 if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl;
52 }
53 }
54
55 // sync
56 if(!bestPV)
57 return false;
58 else {
59 vtx = bestPV;
60 return true;
61 }
62 }
63 //----------------------------------------------------------------------------------------
64 void writeEntries(ControlFlags ctrl, unsigned total_unskimmed)
65 {
66 ofstream ofs;
67 TString ofname(ctrl.outputfile);
68 ofs.open(ofname.ReplaceAll(".root",".nevents"));
69 assert(ofs.is_open());
70 ofs << total_unskimmed << endl;
71 ofs.close();
72 }
73 //----------------------------------------------------------------------------------------
74 void writeEntries(FOFlags foctrl, unsigned total_unskimmed)
75 {
76 ControlFlags ctrl;
77 ctrl.outputfile = foctrl.outputfile;
78 writeEntries(ctrl,total_unskimmed);
79 }
80 //----------------------------------------------------------------------------------------
81 void getEATargets(ControlFlags &ctrl, mithep::MuonTools::EMuonEffectiveAreaTarget &eraMu, mithep::ElectronTools::EElectronEffectiveAreaTarget &eraEle)
82 {
83 if( !ctrl.mc && ctrl.era == 2011 ) {
84 eraMu = mithep::MuonTools::kMuEAData2011;
85 eraEle = mithep::ElectronTools::kEleEAData2011;
86 } else if( !ctrl.mc && ctrl.era == 2012 ) {
87 eraMu = mithep::MuonTools::kMuEAData2012;
88 eraEle = mithep::ElectronTools::kEleEAData2012;
89 } else if( ctrl.mc && ctrl.era == 2011 ) {
90 eraMu = mithep::MuonTools::kMuEAFall11MC;
91 eraEle = mithep::ElectronTools::kEleEAFall11MC;
92 } else if( ctrl.mc && ctrl.era == 2012 ) {
93 eraMu = mithep::MuonTools::kMuEAData2012;
94 eraEle = mithep::ElectronTools::kEleEAData2012;
95 } else {
96 cerr << "unknown era for effective areas ... quitting." << endl;
97 exit(1);
98 }
99 }
100 //----------------------------------------------------------------------------------------
101 void getEATargets(FOFlags foctrl, mithep::MuonTools::EMuonEffectiveAreaTarget &eraMu, mithep::ElectronTools::EElectronEffectiveAreaTarget &eraEle)
102 {
103 ControlFlags ctrl;
104 ctrl.mc = foctrl.mc;
105 ctrl.era = foctrl.era;
106 getEATargets(ctrl,eraMu,eraEle);
107 }
108 //----------------------------------------------------------------------------
109 unsigned makePFnoPUArray(mithep::Array<PFCandidate> * fPFCandidates,
110 vector<bool> &pfNoPileUpflag,
111 const mithep::Array<mithep::Vertex> * vtxArr )
112 //----------------------------------------------------------------------------
113 {
114 for(UInt_t i = 0; i < fPFCandidates->GetEntries(); i++) {
115 const PFCandidate *pf = fPFCandidates->At(i);
116 assert(pf);
117
118 if(pf->PFType() == PFCandidate::eHadron) {
119 if(pf->HasTrackerTrk() &&
120 vtxArr->At(0)->HasTrack(pf->TrackerTrk()) &&
121 vtxArr->At(0)->TrackWeight(pf->TrackerTrk()) > 0) {
122
123 pfNoPileUpflag.push_back(1);
124
125 } else {
126
127 Bool_t vertexFound = kFALSE;
128 const Vertex *closestVtx = 0;
129 Double_t dzmin = 10000;
130
131 // loop over vertices
132 for(UInt_t j = 0; j < vtxArr->GetEntries(); j++) {
133 const Vertex *vtx = vtxArr->At(j);
134 assert(vtx);
135
136 if(pf->HasTrackerTrk() &&
137 vtx->HasTrack(pf->TrackerTrk()) &&
138 vtx->TrackWeight(pf->TrackerTrk()) > 0) {
139 vertexFound = kTRUE;
140 closestVtx = vtx;
141 break;
142 }
143 Double_t dz = fabs(pf->SourceVertex().Z() - vtx->Z());
144 if(dz < dzmin) {
145 closestVtx = vtx;
146 dzmin = dz;
147 }
148 }
149
150 // if(fCheckClosestZVertex) {
151 if(1) {
152 // Fallback: if track is not associated with any vertex,
153 // associate it with the vertex closest in z
154 if(vertexFound || closestVtx != vtxArr->At(0)) {
155 // pfPileUp->Add(pf);
156 pfNoPileUpflag.push_back(0);
157 } else {
158 pfNoPileUpflag.push_back(1);
159 }
160 } else {
161 if(vertexFound && closestVtx != vtxArr->At(0)) {
162 // pfPileUp->Add(pf);
163 pfNoPileUpflag.push_back(0);
164 } else {
165 // PFCandidate * pfNoPileUp->AddNew(); // Ridiculous but that's how it is
166 pfNoPileUpflag.push_back(1);
167 }
168 }
169 } //hadron & trk stuff
170 } else { // hadron
171 // PFCandidate * ptr = pfNoPileUp->AddNew();
172 pfNoPileUpflag.push_back(1);
173 }
174 } // Loop over PF candidates
175
176 return pfNoPileUpflag.size();
177 }
178
179 //--------------------------------------------------------------------------------
180 void setEffiencyWeights(unsigned era, EventData & evtdat, WeightStruct & weights )
181 //--------------------------------------------------------------------------------
182 {
183 vector<SimpleLepton> lepvec = evtdat.Z1leptons;
184 lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end());
185 double w_offline=-1, werr_offline=0;
186 double w_online=-1, werr_online=0;
187
188 vector< pair <double,double> > wlegs; // pair here is eff & err
189 vector< pair <float,float> > mvec; // pair here is eta & pt
190
191 // vector< pair <float,float> > evec; // pair here is eta & pt
192 // now deal with medium vs loose
193 vector< pair< bool, pair <float,float> > > evec; // pair here is eta & pt
194
195 for( int k=0; k<lepvec.size(); k++ ) {
196 if( abs(lepvec[k].type) == 13 ) {
197 mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
198 wlegs.push_back( muonPerLegOfflineEfficiencyWeight( era, fabs(lepvec[k].vec.Eta()),
199 lepvec[k].vec.Pt() ) );
200 } else {
201
202 // now deal with medium vs loose
203 // evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) );
204
205 std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt());
206 evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) );
207
208 // wlegs.push_back( elePerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()),
209 // lepvec[k].vec.Pt() ) );
210
211 wlegs.push_back( elePerLegOfflineEfficiencyWeight( era,
212 fabs(lepvec[k].vec.Eta()),
213 lepvec[k].vec.Pt(),
214 lepvec[k].isTight ) );
215
216 }
217 }
218
219 pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs );
220 weights.woff = offpair.first;
221 weights.werroff = offpair.second;
222
223 pair<double,double> onpair = getOnlineEfficiencyWeight( mvec, evec );
224 weights.won = onpair.first;
225 weights.werron = onpair.second;
226 }
227 //----------------------------------------------------------------------------
228 void initEvtRhoMap( map<unsigned, float> & evtrhoMap )
229 //----------------------------------------------------------------------------
230 {
231 unsigned evt;
232 float rho;
233
234 cout << "initialzing EvtRhoMap ";
235 //FILE * f = fopen("evtrho.dat", "r");
236 // FILE * f = fopen("allrho.cali.uniq.dat", "r");
237 // FILE * f = fopen("/data/blue/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
238 assert(fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat","r")); FILE * f = fopen("/scratch/dkralph/khahn/CMSSW_4_4_2/src/MitHzz4l/allrhoAA.cali.uniq.dat", "r");
239 int lcount=0;
240 while( fscanf( f, "%u:%f", &evt, &rho ) ) {
241 if( feof(f) ) break;
242 evtrhoMap[evt] = rho;
243 lcount++;
244 if( !(lcount%1000) ) cout << "."; flush(cout);
245 }
246 cout << " done" << endl;
247 }
248
249 void setEra(string fname, ControlFlags &ctrl)
250 {
251 if( ctrl.mc && (fname.find("f11-") != string::npos) ) ctrl.era=2011 ;
252 if( ctrl.mc && (fname.find("s12-") != string::npos) ) ctrl.era=2012 ;
253 if( !ctrl.mc && (fname.find("r11-") != string::npos) ) ctrl.era=2011 ;
254 if( !ctrl.mc && (fname.find("r12-") != string::npos) ) ctrl.era=2012 ;
255 }
256 //----------------------------------------------------------------------------------------
257 int getNPU(mithep::Array<mithep::PileupInfo> *puArr, int bx)
258 {
259
260 cout << "size of puarr: " << puArr->GetEntries() << endl;
261 // Identify bunch crossing 0 in PU info
262 int ibx0=-1, ibxm=-1, ibxp=-1;
263 cout<< "before: " << ibx0<< "\t" << ibxm<<"\t"<<ibxp<<endl;
264 assert(puArr);
265 for(unsigned i=0; i<puArr->GetEntries(); ++i) {
266 if(puArr->At(i)->GetBunchCrossing() == 0) ibx0=i;
267 if(puArr->At(i)->GetBunchCrossing() ==-1) ibxm=i;
268 if(puArr->At(i)->GetBunchCrossing() == 1) ibxp=i;
269 }
270 cout<< "after: " << ibx0<< "\t" << ibxm<<"\t"<<ibxp<<endl;
271
272 if(bx==0) return puArr->At(ibx0)->GetPU_NumInteractions();
273 else if(bx==-1) return puArr->At(ibxm)->GetPU_NumInteractions();
274 else if(bx== 1) return puArr->At(ibxp)->GetPU_NumInteractions();
275 else assert(0);
276
277 return -1;
278 }