11 |
|
// |
12 |
|
// root headers |
13 |
|
// |
14 |
< |
#include <TFile.h> |
14 |
> |
#include <TFile.h> |
15 |
> |
#include <TRegexp.h> |
16 |
|
#include <TTree.h> |
17 |
|
#include <TChain.h> |
18 |
|
#include <TBranch.h> |
27 |
|
#include "TMVA/Config.h" |
28 |
|
#include "TMVA/MethodBDT.h" |
29 |
|
|
30 |
+ |
#include "FactorizedJetCorrector.h" |
31 |
+ |
#include "JetCorrectorParameters.h" |
32 |
+ |
|
33 |
|
// |
34 |
|
// ntuple format headers |
35 |
|
// |
37 |
|
#include "Electron.h" |
38 |
|
#include "Muon.h" |
39 |
|
#include "Vertex.h" |
40 |
+ |
#include "PFMet.h" |
41 |
+ |
#include "PFJetCol.h" |
42 |
|
#include "PFCandidate.h" |
43 |
|
#include "PFCandidateCol.h" |
44 |
|
#include "PileupInfoCol.h" |
45 |
|
#include "PileupEnergyDensity.h" |
46 |
+ |
#include "MCParticle.h" |
47 |
|
#include "TriggerMask.h" |
48 |
|
#include "TriggerTable.h" |
49 |
+ |
#include "TriggerObject.h" |
50 |
+ |
#include "TriggerObjectRel.h" |
51 |
|
#include "Names.h" |
52 |
|
#include "BaseMod.h" |
53 |
+ |
#include "TriggerObjectsTable.h" |
54 |
+ |
#include "JetIDMVA.h" |
55 |
|
|
56 |
|
// |
57 |
|
// our headers |
59 |
|
#include "ParseArgs.h" |
60 |
|
#include "MuonSelection.h" |
61 |
|
#include "ElectronSelection.h" |
62 |
+ |
#include "ElectronMomentumCorrection.h" |
63 |
+ |
#include "JetSelection.h" |
64 |
|
#include "IsolationSelection.h" |
65 |
< |
|
66 |
< |
//#include "Selection.h" |
54 |
< |
//#include "ReferenceSelection.h" |
55 |
< |
#include "ReferenceSelection2.h" |
65 |
> |
#include "ReferenceSelection.h" |
66 |
> |
#include "MuonMomentumCorrection.h" |
67 |
|
|
68 |
|
#include "TriggerUtils.h" |
69 |
|
#include "PassHLT.h" |
59 |
– |
#include "Angles.h" |
70 |
|
#include "KinematicsStruct.h" |
71 |
|
#include "InfoStruct.h" |
72 |
< |
//#include "GenInfoStruct.h" |
72 |
> |
#include "GenInfoStruct.h" |
73 |
|
#include "WeightStruct.h" |
64 |
– |
//#include "GlobalDataAndFuncs.h" |
65 |
– |
#include "RunLumiRangeMap.h" |
66 |
– |
|
74 |
|
#include "AngleTuple.h" |
75 |
|
#include "FOTuple.h" |
76 |
|
|
77 |
+ |
#include "RunLumiRangeMap.h" |
78 |
|
#include "SampleWeight.h" |
79 |
|
#include "EfficiencyWeightsInterface.h" |
80 |
+ |
#include "SelectionFuncs.h" |
81 |
|
|
82 |
|
#include "SimpleLepton.h" |
83 |
|
|
84 |
+ |
#include "Cintex/Cintex.h" |
85 |
+ |
|
86 |
|
#ifndef CMSSW_BASE |
87 |
|
#define CMSSW_BASE "../" |
88 |
|
#endif |
92 |
|
// |
93 |
|
// globals |
94 |
|
//---------------------------------------------------------------------------- |
95 |
< |
TH1D * hpu; |
95 |
> |
TH1D *hpu_2011, *hpu_2012; |
96 |
|
RunLumiRangeMap rlrm; |
97 |
|
vector<SimpleLepton> failingLeptons ; // for fake estimation |
98 |
|
vector<SimpleLepton> passingLeptons; // for fake estimation |
100 |
|
vector<vector<unsigned> > zcutvec; |
101 |
|
vector<vector<unsigned> > zzcutvec; |
102 |
|
map<unsigned,float> evtrhoMap; |
103 |
+ |
TrigInfo ti; |
104 |
|
vector<string> cutstrs; |
105 |
< |
bool passes_HLT_MC; |
105 |
> |
bool passes_HLT; |
106 |
> |
vector<bool> PFnoPUflag;; |
107 |
> |
map<TString, map<TString,int>* > counts; |
108 |
> |
ElectronMomentumCorrection electron_momentum_correction; |
109 |
> |
MuCorr *muCorr; |
110 |
> |
void addDummyZ2Lepotons(EventData &ret4l); |
111 |
> |
void setNtupleVals(ControlFlags &ctrl, EventData &ret4l, InfoStruct &evtinfo, Angles &angles, KinematicsStruct &kine, |
112 |
> |
JetInfoStruct &ji, WeightStruct &weights, EventHeader *info, Array<Muon> *muonArr, Array<Electron> *electronArr, Array<PFJet> *jetArr, Array<PFMet> *metArr, |
113 |
> |
vector<SimpleLepton> &simpleLeptonJets, FactorizedJetCorrector *fJetCorrector, Array<PileupInfo> *puArr, |
114 |
> |
Array<PileupEnergyDensity> *puDArr, JetIDMVA *fJetIDMVA, Array<Vertex> *vtxArr, bitset<1024> &triggerBits, |
115 |
> |
TriggerTable *hltTable, Array<TriggerObject> *hltObjArr, TriggerObjectsTable *fTrigObjs, AngleTuple *nt); |
116 |
> |
//---------------------------------------------------------------------------- |
117 |
> |
void initRunLumiRangeMap(ControlFlags &ctrl) |
118 |
> |
//---------------------------------------------------------------------------- |
119 |
> |
{ |
120 |
> |
if(ctrl.jsonFile==""){ |
121 |
> |
cout << "WARNING: running with empty jsonFile, so setting noJSON to true" << endl; |
122 |
> |
ctrl.noJSON = true; |
123 |
> |
} |
124 |
> |
rlrm.AddJSONFile(string(ctrl.jsonFile)); |
125 |
> |
} |
126 |
|
|
127 |
|
// |
128 |
|
// prototypes |
129 |
|
//---------------------------------------------------------------------------- |
98 |
– |
bool setPV(ControlFlags,const mithep::Array<mithep::Vertex>*, mithep::Vertex& ); |
99 |
– |
void initPUWeights(); |
100 |
– |
double getPUWeight(unsigned npu); |
101 |
– |
void setEffiencyWeights(EventData&, WeightStruct& ); |
102 |
– |
void initRunLumiRangeMap(); |
103 |
– |
void initEvtRhoMap(map<unsigned,float> &); |
104 |
– |
//---------------------------------------------------------------------------- |
105 |
– |
|
130 |
|
|
131 |
|
// |
132 |
|
// MAIN |
134 |
|
int main(int argc, char** argv) |
135 |
|
//---------------------------------------------------------------------------- |
136 |
|
{ |
113 |
– |
cutstrs.push_back(string("skim0")); //0 |
114 |
– |
cutstrs.push_back(string("skim1")); //1 |
115 |
– |
cutstrs.push_back(string("trigger")); //2 |
116 |
– |
// ------------------------------------------------- |
117 |
– |
cutstrs.push_back(string("Z candidate")); //3 |
118 |
– |
cutstrs.push_back(string("good Z1")); //4 |
119 |
– |
// ------------------------------------------------- |
120 |
– |
cutstrs.push_back(string("4l")); //5 |
121 |
– |
cutstrs.push_back(string("ZZ candidate")); //6 |
122 |
– |
cutstrs.push_back(string("good Z2")); //7 |
123 |
– |
cutstrs.push_back(string("ZZ 20/10")); //8 |
124 |
– |
cutstrs.push_back(string("resonance")); //9 |
125 |
– |
cutstrs.push_back(string("m4l>70")); //10 |
126 |
– |
cutstrs.push_back(string("m4l>100")); //11 |
127 |
– |
|
128 |
– |
|
129 |
– |
|
137 |
|
|
138 |
+ |
vector<TString> cuts; |
139 |
+ |
cuts.push_back("start"); // NOTE: is incremented in ReferenceSelection, so if you turn on JSONS be careful |
140 |
+ |
cuts.push_back("trigger"); |
141 |
+ |
cuts.push_back("sfOsHiPt"); |
142 |
+ |
cuts.push_back("4<mZ2<120"); |
143 |
+ |
cuts.push_back("pt>20,10"); |
144 |
+ |
cuts.push_back("mll>4"); |
145 |
+ |
cuts.push_back("m4l>70"); |
146 |
+ |
cuts.push_back("mZ2>12"); |
147 |
+ |
cuts.push_back("m4l>100"); |
148 |
+ |
cuts.push_back("m4l > 100 && >= 1 jets"); |
149 |
+ |
cuts.push_back("m4l > 100 && 2 jets"); |
150 |
+ |
cuts.push_back("m4l > 100 && fisher > 0.4"); |
151 |
+ |
|
152 |
+ |
for(unsigned icut=0; icut<cuts.size(); icut++) |
153 |
+ |
counts[cuts[icut]] = new map<TString,int>; |
154 |
+ |
|
155 |
+ |
map<TString,map<TString,int>* >::iterator cit; |
156 |
+ |
for(cit=counts.begin(); cit!=counts.end(); cit++) { |
157 |
+ |
map<TString,int> *cts = (*cit).second; |
158 |
+ |
(*cts)["all"] = int(0); |
159 |
+ |
(*cts)["4e"] = int(0); |
160 |
+ |
(*cts)["4m"] = int(0); |
161 |
+ |
(*cts)["2e2m"] = int(0); |
162 |
+ |
} |
163 |
+ |
|
164 |
|
string cmsswpath(CMSSW_BASE); |
165 |
|
cmsswpath.append("/src"); |
166 |
< |
std::bitset<1024> triggerBits; |
166 |
> |
bitset<1024> triggerBits; |
167 |
|
vector<vector<string> > inputFiles; |
168 |
|
inputFiles.push_back(vector<string>()); |
169 |
|
map<string,unsigned> entrymap; // number of unskimmed entries in each file |
137 |
– |
for( int i=0; i<cutstrs.size(); i++ ) cutvec.push_back(0); |
138 |
– |
for( int i=0; i<2; i++ ) { |
139 |
– |
zcutvec.push_back(vector<unsigned>()); |
140 |
– |
for( int j=0; j<cutstrs.size(); j++ ) zcutvec[i].push_back(0); |
141 |
– |
} |
142 |
– |
for( int i=0; i<3; i++ ) { |
143 |
– |
zzcutvec.push_back(vector<unsigned>()); |
144 |
– |
for( int j=0; j<cutstrs.size(); j++ ) zzcutvec[i].push_back(0); |
145 |
– |
} |
146 |
– |
|
170 |
|
|
171 |
|
// |
172 |
|
// args |
173 |
|
//-------------------------------------------------------------------------------------------------------------- |
174 |
|
ControlFlags ctrl; |
175 |
|
parse_args( argc, argv, ctrl ); |
153 |
– |
if( ctrl.inputfiles.empty() &&ctrl.inputfile.empty() ) |
154 |
– |
{ |
155 |
– |
cerr << "usage: applySelection.exe <flags> " << endl; |
156 |
– |
cerr << "\tmandoatory flags : " << endl; |
157 |
– |
cerr << "\t--inputfiles | file containing a list of ntuples to run over" << endl; |
158 |
– |
cerr << "\t--inputfile | a file to run over" << endl; |
159 |
– |
cerr << "\t--outputfile | file to store selected evet" << endl; |
160 |
– |
return 1; |
161 |
– |
} |
176 |
|
ctrl.dump(); |
177 |
< |
|
178 |
< |
|
177 |
> |
cout << "\n\nNOTE: now setting era with arguments... double check this doesnt mess something else up\n\n" << endl; |
178 |
> |
assert(ctrl.doFSR); |
179 |
|
|
180 |
|
// |
181 |
|
// File I/O |
182 |
|
//-------------------------------------------------------------------------------------------------------------- |
183 |
|
TChain * chain = new TChain("Events"); |
184 |
|
TChain * hltchain = new TChain("HLT"); |
185 |
+ |
TChain * runchain = new TChain(Names::gkRunTreeName); |
186 |
|
|
187 |
|
string fname; |
188 |
+ |
unsigned total_unskimmed=0; |
189 |
|
if( !ctrl.inputfiles.empty() ) { |
190 |
|
ifstream f(ctrl.inputfiles.c_str()); |
191 |
< |
while (f >> fname) { |
191 |
> |
while (f >> fname) { |
192 |
|
if( !(strncmp( fname.c_str(), "#", 1 ) ) ) continue; // skip commented lines |
193 |
+ |
if(fname.find("/store/") != string::npos) fname = "root://eoscms//"+fname; |
194 |
|
cout << "adding inputfile : " << fname.c_str() << endl; |
195 |
|
entrymap[string(fname.c_str())] = unskimmedEntries(fname.c_str()); |
196 |
< |
cout << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl; |
196 |
> |
cerr << "unskimmed entries: " << entrymap[string(fname.c_str())] << endl; |
197 |
> |
total_unskimmed += entrymap[string(fname.c_str())]; |
198 |
|
chain->AddFile(fname.c_str()); |
199 |
|
hltchain->AddFile(fname.c_str()); |
200 |
+ |
runchain->AddFile(fname.c_str()); |
201 |
|
} |
202 |
|
} else { |
203 |
+ |
if(ctrl.inputfile.find("/store/") != string::npos) ctrl.inputfile = "root://eoscms//"+ctrl.inputfile; |
204 |
|
cout << "adding inputfile : " << ctrl.inputfile.c_str() << endl; |
205 |
< |
unsigned tmpent = unskimmedEntries(ctrl.inputfile.c_str()); |
206 |
< |
cout << "tmpent: " << tmpent << endl; |
207 |
< |
entrymap[string(ctrl.inputfile.c_str())] = unskimmedEntries(ctrl.inputfile.c_str()); |
208 |
< |
cout << "unskimmed entries: " << entrymap[string(ctrl.inputfile.c_str())] << endl; |
205 |
> |
unsigned unsk_ents = unskimmedEntries(ctrl.inputfile.c_str()); |
206 |
> |
entrymap[string(ctrl.inputfile.c_str())] = unsk_ents; |
207 |
> |
cerr << "unskimmed entries: " << unsk_ents << endl; |
208 |
> |
total_unskimmed += unsk_ents; |
209 |
|
chain->AddFile(ctrl.inputfile.c_str()); |
210 |
|
hltchain->AddFile(ctrl.inputfile.c_str()); |
211 |
+ |
runchain->AddFile(ctrl.inputfile.c_str()); |
212 |
|
} |
213 |
+ |
// write the total number of unskimmed events that went into making this output file to a text file |
214 |
+ |
writeEntries(ctrl,total_unskimmed); |
215 |
+ |
ctrl.totalMC = total_unskimmed; |
216 |
|
|
217 |
|
const char * ofname; |
218 |
|
if( strcmp( ctrl.outputfile.c_str(), "") ) { |
224 |
|
string xspath = (cmsswpath+string("/MitPhysics/data/xs.dat")); |
225 |
|
cout << "xspath: " << xspath.c_str() << endl; |
226 |
|
SimpleTable xstab(xspath.c_str()); |
227 |
+ |
if(ctrl.mc) { |
228 |
+ |
if(ctrl.use_xs_weights) getWeight(xstab,entrymap,chain); // test *now* whether we've got the xs in here, so it fails before it creates the output |
229 |
+ |
} |
230 |
|
|
231 |
+ |
//root will write to this temporary file instead of to memory and the selection will run faster |
232 |
|
|
233 |
|
// |
234 |
|
// Setup |
235 |
|
//-------------------------------------------------------------------------------------------------------------- |
236 |
|
Angles angles; |
237 |
< |
KinematicsStruct kinematics; |
237 |
> |
KinematicsStruct kine; |
238 |
|
InfoStruct evtinfo; |
239 |
|
WeightStruct weights; |
240 |
< |
// GenInfoStruct geninfo; |
240 |
> |
GenInfoStruct geninfo; |
241 |
> |
JetInfoStruct ji; |
242 |
> |
vector<SimpleLepton> simpleLeptonJets; |
243 |
|
|
244 |
|
AngleTuple nt( (const char*)ofname, (const char*)"zznt"); |
245 |
|
nt.makeAngleBranch(angles); |
246 |
< |
nt.makeKinematicsBranch(kinematics); |
246 |
> |
nt.makeKinematicsBranch(kine); |
247 |
|
nt.makeInfoBranch(evtinfo); |
248 |
+ |
nt.makeJetInfoBranch(ji); |
249 |
+ |
//nt.makeJetAngleBranches("DataStruct/data/jetAngleVars.txt"); |
250 |
|
|
251 |
< |
FOTuple foTree( nt.getFile(), (const char*)"FOtree"); |
252 |
< |
foTree.makeFailedLeptonBranch(failingLeptons); |
253 |
< |
foTree.makeZLeptonBranch(passingLeptons); |
251 |
> |
FOTuple *foTree=0; |
252 |
> |
if(ctrl.fakeScheme != "none") { |
253 |
> |
foTree = new FOTuple( nt.getFile(), (const char*)"FOtree"); |
254 |
> |
foTree->makeInfoBranch(evtinfo); |
255 |
> |
foTree->makeFailedLeptonBranch(failingLeptons); |
256 |
> |
foTree->makePassedLeptonBranch(passingLeptons); |
257 |
> |
} |
258 |
|
|
259 |
< |
TH1F * h_w_hpt; |
259 |
> |
int genSampleType(0); |
260 |
|
if(ctrl.mc) { |
261 |
|
nt.makeWeightBranch(weights); |
262 |
< |
// nt.makeGenInfoBranch(geninfo); |
263 |
< |
initEfficiencyWeights(); |
264 |
< |
initPUWeights(); |
265 |
< |
/* |
266 |
< |
// Higgs only, pt reweighting |
267 |
< |
if( ctrl.mH <= 140 ) { |
268 |
< |
char wptname[256]; |
233 |
< |
sprintf( wptname, "/data/blue/pharris/Flat/ntupler/root/weight_ptH_%i.root", ctrl.mH ); |
234 |
< |
TFile * f_hpt = new TFile(wptname); |
235 |
< |
f_hpt->Print(); |
236 |
< |
sprintf(wptname, "weight_hqt_fehipro_fit_%i", ctrl.mH); |
237 |
< |
h_w_hpt = (TH1F*)(f_hpt->FindObjectAny(wptname)); |
238 |
< |
h_w_hpt->Print(); |
262 |
> |
if(ctrl.fillGen) { |
263 |
> |
genSampleType = getGenSampleType(ofname); |
264 |
> |
if(genSampleType == -1) { |
265 |
> |
cout << "\ngenSampleType returned as -1, setting fillGen to false" << endl; |
266 |
> |
ctrl.fillGen = false; |
267 |
> |
} else |
268 |
> |
nt.makeGenInfoBranch(geninfo); |
269 |
|
} |
270 |
< |
*/ |
270 |
> |
if(ctrl.use_eff_weights) initEfficiencyWeights(); |
271 |
|
} else { |
272 |
< |
initRunLumiRangeMap(); |
272 |
> |
if(!(ctrl.noJSON) ) |
273 |
> |
initRunLumiRangeMap(ctrl); |
274 |
|
} |
275 |
|
|
276 |
< |
// initMuonIDMVA(); |
246 |
< |
initMuonIsoMVA(); |
247 |
< |
initElectronIDMVA(); |
248 |
< |
initElectronIsoMVA(); |
276 |
> |
initElectronIDMVAV1(); |
277 |
|
initTrigger(); |
278 |
+ |
initAnalysisTriggers(ti); |
279 |
+ |
initMela(ctrl); |
280 |
+ |
initMassErrors(); |
281 |
+ |
ti.dump(); |
282 |
+ |
|
283 |
+ |
TFile tmp_file("tmp_file.root"); |
284 |
+ |
|
285 |
+ |
if(ctrl.do_eregression) { |
286 |
+ |
string regPath("EGamma/EGammaAnalysisTools/data/eleEnergyRegWeights_V1.root"); |
287 |
+ |
// if(TString(getenv("HOSTNAME")).Contains(TRegexp("t3[bd][te][cs][hk]"))) |
288 |
+ |
regPath = cmsswpath + "/" + regPath; |
289 |
+ |
electron_momentum_correction.initialize_regression(regPath); |
290 |
+ |
} |
291 |
|
|
292 |
< |
// hack |
293 |
< |
initEvtRhoMap(evtrhoMap); |
292 |
> |
if(ctrl.do_escale) { |
293 |
> |
if(ctrl.mc) { |
294 |
> |
if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Fall11"; |
295 |
> |
else ctrl.smearing_scaling_dataset = "Summer12_DR53X_HCP2012"; |
296 |
> |
} else { |
297 |
> |
if(ctrl.era == 2011) ctrl.smearing_scaling_dataset = "Jan16ReReco"; |
298 |
> |
else ctrl.smearing_scaling_dataset = "Moriond2013"; |
299 |
> |
} |
300 |
> |
} |
301 |
> |
if(ctrl.correct_muon_momentum) { |
302 |
> |
muCorr = new MuCorr; |
303 |
> |
muCorr->corr2011 = new rochcor; |
304 |
> |
muCorr->corr2012 = new rochcor2012; |
305 |
> |
} |
306 |
> |
|
307 |
> |
vector<TString> fJetCorrParsv; |
308 |
> |
vector<JetCorrectorParameters> correctionParameters; |
309 |
> |
FactorizedJetCorrector *fJetCorrector; |
310 |
> |
JetIDMVA *fJetIDMVA; |
311 |
> |
initJets(ctrl, cmsswpath, fJetCorrParsv, correctionParameters, fJetCorrector, fJetIDMVA); |
312 |
|
|
254 |
– |
|
313 |
|
// |
314 |
|
// Setup tree I/O |
315 |
|
//-------------------------------------------------------------------------------------------------------------- |
317 |
|
TTree *eventTree=0; |
318 |
|
string currentFile(""); |
319 |
|
|
320 |
< |
mithep::EventHeader *info = new mithep::EventHeader(); |
321 |
< |
// mithep::TGenInfo *ginfo = new mithep::TGenInfo(); |
322 |
< |
mithep::Array<mithep::Electron> *electronArr = new mithep::Array<mithep::Electron>(); |
323 |
< |
mithep::Array<mithep::Muon> *muonArr = new mithep::Array<mithep::Muon>(); |
324 |
< |
mithep::Array<mithep::Vertex> *vtxArr = new mithep::Array<mithep::Vertex>(); |
325 |
< |
mithep::Array<mithep::PFCandidate> *pfArr = new mithep::Array<mithep::PFCandidate>(); |
326 |
< |
mithep::Array<mithep::PileupInfo> *puArr = new mithep::Array<mithep::PileupInfo>(); |
327 |
< |
mithep::Array<mithep::PileupEnergyDensity> *puDArr = new mithep::Array<mithep::PileupEnergyDensity>(); |
328 |
< |
mithep::Array<mithep::Track> *trkArr = new mithep::Array<mithep::Track>(); |
329 |
< |
mithep::TriggerMask *trigMask = new mithep::TriggerMask(); |
330 |
< |
mithep::TriggerTable *hltTable = new mithep::TriggerTable(); |
331 |
< |
vector<string> *hltTableStrings = new vector<string>(); |
332 |
< |
|
333 |
< |
TString fElectronName(Names::gkElectronBrn); |
334 |
< |
TString fMuonName(Names::gkMuonBrn); |
335 |
< |
TString fInfoName(Names::gkEvtHeaderBrn); |
336 |
< |
TString fPrimVtxName(Names::gkPVBrn); |
337 |
< |
TString fPFCandidateName(Names::gkPFCandidatesBrn); |
338 |
< |
TString fPileupInfoName(Names::gkPileupInfoBrn); |
339 |
< |
TString fPileupEnergyDensityName(Names::gkPileupEnergyDensityBrn); |
340 |
< |
TString fTrackName(Names::gkTrackBrn); |
283 |
< |
TString fTriggerMaskName(Names::gkHltBitBrn); |
284 |
< |
TString fTriggerTableName(Names::gkHltTableBrn); |
320 |
> |
UInt_t fNMaxTriggers = TRIGGER_BIG_NUMBER; |
321 |
> |
EventHeader *info = new EventHeader(); |
322 |
> |
Array<PFMet> *metArr = new Array<PFMet>(); |
323 |
> |
Array<Electron> *electronArr = new Array<Electron>(); |
324 |
> |
Array<Muon> *muonArr = new Array<Muon>(); |
325 |
> |
Array<Vertex> *vtxArr = new Array<Vertex>(); |
326 |
> |
Array<PFCandidate> *pfArr = new Array<PFCandidate>(); |
327 |
> |
Array<PFJet> *jetArr = new Array<PFJet>(); |
328 |
> |
Array<PileupInfo> *puArr = new Array<PileupInfo>(); |
329 |
> |
Array<PileupEnergyDensity> *puDArr = new Array<PileupEnergyDensity>(); |
330 |
> |
Array<Track> *trkArr = new Array<Track>(); |
331 |
> |
Array<MCParticle> *mcArr = new Array<MCParticle>(); |
332 |
> |
MCEventInfo *mcEvtInfo = new MCEventInfo(); |
333 |
> |
TriggerMask *trigMask = new TriggerMask(); |
334 |
> |
TriggerTable *hltTable = new TriggerTable(fNMaxTriggers); |
335 |
> |
vector<string> *hltTableStrings = new vector<string>(); |
336 |
> |
vector<string> *hltLabelStrings = new vector<string>(); |
337 |
> |
Array<TriggerObject> *hltObjArr = new Array<TriggerObject>(); |
338 |
> |
Array<TriggerObjectRel> *hltRelsArr = new Array<TriggerObjectRel>(); |
339 |
> |
TriggerObjectsTable *fTrigObjs = new TriggerObjectsTable(hltTable,fNMaxTriggers); |
340 |
> |
RunInfo *runInfo = new RunInfo(); |
341 |
|
|
342 |
< |
chain->SetBranchAddress(fInfoName, &info); |
343 |
< |
chain->SetBranchAddress(fElectronName, &electronArr); |
344 |
< |
chain->SetBranchAddress(fMuonName, &muonArr); |
345 |
< |
chain->SetBranchAddress(fPrimVtxName, &vtxArr); |
346 |
< |
chain->SetBranchAddress(fPFCandidateName, &pfArr); |
347 |
< |
if( ctrl.mc ) chain->SetBranchAddress(fPileupInfoName, &puArr); |
348 |
< |
chain->SetBranchAddress(fPileupEnergyDensityName, &puDArr); |
349 |
< |
chain->SetBranchAddress(fTrackName, &trkArr); |
350 |
< |
chain->SetBranchAddress(fTriggerMaskName, &trigMask); |
351 |
< |
cout << "hlttable: " << fTriggerTableName << endl; |
352 |
< |
hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings); |
353 |
< |
|
354 |
< |
mithep::Vertex vtx; // best primary vertex in the event |
355 |
< |
|
356 |
< |
// ginfo = NULL; |
357 |
< |
// if( ctrl.mc ) { chain->SetBranchAddress("Gen", &ginfo);} |
358 |
< |
|
359 |
< |
int count=0, pass=0; |
360 |
< |
float passcorr=0., passcorr_errup=0., passcorr_errdown=0.; |
305 |
< |
float denom[3]={0.,0.,0.}; |
306 |
< |
float numer[3]={0.,0.,0.}; |
307 |
< |
float numercorr[3]={0.,0.,0.}; |
342 |
> |
chain->SetBranchAddress(Names::gkEvtHeaderBrn, &info); |
343 |
> |
chain->SetBranchAddress("PFMet", &metArr); |
344 |
> |
chain->SetBranchAddress(Names::gkElectronBrn, &electronArr); |
345 |
> |
chain->SetBranchAddress(Names::gkMuonBrn, &muonArr); |
346 |
> |
chain->SetBranchAddress(Names::gkPVBrn, &vtxArr); |
347 |
> |
chain->SetBranchAddress(Names::gkPFCandidatesBrn, &pfArr); |
348 |
> |
chain->SetBranchAddress(Names::gkPFJetBrn, &jetArr); |
349 |
> |
if( ctrl.mc ) { |
350 |
> |
chain->SetBranchAddress(Names::gkPileupInfoBrn, &puArr); |
351 |
> |
chain->SetBranchAddress(Names::gkMCPartBrn, &mcArr); |
352 |
> |
chain->SetBranchAddress(Names::gkMCEvtInfoBrn, &mcEvtInfo); |
353 |
> |
} |
354 |
> |
chain->SetBranchAddress(Names::gkPileupEnergyDensityBrn, &puDArr); |
355 |
> |
chain->SetBranchAddress(Names::gkTrackBrn, &trkArr); |
356 |
> |
chain->SetBranchAddress(Names::gkHltBitBrn, &trigMask); |
357 |
> |
chain->SetBranchAddress(Names::gkHltObjBrn, &hltObjArr); |
358 |
> |
chain->SetBranchAddress("HLTObjectsRelation", &hltRelsArr); |
359 |
> |
hltchain->SetBranchAddress(Names::gkHltTableBrn, &hltTableStrings); |
360 |
> |
hltchain->SetBranchAddress(Names::gkHltLabelBrn, &hltLabelStrings); |
361 |
|
|
362 |
< |
// only 1 HLT table / file ??? |
363 |
< |
hltchain->GetEntry(0); |
364 |
< |
cerr << "here... " << endl; |
362 |
> |
runchain->SetBranchAddress(Names::gkRunInfoBrn, &runInfo); |
363 |
> |
|
364 |
> |
Vertex vtx; // best primary vertex in the event |
365 |
|
|
366 |
< |
int imax = chain->GetEntries(); |
314 |
< |
cout << "nEntries: " << imax << endl; |
366 |
> |
hltchain->GetEntry(0); |
367 |
|
|
368 |
+ |
int imax = (ctrl.n_events_to_process < 0) ? chain->GetEntries() : ctrl.n_events_to_process; |
369 |
+ |
cerr << "processing: " << imax << " ( / " << chain->GetEntries() << " )" << endl; |
370 |
|
|
371 |
< |
// |
372 |
< |
// Loop !!!!!!!!! |
319 |
< |
//-------------------------------------------------------------------------------------------------------------- |
320 |
< |
for(UInt_t ientry=0; ientry<imax; ientry++) |
321 |
< |
{ |
322 |
< |
chain->GetEntry(ientry); |
323 |
< |
if(!(ientry%1000)) cerr << "entry: " << ientry << endl; |
324 |
< |
setPV( ctrl, vtxArr, vtx); |
371 |
> |
if(ctrl.mc && ctrl.debugOffMc) // necessary for batch mode -- we need debug output for, data but size is too large to store MC output |
372 |
> |
ctrl.debug = false; |
373 |
|
|
374 |
< |
cout << "origrho: " << puDArr->At(0)->Rho() |
375 |
< |
<< "\torigrholoweta: " << puDArr->At(0)->RhoLowEta() |
374 |
> |
int lastHltEntry=-2; |
375 |
> |
for(UInt_t ientry=0; ientry<imax; ientry++) { |
376 |
> |
chain->GetEntry(ientry); |
377 |
> |
if(!(ientry%1000)) cerr << "entry: " << ientry << " ( / " << imax << " )" << endl; |
378 |
> |
|
379 |
> |
//if(info->EvtNum() != 217502611) continue; |
380 |
> |
|
381 |
> |
increment(counts,"start"); |
382 |
> |
if( ctrl.debug ) { |
383 |
> |
cout << "-----------------------------------------------------------------" << endl; |
384 |
> |
cout << "-----------------------------------------------------------------" << endl; |
385 |
> |
cout << "Run: " << info->RunNum() |
386 |
> |
<< "\tEvt: " << info->EvtNum() |
387 |
> |
<< "\tLumi: " << info->LumiSec() |
388 |
> |
<< endl; |
389 |
> |
cerr << "Run: " << info->RunNum() |
390 |
> |
<< "\tEvt: " << info->EvtNum() |
391 |
> |
<< "\tLumi: " << info->LumiSec() |
392 |
|
<< endl; |
393 |
+ |
cout << "-----------------------------------------------------------------" << endl; |
394 |
+ |
} |
395 |
|
|
396 |
+ |
runchain->GetEvent(info->RunEntry()); |
397 |
+ |
hltchain->GetEntry(runInfo->HltEntry()); |
398 |
+ |
hltTable->Clear(); |
399 |
+ |
fillTriggerNames(hltTable, hltTableStrings ); |
400 |
+ |
// if(ctrl.debug && (lastHltEntry != runInfo->HltEntry())) hltTable->Print(); |
401 |
+ |
lastHltEntry = runInfo->HltEntry(); |
402 |
+ |
fillTriggerBits( hltTable, trigMask, triggerBits ); |
403 |
+ |
setHLTObjectRelations( hltObjArr, hltRelsArr, hltTableStrings, hltLabelStrings, fTrigObjs ); |
404 |
+ |
// printTriggerObjs( hltTable, fTrigObjs); |
405 |
+ |
|
406 |
+ |
string fname = string(chain->GetFile()->GetEndpointUrl()->GetFile()); |
407 |
+ |
if(ctrl.debug) cout << "era is " << ctrl.era << endl; |
408 |
+ |
if( ctrl.era == 0 ) { |
409 |
+ |
setEra( fname, ctrl ); |
410 |
+ |
if(ctrl.debug) cout << "era is now " << ctrl.era << endl; |
411 |
+ |
} |
412 |
|
|
413 |
< |
// |
414 |
< |
// data/MC |
415 |
< |
// |
416 |
< |
if(ctrl.mc) { |
417 |
< |
// |
418 |
< |
// xsec & PU weights |
419 |
< |
// |
420 |
< |
// weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC; |
421 |
< |
int npu = -1; |
422 |
< |
for(int i=0;i<puArr->GetEntries();i++) { |
423 |
< |
if(puArr->At(i)->GetBunchCrossing() == 0) |
424 |
< |
npu = puArr->At(i)->GetPU_NumInteractions(); |
425 |
< |
} |
426 |
< |
assert(npu>=0); |
427 |
< |
evtinfo.npu; |
428 |
< |
weights.npuw = getPUWeight(evtinfo.npu); |
429 |
< |
// cout << "weight: " << weights.w << endl; |
430 |
< |
// |
431 |
< |
// trigger |
432 |
< |
// |
433 |
< |
if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) { |
434 |
< |
currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile()); |
435 |
< |
hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings); |
436 |
< |
hltchain->GetEntry(0); |
437 |
< |
hltTable->Clear(); |
356 |
< |
fillTriggerBits(hltTable, hltTableStrings ); |
357 |
< |
} |
358 |
< |
if( ctrl.debug ) cout << "file is : " << currentFile << endl; |
359 |
< |
fillTriggerBits( hltTable, trigMask, triggerBits ); |
360 |
< |
passes_HLT_MC = true; // passed to apply as extern global, just for sync |
361 |
< |
if( !passHLTMC(triggerBits, info->RunNum(), info->EvtNum() ) ) { |
362 |
< |
passes_HLT_MC = false; |
363 |
< |
} |
364 |
< |
} else { |
365 |
< |
// |
366 |
< |
// JSON |
367 |
< |
// |
368 |
< |
/* |
413 |
> |
// pfNoPU |
414 |
> |
PFnoPUflag.clear(); |
415 |
> |
int pfnopu_size = makePFnoPUArray( pfArr, PFnoPUflag, vtxArr ); |
416 |
> |
assert(pfnopu_size == pfArr->GetEntries()); |
417 |
> |
|
418 |
> |
// trigger |
419 |
> |
if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) { |
420 |
> |
currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile()); |
421 |
> |
} |
422 |
> |
if( ctrl.debug ) cout << "file is : " << currentFile << endl; |
423 |
> |
passes_HLT = passHLT(ctrl, triggerBits, ti, info->RunNum(), hltTable, fTrigObjs); // NOTE: run value is not used if it's MC |
424 |
> |
|
425 |
> |
// data/MC |
426 |
> |
if(ctrl.mc) { |
427 |
> |
if(ctrl.fillGen) fillGenInfo( mcArr, mcEvtInfo, geninfo, genSampleType, ctrl); |
428 |
> |
// NOTE: xs weight and npuw are reset in merge.cc |
429 |
> |
if(ctrl.use_xs_weights) |
430 |
> |
weights.w = getWeight(xstab,entrymap,chain)/ctrl.totalMC; |
431 |
> |
else |
432 |
> |
weights.w = 1; |
433 |
> |
weights.npu = getNPU(puArr, 0); |
434 |
> |
weights.npuw = getPUWeight(ctrl.era, fname, weights.npu); |
435 |
> |
} else { |
436 |
> |
// JSON |
437 |
> |
if(!(ctrl.noJSON) ) { |
438 |
|
RunLumiRangeMap::RunLumiPairType rl(info->RunNum(), info->LumiSec()); |
439 |
|
if( !(rlrm.HasRunLumi(rl)) ) { |
440 |
< |
if( ctrl.debug ) cout << "\tfails JSON" << endl; |
440 |
> |
if( ctrl.debug ) cout << "fails JSON" << endl; |
441 |
|
continue; |
442 |
|
} |
374 |
– |
*/ |
375 |
– |
// |
376 |
– |
// trigger |
377 |
– |
// |
378 |
– |
if( string(chain->GetFile()->GetEndpointUrl()->GetFile()) != currentFile ) { |
379 |
– |
currentFile = string(chain->GetFile()->GetEndpointUrl()->GetFile()); |
380 |
– |
hltchain->SetBranchAddress(fTriggerTableName, &hltTableStrings); |
381 |
– |
hltchain->GetEntry(0); |
382 |
– |
hltTable->Clear(); |
383 |
– |
fillTriggerBits(hltTable, hltTableStrings ); |
384 |
– |
} |
385 |
– |
if( ctrl.debug ) cout << "file is : " << currentFile << endl; |
386 |
– |
/* |
387 |
– |
fillTriggerBits( hltTable, trigMask, triggerBits ); |
388 |
– |
if( !passHLT(triggerBits, info->RunNum(), info->EvtNum() ) ) { |
389 |
– |
if( ctrl.debug ) cout << "\tfails trigger ... " << endl; |
390 |
– |
continue; |
391 |
– |
} |
392 |
– |
*/ |
443 |
|
} |
444 |
+ |
} |
445 |
|
|
446 |
< |
// |
447 |
< |
// lepton/kinematic selection ... |
448 |
< |
// |
449 |
< |
failingLeptons.clear(); |
450 |
< |
passingLeptons.clear(); |
451 |
< |
EventData ret4l = |
452 |
< |
|
453 |
< |
/* |
454 |
< |
apply_HZZ4L_reference2_selection(ctrl, info, |
455 |
< |
vtx, |
456 |
< |
pfArr, |
457 |
< |
puDArr, |
458 |
< |
electronArr, |
459 |
< |
&electronReferencePreSelection, |
460 |
< |
&electronReferenceIDMVASelection, |
461 |
< |
&electronReferenceIsoSelection, |
462 |
< |
muonArr, |
463 |
< |
&muonReferencePreSelection, |
464 |
< |
&muonIDPFSelection, |
414 |
< |
&muonReferenceIsoSelection); |
415 |
< |
*/ |
416 |
< |
apply_HZZ4L_reference2_selection(ctrl, info, |
417 |
< |
vtx, |
418 |
< |
pfArr, |
419 |
< |
puDArr, |
420 |
< |
electronArr, |
421 |
< |
&electronReferencePreSelection, |
422 |
< |
&electronReferenceIDMVASelection, |
423 |
< |
&electronIsoMVASelection, |
424 |
< |
muonArr, |
425 |
< |
&muonReferencePreSelection, |
426 |
< |
&muonIDPFSelection, |
427 |
< |
&muonIsoMVASelection); |
446 |
> |
// |
447 |
> |
// lepton/kinematic selection ... |
448 |
> |
// |
449 |
> |
failingLeptons.clear(); |
450 |
> |
passingLeptons.clear(); |
451 |
> |
EventData ret4l = apply_HZZ4L_reference_selection(ctrl, info, |
452 |
> |
vtxArr, |
453 |
> |
pfArr, |
454 |
> |
puDArr, |
455 |
> |
electronArr, |
456 |
> |
&electronReferencePreSelection, |
457 |
> |
&electronReferenceIDMVASelectionV1, // not used anymore! (Blame the damn regression implementation, not me!) |
458 |
> |
&electronReferenceIsoSelection, |
459 |
> |
muonArr, |
460 |
> |
&muonReferencePreSelection, |
461 |
> |
&muonIDPFSelection, |
462 |
> |
&muonReferenceIsoSelection); |
463 |
> |
|
464 |
> |
if( ctrl.debug ) cout << endl; |
465 |
|
|
466 |
+ |
if(ctrl.fakeScheme == "none") { // fill angletuple if event passes full selection |
467 |
+ |
if( ret4l.status.pass() ) { |
468 |
+ |
setNtupleVals(ctrl, ret4l, evtinfo, angles, kine, ji, weights, info, muonArr, electronArr, jetArr, metArr, simpleLeptonJets, fJetCorrector, puArr, puDArr, fJetIDMVA, vtxArr, triggerBits, hltTable, hltObjArr, fTrigObjs, &nt); |
469 |
+ |
nt.Fill(); |
470 |
+ |
|
471 |
+ |
int theZ1type = kine.l1type; |
472 |
+ |
int theZ2type = kine.l3type; |
473 |
|
|
474 |
+ |
if(ctrl.debug) cout << "simpleLeptonJets.size() = " << simpleLeptonJets.size() << endl; |
475 |
|
|
476 |
< |
if( ctrl.debug ) cout << endl; |
477 |
< |
|
478 |
< |
|
479 |
< |
int Z1type=0, Z2type=0; |
480 |
< |
if( ret4l.status.selectionBits.to_ulong() >= PASS_ZCANDIDATE ) { |
481 |
< |
if( abs(ret4l.Z1leptons[0].type) == 11 ) Z1type = 11; |
482 |
< |
if( abs(ret4l.Z1leptons[0].type) == 13 ) Z1type = 13; |
483 |
< |
} |
484 |
< |
if( ret4l.status.selectionBits.to_ulong() >= PASS_ZZCANDIDATE ) { |
485 |
< |
if( abs(ret4l.Z2leptons[0].type) == 11 ) Z2type = 11; |
486 |
< |
if( abs(ret4l.Z2leptons[0].type) == 13 ) Z2type = 13; |
487 |
< |
} |
488 |
< |
|
444 |
< |
cout << "bits: " << ret4l.status.selectionBits << "\t" |
445 |
< |
<< "Z1: " << Z1type << "\t" |
446 |
< |
<< "Z2: " << Z2type << endl; |
447 |
< |
cout << endl; |
448 |
< |
|
449 |
< |
|
450 |
< |
if( ret4l.status.pass() ) { |
451 |
< |
fillAngles(ret4l,angles); |
452 |
< |
fillKinematics(ret4l,kinematics); |
453 |
< |
fillEventInfo(info,evtinfo); |
454 |
< |
if( ctrl.mc) { |
455 |
< |
// fillGenInfo(ginfo,geninfo); |
456 |
< |
setEffiencyWeights(ret4l, weights); |
457 |
< |
if(ctrl.debug) |
458 |
< |
cout << "w: " << weights.w << "\t" |
459 |
< |
<< "won: " << weights.won << "\t" |
460 |
< |
<< "woff: " << weights.woff << endl; |
476 |
> |
if(kine.m4l > 100){ |
477 |
> |
if(simpleLeptonJets.size() >= 1) { |
478 |
> |
increment(counts,"m4l > 100 && >= 1 jets",theZ1type,theZ2type); |
479 |
> |
if(simpleLeptonJets.size() == 2) { |
480 |
> |
increment(counts,"m4l > 100 && 2 jets",theZ1type,theZ2type); |
481 |
> |
SimpleLepton j1 = simpleLeptonJets[0]; |
482 |
> |
SimpleLepton j2 = simpleLeptonJets[1]; |
483 |
> |
double deta = abs(j1.vec.Eta()-j2.vec.Eta()); |
484 |
> |
double mjj = (j1.vec + j2.vec).M(); |
485 |
> |
double fisher = 0.09407 * deta +0.00041581*mjj; |
486 |
> |
if(fisher > 0.4) increment(counts, "m4l > 100 && fisher > 0.4",theZ1type,theZ2type); |
487 |
> |
} |
488 |
> |
} |
489 |
|
} |
490 |
+ |
|
491 |
|
|
492 |
< |
/* |
493 |
< |
// only for Higgs < 140 |
494 |
< |
geninfo.weight *= h_w_hpt->GetBinContent(h_w_hpt->FindBin(geninfo.pt_zz)); |
495 |
< |
angles.bdt = reader->EvaluateMVA("BDTG"); |
496 |
< |
*/ |
497 |
< |
nt.Fill(); |
498 |
< |
|
499 |
< |
cerr << "PASS:: " |
500 |
< |
<< "\trun: " << evtinfo.run |
501 |
< |
<< "\tevt: " << evtinfo.evt |
502 |
< |
<< "\tlumi: " << evtinfo.lumi |
503 |
< |
<< "\tcostheta1: " << angles.costheta1 |
475 |
< |
<< "\tcostheta2: " << angles.costheta2 |
476 |
< |
<< "\tcostheta*: " << angles.costhetastar |
477 |
< |
<< "\tbdt: " << angles.bdt |
478 |
< |
<< endl; |
479 |
< |
pass++; |
480 |
< |
// if( pass > 3 ) break; |
481 |
< |
|
482 |
< |
} else if( ret4l.status.selectionBits.test(PASS_GOODZ1) ) { // save for fakes ... |
483 |
< |
// if(ctrl.debug) { |
484 |
< |
cout << "failingLeptons : " << failingLeptons.size() << endl; |
485 |
< |
for( int f=0; f<failingLeptons.size(); f++ ) { |
486 |
< |
SimpleLepton l = failingLeptons[f]; |
487 |
< |
cout << "f: " << f << "\t" |
488 |
< |
<< "type: " << l.type << "\t" |
489 |
< |
<< "pT: " << l.vec.Pt() << "\t" |
490 |
< |
<< "eta: " << l.vec.Eta() |
491 |
< |
<< endl; |
492 |
< |
// } |
492 |
> |
} else { |
493 |
> |
if(ctrl.debug) cout << "failing with code : " << ret4l.status.selectionBits << endl; |
494 |
> |
} |
495 |
> |
} else { // fill FOtuple and angletuple if event has a good Z and one or two extra leptons |
496 |
> |
assert(ctrl.fakeScheme=="ZPlusF" || ctrl.fakeScheme=="ZPlusFF"); |
497 |
> |
if(ret4l.status.selectionBits.test(PASS_GOODZ1) ) { |
498 |
> |
unsigned nleps = passingLeptons.size() + failingLeptons.size(); |
499 |
> |
if( (nleps>=3 && ctrl.fakeScheme=="ZPlusF" ) || (nleps>=4 && ctrl.fakeScheme=="ZPlusFF" ) ) { |
500 |
> |
addDummyZ2Lepotons(ret4l); |
501 |
> |
setNtupleVals(ctrl, ret4l, evtinfo, angles, kine, ji, weights, info, muonArr, electronArr, jetArr, metArr, simpleLeptonJets, fJetCorrector, puArr, puDArr, fJetIDMVA, vtxArr, triggerBits, hltTable, hltObjArr, fTrigObjs, &nt); |
502 |
> |
nt.Fill(); |
503 |
> |
foTree->Fill(); |
504 |
|
} |
505 |
< |
foTree.Fill(); |
505 |
> |
} else { |
506 |
> |
if(ctrl.debug) cout << "failing PASS_GOODZ1" << endl; |
507 |
|
} |
508 |
< |
} |
497 |
< |
|
498 |
< |
|
499 |
< |
foTree.getFile()->cd(); |
500 |
< |
foTree.getTree()->Write(); |
501 |
< |
nt.WriteClose(); |
502 |
< |
|
503 |
< |
cout << endl; |
504 |
< |
cout << endl; |
505 |
< |
for( int i=0; i<cutvec.size(); i++ ) { |
506 |
< |
cout << "cut: " << i << "\t" |
507 |
< |
<< "pass: " << cutvec[i]; |
508 |
< |
|
509 |
< |
if( i>PASS_TRIGGER&&i<=PASS_GOODZ1 ) |
510 |
< |
cout << "\t2e: " << zcutvec[0][i] |
511 |
< |
<< "\t2m: " << zcutvec[1][i]; |
512 |
< |
|
513 |
< |
if( i>PASS_ZCANDIDATE ) |
514 |
< |
cout << "\t4e: " << zzcutvec[0][i] |
515 |
< |
<< "\t4m: " << zzcutvec[1][i] |
516 |
< |
<< "\t2e2m: " << zzcutvec[2][i]; |
508 |
> |
} |
509 |
|
|
510 |
< |
cout << "\t" << cutstrs[i] << endl; |
510 |
> |
hltTable->Delete(); |
511 |
|
|
520 |
– |
cout << endl; |
512 |
|
} |
513 |
|
|
514 |
< |
} |
515 |
< |
|
516 |
< |
//---------------------------------------------------------------------------- |
526 |
< |
// |
527 |
< |
// Get primary vertices |
528 |
< |
// Assumes primary vertices are ordered by sum-pT^2 (as should be in CMSSW) |
529 |
< |
// NOTE: if no PV is found from fitting tracks, the beamspot is used |
530 |
< |
// |
531 |
< |
//---------------------------------------------------------------------------- |
532 |
< |
bool setPV(ControlFlags ctrl, |
533 |
< |
const mithep::Array<mithep::Vertex> * vtxArr, |
534 |
< |
mithep::Vertex & vtx) |
535 |
< |
//---------------------------------------------------------------------------- |
536 |
< |
{ |
537 |
< |
const Vertex *bestPV = 0; |
538 |
< |
float best_sumpt=-1; |
539 |
< |
|
540 |
< |
// good PV requirements |
541 |
< |
const UInt_t fMinNTracksFit = 0; |
542 |
< |
const Double_t fMinNdof = 4; |
543 |
< |
const Double_t fMaxAbsZ = 24; |
544 |
< |
const Double_t fMaxRho = 2; |
545 |
< |
|
546 |
< |
for(int i=0; i<vtxArr->GetEntries(); ++i) { |
547 |
< |
const Vertex *pv = (Vertex*)(vtxArr->At(i)); |
548 |
< |
if( ctrl.debug ) cout << "vertex :: " << i << "\tntrks: " << pv->NTracks() << endl; |
549 |
< |
|
550 |
< |
// Select best PV for corrected d0; if no PV passing cuts, the first PV in the collection will be used |
551 |
< |
// if(!pv->IsValid()) continue; |
552 |
< |
if(pv->NTracksFit() < fMinNTracksFit) continue; |
553 |
< |
if(pv->Ndof() < fMinNdof) continue; |
554 |
< |
if(fabs(pv->Z()) > fMaxAbsZ) continue; |
555 |
< |
if(pv->Position().Rho() > fMaxRho) continue; |
556 |
< |
|
557 |
< |
// take the first ... |
558 |
< |
bestPV = pv; |
559 |
< |
break; |
560 |
< |
|
561 |
< |
// this never reached ... |
562 |
< |
float tmp_sumpt=0; |
563 |
< |
for( int t=0; t<pv->NTracks(); t++ ) |
564 |
< |
tmp_sumpt += pv->Trk(t)->Pt(); |
565 |
< |
|
566 |
< |
if( tmp_sumpt > best_sumpt ) { |
567 |
< |
bestPV = pv; |
568 |
< |
best_sumpt = tmp_sumpt; |
569 |
< |
if( ctrl.debug) cout << "new PV set, pt : " << best_sumpt << endl; |
570 |
< |
} |
514 |
> |
if(ctrl.fakeScheme != "none") { |
515 |
> |
foTree->getFile()->cd(); |
516 |
> |
foTree->getTree()->Write(); |
517 |
|
} |
518 |
< |
|
573 |
< |
if(!bestPV) bestPV = vtxArr->At(0); |
574 |
< |
vtx.SetPosition(bestPV->X(),bestPV->Y(),bestPV->Z()); |
575 |
< |
vtx.SetErrors(bestPV->XErr(),bestPV->YErr(),bestPV->ZErr()); |
576 |
< |
return true; |
577 |
< |
}; |
578 |
< |
|
579 |
< |
|
580 |
< |
|
581 |
< |
//---------------------------------------------------------------------------- |
582 |
< |
void setEffiencyWeights(EventData & evtdat, WeightStruct & weights ) |
583 |
< |
//---------------------------------------------------------------------------- |
584 |
< |
{ |
585 |
< |
vector<SimpleLepton> lepvec = evtdat.Z1leptons; |
586 |
< |
lepvec.insert(lepvec.end(), evtdat.Z2leptons.begin(), evtdat.Z2leptons.end()); |
587 |
< |
double w_offline=-1, werr_offline=0; |
588 |
< |
double w_online=-1, werr_online=0; |
589 |
< |
|
590 |
< |
vector< pair <double,double> > wlegs; // pair here is eff & err |
591 |
< |
vector< pair <float,float> > mvec; // pair here is eta & pt |
518 |
> |
nt.WriteClose(); |
519 |
|
|
520 |
< |
// vector< pair <float,float> > evec; // pair here is eta & pt |
521 |
< |
// now deal with medium vs loose |
522 |
< |
vector< pair< bool, pair <float,float> > > evec; // pair here is eta & pt |
523 |
< |
|
597 |
< |
for( int k=0; k<lepvec.size(); k++ ) { |
598 |
< |
if( abs(lepvec[k].type) == 13 ) { |
599 |
< |
mvec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) ); |
600 |
< |
wlegs.push_back( muonPerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()), |
601 |
< |
lepvec[k].vec.Pt() ) ); |
602 |
< |
} else { |
603 |
< |
|
604 |
< |
// now deal with medium vs loose |
605 |
< |
// evec.push_back( std::pair<float,float> (fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()) ); |
606 |
< |
|
607 |
< |
std::pair<float,float> tmppair(fabs(lepvec[k].vec.Eta()), lepvec[k].vec.Pt()); |
608 |
< |
evec.push_back( std::pair<bool, std::pair<float,float> > (lepvec[k].isTight, tmppair) ); |
609 |
< |
|
610 |
< |
// wlegs.push_back( elePerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()), |
611 |
< |
// lepvec[k].vec.Pt() ) ); |
612 |
< |
|
613 |
< |
wlegs.push_back( elePerLegOfflineEfficiencyWeight( fabs(lepvec[k].vec.Eta()), |
614 |
< |
lepvec[k].vec.Pt(), |
615 |
< |
lepvec[k].isTight ) ); |
616 |
< |
|
617 |
< |
} |
520 |
> |
for(unsigned icut=0; icut<cuts.size(); icut++) { |
521 |
> |
map<TString,int> thisCount(*(counts[cuts[icut]])); |
522 |
> |
cerr << setw(25) << cuts[icut] << "\t" << setw(9) << thisCount["all"] << " (" |
523 |
> |
<< setw(9) << thisCount["4e"] << setw(9) << thisCount["4m"] << setw(9) << thisCount["2e2m"] << ")" << endl; |
524 |
|
} |
619 |
– |
|
620 |
– |
pair<double,double> offpair = getOfflineEfficiencyWeight( wlegs ); |
621 |
– |
weights.woff = offpair.first; |
622 |
– |
weights.werroff = offpair.second; |
623 |
– |
|
624 |
– |
pair<double,double> onpair = getOnlineEfficiencyWeight( mvec, evec ); |
625 |
– |
weights.won = onpair.first; |
626 |
– |
weights.werron = onpair.second; |
627 |
– |
} |
525 |
|
|
526 |
+ |
tmp_file.Close(); |
527 |
|
|
630 |
– |
//---------------------------------------------------------------------------- |
631 |
– |
void initRunLumiRangeMap() |
632 |
– |
//---------------------------------------------------------------------------- |
633 |
– |
{ |
634 |
– |
/* |
635 |
– |
rlrm.AddJSONFile(std::string("./data/Cert_136033-149442_7TeV_Apr21ReReco_Collisions10_JSON.txt")); |
636 |
– |
// rlrm.AddJSONFile(std::string("./data/Cert_160404-173244_7TeV_PromptReco_Collisions11_JSON_v2.txt")); |
637 |
– |
rlrm.AddJSONFile(std::string("./data/Cert_160404-178078_7TeV_PromptReco_Collisions11_JSON.txt")); |
638 |
– |
rlrm.AddJSONFile(std::string("./data/Cert_160404-163869_7TeV_May10ReReco_Collisions11_JSON_v3.txt")); |
639 |
– |
rlrm.AddJSONFile(std::string("./data/Cert_170249-172619_7TeV_ReReco5Aug_Collisions11_JSON.txt")); |
640 |
– |
// r11b |
641 |
– |
rlrm.AddJSONFile(std::string("./data/Cert_160404-180252_7TeV_PromptReco_Collisions11_JSON.txt")); |
642 |
– |
*/ |
643 |
– |
}; |
644 |
– |
|
645 |
– |
|
646 |
– |
//---------------------------------------------------------------------------- |
647 |
– |
void initPUWeights() |
648 |
– |
//---------------------------------------------------------------------------- |
649 |
– |
{ |
650 |
– |
TFile * puf = new TFile("data/PileupReweighting.Summer11DYmm_To_Full2011.root"); |
651 |
– |
hpu = (TH1D*)(puf->Get("puWeights")); |
528 |
|
} |
529 |
< |
|
530 |
< |
//---------------------------------------------------------------------------- |
531 |
< |
double getPUWeight(unsigned npu) |
656 |
< |
//---------------------------------------------------------------------------- |
529 |
> |
//---------------------------------------------------------------------------------------- |
530 |
> |
void addDummyZ2Lepotons(EventData &ret4l) |
531 |
> |
// for fake ntuples, we don't know which are the Z2 leptons yet, but need to fill something so we don't seg fault |
532 |
|
{ |
533 |
< |
return hpu->GetBinContent(hpu->FindBin(npu)); |
533 |
> |
SimpleLepton dum1,dum2; |
534 |
> |
dum1.vec.SetPtEtaPhiM(1,2,3,4); |
535 |
> |
dum2.vec.SetPtEtaPhiM(5,6,7,8); |
536 |
> |
ret4l.Z2leptons.push_back(dum1); |
537 |
> |
ret4l.Z2leptons.push_back(dum2); |
538 |
|
} |
539 |
< |
|
540 |
< |
//---------------------------------------------------------------------------- |
541 |
< |
void initEvtRhoMap( map<unsigned, float> & evtrhoMap ) |
542 |
< |
//---------------------------------------------------------------------------- |
539 |
> |
//---------------------------------------------------------------------------------------- |
540 |
> |
void setNtupleVals(ControlFlags &ctrl, EventData &ret4l, InfoStruct &evtinfo, Angles &angles, KinematicsStruct &kine, |
541 |
> |
JetInfoStruct &ji, WeightStruct &weights, EventHeader *info, Array<Muon> *muonArr, Array<Electron> *electronArr, Array<PFJet> *jetArr, Array<PFMet> *metArr, |
542 |
> |
vector<SimpleLepton> &simpleLeptonJets, FactorizedJetCorrector *fJetCorrector, Array<PileupInfo> *puArr, |
543 |
> |
Array<PileupEnergyDensity> *puDArr, JetIDMVA *fJetIDMVA, Array<Vertex> *vtxArr, bitset<1024> &triggerBits, |
544 |
> |
TriggerTable *hltTable, Array<TriggerObject> *hltObjArr, TriggerObjectsTable *fTrigObjs, AngleTuple *nt) |
545 |
|
{ |
546 |
< |
unsigned evt; |
666 |
< |
float rho; |
546 |
> |
fillKinematics(ret4l,kine); |
547 |
|
|
548 |
< |
FILE * f = fopen("evtrho.dat", "r"); |
549 |
< |
while( fscanf( f, "%u:%f", &evt, &rho ) ) { |
550 |
< |
evtrhoMap[evt] = rho; |
548 |
> |
//the last two arguments are usePt and useY respectively |
549 |
> |
fillMela(kine,0,0); |
550 |
> |
|
551 |
> |
fillMassErrors(ctrl,ret4l,muonArr,electronArr,kine); |
552 |
> |
|
553 |
> |
TLorentzVector pfmet; pfmet.SetPxPyPzE(metArr->At(0)->Mex(),metArr->At(0)->Mey(),0,0); |
554 |
> |
fillEventInfo( info, pfmet, evtinfo, ctrl.mc ? getNPU(puArr) : 0, vtxArr->GetEntries()); |
555 |
> |
fillJets(ctrl, jetArr, simpleLeptonJets, fJetCorrector, puDArr, fJetIDMVA, vtxArr, ret4l); |
556 |
> |
fillJetInfo(simpleLeptonJets, ji, ctrl); |
557 |
> |
//fillJetAngleBranches(ji, ret4l, nt, ctrl); |
558 |
> |
evtinfo.status = 0; |
559 |
> |
evtinfo.status = setHltBits(ti, triggerBits); |
560 |
> |
|
561 |
> |
if(ctrl.fakeScheme == "none") { |
562 |
> |
fillAngles(ret4l,angles); |
563 |
> |
|
564 |
> |
if( ctrl.mc) { |
565 |
> |
if(ctrl.use_eff_weights) setEfficiencyWeights(ctrl, ctrl.era, ret4l, triggerBits, hltTable, hltObjArr, fTrigObjs, weights); |
566 |
> |
if(ctrl.debug) |
567 |
> |
cout << "w: " << weights.w << "\t" << "npuw: " << weights.npuw << "\t" << "won: " << weights.won << "\t" << "woff: " << weights.woff << endl; |
568 |
> |
} |
569 |
> |
|
570 |
> |
if(ctrl.debug) { |
571 |
> |
stringstream ss; |
572 |
> |
ss << "PASS:: " << "\trun: " << evtinfo.run << "\tevt: " << evtinfo.evt << "\tlumi: " << evtinfo.lumi << "\tchannel: " |
573 |
> |
<< kine.channel << "\tm4l: " << kine.m4l << "\tmZ1: " << kine.mZ1 << "\tmZ2: " << kine.mZ2 << endl; |
574 |
> |
cerr << ss.str(); |
575 |
> |
cout << ss.str(); |
576 |
> |
} |
577 |
> |
} else { |
578 |
> |
// fillAngles(ret4l,angles); // at this point we don't know which four leptons will constitute the ZZ system, so doesn't make sense to fill the angles |
579 |
|
} |
672 |
– |
|
580 |
|
} |