ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/mstein/triggerStudy/triggerStudy.C
Revision: 1.6
Committed: Mon Dec 17 10:44:45 2012 UTC (12 years, 4 months ago) by mstein
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.5: +61 -29 lines
Log Message:
updated paths

File Contents

# Content
1 /////////////////////////////////////////////////////////
2 // Trigger Study
3 /////////////////////////////////////////////////////////
4 //to run, execute: root run.C
5 //
6 //In this file, the implementation of the T&P method and
7 //the evaluation of the trigger efficiency is performed.
8 //some documentation can be found here:
9 //https://twiki.cern.ch/twiki/bin/viewauth/CMSsandbox/MatthiasSteinTriggerStudy
10
11
12 #include <iostream>
13 #include <fstream>
14 #include <map>
15 #include <vector>
16
17 #include <Math/LorentzVector.h>
18
19 #include "TString.h"
20 #include "TH1.h"
21 #include "TH1D.h"
22 #include "TTree.h"
23 #include "TFile.h"
24
25 #include "plotMaker.h"
26 #include "tools.h"
27 //--------------------------------------------------- declare some special types
28 #ifdef __MAKECINT__
29 #pragma link C++ class pair<string,bool>+;
30 #pragma link C++ class pair<string,string>+;
31 #pragma link C++ class map<string,bool>+;
32 #pragma link C++ class map<string,string>+;
33 #pragma link C++ class ROOT::Math::PtEtaPhiM4D<float>+;
34 #pragma link C++ class ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> >+;
35 #pragma link C++ class vector<ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> > >+;
36 #endif
37
38 typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiM4D<float> > LorentzM;
39
40 //--------------------------------------------------- Structure to store the whole event information in order to pass it to functions
41 struct eventInfo{
42 vector<string> *DESYtriggerElMatchedTriggerFilter;
43 vector<string> *DESYtriggerMuMatchedTriggerFilter;
44 map<string, string> *DESYtriggerNameMap;
45 vector<int> *electronCharge;
46 vector<LorentzM> *electronP4;
47 int Event;
48 double HT;
49 vector<LorentzM> *jetP4;
50 double MET;
51 vector<int> *muonCharge;
52 vector<LorentzM> *muonP4;
53 map<string, int> *prescaled;
54 int Run;
55 int PUInter;
56 double PUWeight;
57 map<string, bool> *triggered;
58 vector<vector<TString> > triggerFilterEl;
59 vector<vector<TString> > triggerFilterMu;
60 double Weight;
61 };
62
63 //--------------------------------------------------- functions
64 // void init(TString inputFileName);
65 void handleTriggerMatchingInfo(eventInfo& evt);
66 void initBranches(TTree *tree, eventInfo &event);
67 void investigateTriggerFilterLep(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName, TString label);
68 void Loop(TString inputFileName, TString outputFileName);
69 void makeEfficiencyPlots(map<TString, plotMaker> &plot);
70 void plotHistograms(eventInfo& evt, map<TString, plotMaker> &plot);
71 void printElectronInfo(eventInfo& evt, ostream& os=cout);
72 void printEventInfo(eventInfo& evt, ostream& os=cout);
73 void printMuonInfo(eventInfo& evt, ostream& os=cout);
74 void setEfficiencyErrors(TH1D* &eff, TH1D* &h_passed, TH1D* &h_failed);
75
76 //--------------------------------------------------- global Variables
77 vector<string> triggerFilterToCheckMu;
78 vector<string> triggerFilterToCheckEl;
79 //----------------------------------------
80 const Int_t binPt = 25;
81 const Double_t lowPt = 0.;
82 const Double_t upPt = 250.;
83 //----------------------------------------
84 const Int_t binEta = 50;
85 const Double_t lowEta = -2.5;
86 const Double_t upEta = 2.5;
87 //----------------------------------------
88 const Int_t binPU = 64;
89 const Double_t lowPU = 0.;
90 const Double_t upPU = 64.;
91 //----------------------------------------
92 const Int_t binMinv = 40;
93 const Double_t lowMinv = 50.;
94 const Double_t upMinv = 130.;
95
96
97 //================================================================================= triggerStudy
98 //main function which will be called
99 //here, it is possible to run over multiple files and to estimate trigger efficiencies for different samples
100 void triggerStudy(){
101
102 //--------------------------------------------------- define the muon TriggerFilters to be checked
103 triggerFilterToCheckMu.push_back("hltL1sMu16Eta2p1");
104 triggerFilterToCheckMu.push_back("hltL2fL1sMu16Eta2p1L1f0L2Filtered16Q");
105 triggerFilterToCheckMu.push_back("hltL3crIsoL1sMu16Eta2p1L1f0L2f16QL3f20L3crIsoFiltered10");
106 triggerFilterToCheckMu.push_back("hltIsoMu202p1TriCentralPFJet30MuCleaned");
107 triggerFilterToCheckMu.push_back("hltIsoMu202p1TriCentralPFJet303020MuCleaned");
108 triggerFilterToCheckMu.push_back("hltDiMuonGlb22Trk8DzFiltered0p2");
109 triggerFilterToCheckMu.push_back("hltL3crIsoL1sMu14Eta2p1L1f0L2f16QL3f20L3crIsoRhoFiltered0p15");
110 //--------------------------------------------------- define the electron TriggerFilters to be checked
111 // triggerFilterToCheckEl.push_back("hltL1sL1SingleEG20orL1SingleEG22");
112 // triggerFilterToCheckEl.push_back("hltEle25CaloIdVTTrkIdTCaloIsoTTrkIsoTTrackIsoFilter");
113 // triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFJet30EleCleaned");
114 // triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFNoPUJet30EleCleaned");
115 // triggerFilterToCheckEl.push_back("hltEle25CaloIdVTCaloIsoTTrkIdTTrkIsoTTriCentralPFNoPUJet303020EleCleaned");
116
117 //--------------------------------------------------- define the sample(s) to loop over and the output file
118 //these are small test files of ONE single file af a sample
119 // Loop("/scratch/hh/current/cms/user/mstein/output/TTJets_mu_tree.root", "out_TTJets_mu_tree.root");
120 // Loop("/scratch/hh/current/cms/user/mstein/output/TTJets_el_tree.root", "out_TTJets_el_tree.root");
121 // Loop("/scratch/hh/current/cms/user/mstein/output/DYJets_mu_tree.root", "out_DYJets_mu_tree.root");
122 // Loop("/scratch/hh/current/cms/user/mstein/output/DYJets_el_tree.root", "out_DYJets_el_tree.root");
123 // Loop("/scratch/hh/current/cms/user/mstein/output/SingleMu_mu_tree.root", "out_SingleMu_mu_tree.root");
124 // Loop("/scratch/hh/current/cms/user/mstein/output/SingleElectron_el_tree.root", "out_SingleElectron_el_tree.root");
125 //these are the whole samples
126 // Loop("/scratch/hh/current/cms/user/mstein/output/TTJets/SUMMER12/TTJets_SUMMER12_TrigStudy-mu_noTail_Tree.root", "out_TTJets_SUMMER12_TrigStudy-mu_noTail_Tree.root");
127 // Loop("/scratch/hh/current/cms/user/mstein/output/TTJets/SUMMER12/TTJets_SUMMER12_TrigStudy-el_noTail_Tree.root", "out_TTJets_SUMMER12_TrigStudy-el_noTail_Tree.root");
128 Loop("/scratch/hh/current/cms/user/mstein/output/DYJets/M-50/DYJets_M-50_TrigStudy-mu_noTail_Tree.root", "out_DYJets_M-50_TrigStudy-mu_noTail_Tree.root");
129 // Loop("/scratch/hh/current/cms/user/mstein/output/DYJets/M-50/DYJets_M-50_TrigStudy-el_noTail_Tree.root", "out_DYJets_M-50_TrigStudy-el_noTail_Tree.root");
130 // Loop("/scratch/hh/current/cms/user/mstein/output/SingleElectron/Run2012B-13Jul2012-v1/SingleElectron_Run2012B-13Jul2012-v1_TrigStudy-el_noTail_Tree.root", "out_SingleElectron_Run2012B-13Jul2012-v1_TrigStudy-el_noTail_Tree.root");
131 // Loop("/scratch/hh/current/cms/user/mstein/output/SingleMu/Run2012B-13Jul2012-v1/SingleMu_Run2012B-13Jul2012-v1_TrigStudy-mu_noTail_Tree.root", "out_SingleMu_Run2012B-13Jul2012-v1_TrigStudy-mu_noTail_Tree.root");
132
133
134 }
135
136
137 //================================================================================= plotHistograms
138 void plotHistograms(eventInfo& evt, map<TString, plotMaker> &plot){
139 plot["event"].addh1d("muMultiplicity", "muMultiplicity", "muon multiplicity", "events", 5, -0.5, 4.5, evt.muonP4->size(), evt.Weight);
140 plot["event"].addh1d("elMultiplicity", "elMultiplicity", "electron multiplicity", "events", 5, -0.5, 4.5, evt.electronP4->size(), evt.Weight);
141 plot["event"].addh1d("jetMultiplicity", "jetMultiplicity", "jet multiplicity", "events", 10, -0.5, 9.5, evt.jetP4->size(), evt.Weight);
142 if(evt.electronP4->size()==2){
143 Double_t minv = (evt.electronP4->at(0) + evt.electronP4->at(1)).mass();
144 plot["event"].addh1d("Minv", "Minv", "M_{inv} [GeV/c^{2}]", "events", binMinv, lowMinv, upMinv, minv, evt.Weight);
145 }
146
147
148 for(Int_t i=0,N=triggerFilterToCheckMu.size(); i<N; ++i){
149 investigateTriggerFilterLep(evt, plot, triggerFilterToCheckMu.at(i), "muon");
150 }
151 for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
152 investigateTriggerFilterLep(evt, plot, triggerFilterToCheckEl.at(i), "electron");
153 }
154 }
155
156
157 //================================================================================= investigateTriggerFilterMu
158 //check, if the positive lepton can be matched to the trigger filter
159 //if yes, take this lepton as Tag
160 //then probe the other lepton
161 void investigateTriggerFilterLep(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName, TString label){
162 vector<LorentzM>* lepP4;
163 vector<int>* lepCharge;
164 vector<vector<TString> > triggerFilterLep;
165 //--------------------------------------------------- chose input variables according to label
166 if(label = "muon"){
167 lepP4 = (evt.muonP4);
168 lepCharge = (evt.muonCharge);
169 triggerFilterLep = (evt.triggerFilterMu);
170 }
171 else if(label = "electron"){
172 lepP4 = (evt.electronP4);
173 lepCharge = (evt.electronCharge);
174 triggerFilterLep = (evt.triggerFilterEl);
175 }
176 else{
177 cout<<"WARNING: Label for 'investigateTriggerFilter' is neither 'muon' nor 'electron'!"<<endl;
178 cout<<"label = " << label << endl;
179 return;
180 }
181 //--------------------------------------------------- do some consistency check - just to be sure
182 if(lepP4->size() != 2){
183 cout<<"WARNING: The size of the " << label << "vector is NOT two!"<<endl;
184 cout<<"lepP4->size() = " << lepP4->size() << endl;
185 return;
186 }
187 if(lepP4->size() != triggerFilterLep.size()){
188 cout<<"WARNING: The size of the "<< label << " vector and Matching-info vector is different!"<<endl;
189 cout<<"lepP4->size() = " << lepP4->size() << endl;
190 cout<<"triggerFilterLep.size() = " << triggerFilterLep.size() << endl;
191 return;
192 }
193 if(lepP4->size() != lepCharge->size()){
194 cout<<"WARNING: The size of the "<< label << "muon vector and muon-charge vector is different!"<<endl;
195 cout<<"lepP4->size() = " << lepP4->size() << endl;
196 cout<<"lepCharge->size() = " << lepCharge->size() << endl;
197 return;
198 }
199 // if(evt.electronP4->size() != evt.triggerFilterEl.size()){
200 // cout<<"WARNING: The size of the electron vector and Matching-info vector is different!"
201 // cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
202 // cout<<"evt.triggerFilterEl->size() = " << evt.triggerFilterEl->size() << endl;
203 // return;
204 // }
205 // if(evt.electronP4->size() != evt.electronP4.size()){
206 // cout<<"WARNING: The size of the electron vector and muon-charge vector is different!"
207 // cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
208 // cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
209 // return;
210 // }
211 //--------------------------------------------------- check which lepton has positive and which negative charge
212 Int_t index_pos = 0;
213 Int_t index_neg = 1;
214 if(lepCharge->at(0)==-1){
215 index_pos = 0;
216 index_pos = 1;
217 }
218 //--------------------------------------------------- check which lepton has been matched to the trigger filter
219 vector<bool> matched;
220 for(Int_t i=0,N=lepP4->size(); i<N; ++i){
221 matched.push_back(false);
222 for(Int_t j=0,M=triggerFilterLep.at(i).size(); j<M; ++j){
223 if(filterName == triggerFilterLep.at(i).at(j)){
224 matched.back() = true;
225 break;
226 }
227 }
228 }
229 //--------------------------------------------------- do the plotting according to the outcome, which lepton could be matched
230
231 //if the positive lepton IS NOT matched => no Tag found => return
232 if(!matched.at(index_pos)) return;
233
234 //if the positive lepton IS matched => plot Tag and Probe, and passed and failed
235 //plot the Tag
236 plot[filterName].addh1d(filterName+"_pt_tag", filterName+"_pt_tag", "p_{T}(tag) [GeV/c]", "events", binPt, lowPt, upPt, lepP4->at(index_pos).Pt(), evt.Weight);
237 plot[filterName].addh1d(filterName+"_eta_tag", filterName+"_eta_tag", "#eta(tag)", "events", binEta, lowEta, upEta, lepP4->at(index_pos).eta(), evt.Weight);
238 //plot the Probe
239 plot[filterName].addh1d(filterName+"_pt_probe", filterName+"_pt_probe", "p_{T}(probe) [GeV/c]", "events", binPt, lowPt, upPt, lepP4->at(index_neg).pt(), evt.Weight);
240 plot[filterName].addh1d(filterName+"_eta_probe", filterName+"_eta_probe", "#eta(probe)", "events", binEta, lowEta, upEta, lepP4->at(index_neg).eta(), evt.Weight);
241 //plot the passed
242 if(matched.at(index_neg)){
243 Double_t minv = (lepP4->at(index_pos) + lepP4->at(index_neg)).mass();
244 plot[filterName].addh1d(filterName+"_pt_passed", filterName+"_pt_passed", "p_{T}(passed) [GeV/c]", "events", binPt, lowPt, upPt, lepP4->at(index_neg).pt(), evt.Weight);
245 plot[filterName].addh1d(filterName+"_eta_passed", filterName+"_eta_passed", "#eta(passed)", "events", binEta, lowEta, upEta, lepP4->at(index_neg).eta(), evt.Weight);
246 plot[filterName].addh1d(filterName+"_PUVertices_passed", filterName+"_PUVertices_passed", "number of primary vertices (passed)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
247 plot[filterName].addh1d(filterName+"_Minv_passed", filterName+"_Minv_passed", "M_{inv} [GeV/c^{2}]", "events", binMinv, lowMinv, upMinv, minv, evt.Weight);
248 }
249 else{//plot the failed
250 Double_t minv = (lepP4->at(index_pos) + lepP4->at(index_neg)).mass();
251 plot[filterName].addh1d(filterName+"_pt_failed", filterName+"_pt_failed", "p_{T}(failed) [GeV/c]", "events", binPt, lowPt, upPt, lepP4->at(index_neg).pt(), evt.Weight);
252 plot[filterName].addh1d(filterName+"_eta_failed", filterName+"_eta_failed", "#eta(failed)", "events", binEta, lowEta, upEta, lepP4->at(index_neg).eta(), evt.Weight);
253 plot[filterName].addh1d(filterName+"_PUVertices_failed", filterName+"_PUVertices_failed", "number of primary vertices (failed)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
254 plot[filterName].addh1d(filterName+"_Minv_failed", filterName+"_Minv_failed", "M_{inv} [GeV/c^{2}]", "events", binMinv, lowMinv, upMinv, minv, evt.Weight);
255 }
256 }
257
258
259 //================================================================================= makeEfficiencyPlots
260 //eff = passed/(passed+failed)
261 void makeEfficiencyPlots(map<TString, plotMaker> &plot){
262 cout<<"----------> Create efficiency histograms"<<endl;
263
264 vector<string> allTriggerFiltersToCheck = triggerFilterToCheckMu;
265 for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
266 allTriggerFiltersToCheck.push_back(triggerFilterToCheckEl.at(i));
267 }
268 vector<string> histNames;
269 histNames.push_back("_pt_passed");
270 histNames.push_back("_eta_passed");
271 histNames.push_back("_PUVertices_passed");
272
273 for(Int_t i=0,N=allTriggerFiltersToCheck.size(); i<N; ++i){
274 TString triggerFilter = allTriggerFiltersToCheck.at(i);
275 for(Int_t j=0,M=histNames.size(); j<M; ++j){
276 TString passedName = triggerFilter;
277 TString failedName = triggerFilter;
278 TString effName = triggerFilter;
279 passedName += histNames.at(j);
280 failedName += histNames.at(j);
281 effName += histNames.at(j);
282 failedName.ReplaceAll("passed", "failed");
283 effName.ReplaceAll("passed", "eff");
284
285 if(plot[triggerFilter].isNewh1d(passedName)) continue;
286 if(plot[triggerFilter].isNewh1d(failedName)) continue;
287
288 TH1D *h_passed = (TH1D*)(plot[triggerFilter].h1d[passedName])->Clone();
289 TH1D *h_failed = (TH1D*)(plot[triggerFilter].h1d[failedName])->Clone();
290 TH1D *den = (TH1D*)(h_passed)->Clone();
291 den->Add((TH1D*)(h_failed)->Clone());
292
293 TH1D *eff = (TH1D*)(h_passed)->Clone();
294 eff->Divide(den);
295 setEfficiencyErrors(eff, h_passed, h_failed);
296
297 eff->SetName(effName);
298 eff->SetTitle(effName);
299
300 TString xTitle_eff = eff->GetXaxis()->GetTitle();
301 xTitle_eff.ReplaceAll("(passed)","");
302 eff->GetXaxis()->SetTitle(xTitle_eff);
303 eff->GetYaxis()->SetTitle("efficiency");
304 //eff->SetAxisRange(0.,1.,"Y");
305 plot[triggerFilter].addh1d(eff);
306 }
307 }
308 }
309
310
311 //================================================================================= setEfficiencyErrors
312 //p = passed; f = failed
313 // f(eff) = p/(p+f)
314 // df/dt = 1/(p+f) - p/(p+f)^2
315 // df/dp = - p/(p+f)^2
316 // eff_sig = sqrt(p_sig^2*(df/dp)^2 + f_sig^2*(df/df)^2)
317 void setEfficiencyErrors(TH1D* &eff, TH1D* &h_passed, TH1D* &h_failed){
318 if(h_passed->GetNbinsX() != h_failed->GetNbinsX()){
319 cout<<"WARNING: the number of bins of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
320 cout<<"Number of bins (passed | failed): " << h_passed->GetNbinsX() << " | " << h_failed->GetNbinsX() << endl;
321 }
322 if(h_passed->GetXaxis()->GetXmax() != h_failed->GetXaxis()->GetXmax()){
323 cout<<"WARNING: the upper x-range of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
324 cout<<"Number of bins (passed | failed): " << h_passed->GetXaxis()->GetXmax() << " | " << h_failed->GetXaxis()->GetXmax() << endl;
325 }
326 if(h_passed->GetXaxis()->GetXmin() != h_failed->GetXaxis()->GetXmin()){
327 cout<<"WARNING: the lower x-range of the passed and failed histograms differ! The efficiency calculation might be depricated!"<<endl;
328 cout<<"Number of bins (passed | failed): " << h_passed->GetXaxis()->GetXmin() << " | " << h_failed->GetXaxis()->GetXmin() << endl;
329 }
330 //loop over single bins of histograms an set the efficiency error
331 for(Int_t i=1,N=h_passed->GetNbinsX(); i<=N; ++i){
332 Double_t p = h_passed ->GetBinContent(i);
333 Double_t f = h_failed->GetBinContent(i);
334 Double_t p_sig = h_passed ->GetBinError(i);
335 Double_t f_sig = h_failed->GetBinError(i);
336
337 if(eff->GetBinContent(i)==0.) continue;
338
339 Double_t dedt = 1/(p+f)-p/((p+f)*(p+f));
340 Double_t dedp = -p/((p+f)*(p+f));
341
342 Double_t eff_sig = sqrt(p_sig*p_sig*dedt*dedt + f_sig*f_sig*dedp*dedp);
343 eff->SetBinError(i, eff_sig);
344 }
345 }
346
347
348 //================================================================================= Loop
349 void Loop(TString inputFileName, TString outputFileName){
350 //timer.Start("Loop");
351 TString treeName = "trigStudyTree";
352 map<TString, plotMaker> myPlots;
353
354 //--------------------------------------------------- open file and read in leafs
355 TFile* inputFile = new TFile(inputFileName,"READ");
356 if (!inputFile->IsOpen()) std::cout<<"Could not open file "<<inputFileName<<std::endl;
357 TTree* inputTree= (TTree*)inputFile->Get(treeName);
358
359 eventInfo event;
360 initBranches(inputTree, event);
361 //--------------------------------------------------- Event Loop
362 int nEvents=inputTree->GetEntries();
363 cout<<endl;
364 cout<<"///////////////////////////////////////////////////////////////////////////////"<<endl;
365 cout<<"////////////////////////////// Starting Event Loop ////////////////////////////"<<endl;
366 cout<<"///////////////////////////////////////////////////////////////////////////////"<<endl;
367 cout<<"Investigate File: " << inputFileName << endl;
368 cout<<"Read out TTree: " << treeName << endl;
369 cout<<"Number of events: " << nEvents << endl;
370 cout<<"-------------------------------------------------------------------------------"<<endl;
371 cout<<endl;
372 cout<<"Looping over events..."<<endl;
373 nEvents=50000;
374 for (int iEvent=0; iEvent<nEvents; ++iEvent){
375 tools::progress(iEvent);
376 inputTree->GetEntry(iEvent);
377 handleTriggerMatchingInfo(event);
378 // printEventInfo(event);
379 // printMuonInfo(event);
380 //printElectronInfo(event);
381 plotHistograms(event, myPlots);
382 //tools::progressReset();
383 }
384
385 makeEfficiencyPlots(myPlots);
386 //--------------------------------------------------- Save all histograms in the appropriate folders
387 TString option="RECREATE";
388 for(typename map<TString, plotMaker>::iterator it = myPlots.begin(); it != myPlots.end(); ++it){
389 it->second.savePlots(outputFileName, option, it->first);
390 option="UPDATE";
391 }
392 //timer.Stop("Loop");
393 }
394
395
396 //================================================================================= handleTriggerMatchingInfo
397 // converting the trigger Info to a more practical format: vector<string> --> vector<vector<string> >
398 void handleTriggerMatchingInfo(eventInfo& evt){
399 evt.triggerFilterEl.clear();
400 evt.triggerFilterMu.clear();
401 //make simple check
402
403 if(evt.electronP4->size() != evt.DESYtriggerElMatchedTriggerFilter->size()){
404 cout<<"WARNING: Size of electron vector does not correspond to size of matching info vector!"<<endl;
405 cout<<"evt.electronP4->size() = " << evt.electronP4->size() << endl;
406 cout<<"evt.DESYtriggerElMatchedTriggerFilter->size() = " << evt.DESYtriggerElMatchedTriggerFilter->size() << endl;
407
408 tools::printVector("evt.electronP4", *(evt.electronP4));
409 tools::printVector("evt.DESYtriggerElMatchedTriggerFilter", *(evt.DESYtriggerElMatchedTriggerFilter));
410
411 return;
412 }
413 if(evt.muonP4->size() != evt.DESYtriggerMuMatchedTriggerFilter->size()){
414 cout<<"WARNING: Size of muon vector does not correspond to size of matching info vector!"<<endl;
415 cout<<"evt.muonP4->size() = " << evt.muonP4->size() << endl;
416 cout<<"evt.DESYtriggerMuMatchedTriggerFilter->size() = " << evt.DESYtriggerMuMatchedTriggerFilter->size() << endl;
417 return;
418 }
419
420 for(Int_t i=0,N=evt.electronP4->size(); i<N; ++i){
421 evt.triggerFilterEl.push_back(tools::splitToWords(evt.DESYtriggerElMatchedTriggerFilter->at(i), ";"));
422 //cout<<"evt.DESYtriggerElMatchedTriggerFilter->at("<<i<<") = " << evt.DESYtriggerElMatchedTriggerFilter->at(i) << endl;
423 //tools::printVector<TString>("triggerFilterEl", evt.triggerFilterEl.at(i));
424 }
425 // printMuonInfo(evt);
426 // printElectronInfo(evt);
427
428 //tools::printVector("evt.triggerFilterMu", evt.triggerFilterMu);
429 // tools::printVector("evt.DESYtriggerMuMatchedTriggerFilter", *(evt.DESYtriggerMuMatchedTriggerFilter));
430 // tools::printVector("evt.muonP4", *(evt.muonP4));
431
432 for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
433 evt.triggerFilterMu.push_back(tools::splitToWords(evt.DESYtriggerMuMatchedTriggerFilter->at(i), ";"));
434 //cout<<"evt.DESYtriggerMuMatchedTriggerFilter->at("<<i<<") = " << evt.DESYtriggerMuMatchedTriggerFilter->at(i) << endl;
435 //tools::printVector<TString>("triggerFilterMu", evt.triggerFilterMu.at(i));
436 }
437 //cout<<"size: muonP4 | filterinfo = "
438 //<<evt.muonP4->size() << " | "
439 //<<evt.DESYtriggerMuMatchedTriggerFilter->size()
440 //<<endl;
441 }
442
443
444 //================================================================================= initBranches
445 void initBranches(TTree *tree, eventInfo& evt){
446 //--------------------------------------------------- initialize event info
447 evt.DESYtriggerElMatchedTriggerFilter = 0;
448 evt.DESYtriggerMuMatchedTriggerFilter = 0;
449 evt.DESYtriggerNameMap = 0;
450 evt.electronCharge = 0;
451 evt.electronP4 = 0;
452 evt.Event = 0;
453 evt.HT = 0.;
454 evt.jetP4 = 0;
455 evt.MET = 0.;
456 evt.muonCharge = 0;
457 evt.muonP4 = 0;
458 evt.prescaled = 0;
459 evt.Run = 0;
460 evt.PUInter = 0.;
461 evt.PUWeight = 0.;
462 evt.triggered = 0;
463 evt.Weight = 0.;
464 //--------------------------------------------------- assign branches to event structure
465 tree->SetBranchAddress("DESYtriggerElMatchedTriggerFilter", &(evt.DESYtriggerElMatchedTriggerFilter));
466 tree->SetBranchAddress("DESYtriggerMuMatchedTriggerFilter",&(evt.DESYtriggerMuMatchedTriggerFilter));
467 tree->SetBranchAddress("DESYtriggerNameMap",&(evt.DESYtriggerNameMap));
468 tree->SetBranchAddress("electronCharge",&(evt.electronCharge));
469 tree->SetBranchAddress("electronP4",&(evt.electronP4));
470 tree->SetBranchAddress("Event",&(evt.Event));
471 tree->SetBranchAddress("HT",&(evt.HT));
472 tree->SetBranchAddress("jetP4",&(evt.jetP4));
473 tree->SetBranchAddress("MET",&(evt.MET));
474 tree->SetBranchAddress("muonCharge",&(evt.muonCharge));
475 tree->SetBranchAddress("muonP4",&(evt.muonP4));
476 tree->SetBranchAddress("prescaled",&(evt.prescaled));
477 tree->SetBranchAddress("Run",&(evt.Run));
478 tree->SetBranchAddress("PUInter",&(evt.PUInter));
479 tree->SetBranchAddress("PUWeight",&(evt.PUWeight));
480 tree->SetBranchAddress("triggered",&(evt.triggered));
481 tree->SetBranchAddress("Weight",&(evt.Weight));
482 }
483
484
485 //================================================================================= printEventInfo
486 void printEventInfo(eventInfo& evt, ostream& os){
487 os<<endl;
488 os<<"///////////////////////////////////////////////////////"<<endl;
489 os<<"///////////////////// Event Info //////////////////////"<<endl;
490 os<<"///////////////////////////////////////////////////////"<<endl;
491 os<<endl;
492 os<<"-----------------------------------------"<<endl;
493 os<<"Run = "<<evt.Run<<endl;
494 os<<"Event = "<<evt.Event<<endl;
495 os<<"Weight = "<<evt.Weight<<endl;
496 os<<"PUWeight = "<<evt.PUWeight<<endl;
497 os<<"PUInter = "<<evt.PUInter<<endl;
498 os<<"-----------------------------------------"<<endl;
499 os<<"HT = "<<evt.HT<<endl;
500 os<<"MET = "<<evt.MET<<endl;
501 os<<"-----------------------------------------"<<endl;
502 os<<"No. of Muons = "<<evt.muonP4->size()<<endl;
503 os<<"No. of Electrons = "<<evt.electronP4->size()<<endl;
504 os<<"No. of Jets = "<<evt.jetP4->size()<<endl;
505 os<<"-----------------------------------------"<<endl;
506 os<<"muonCharge.size() = "<<evt.muonCharge->size()<<endl;
507 os<<"electronCharge.size() = "<<evt.electronCharge->size()<<endl;
508 os<<"-----------------------------------------"<<endl;
509 os<<"MuMatchedTriggerFilter.size() = "<<evt.DESYtriggerMuMatchedTriggerFilter->size()<<endl;
510 os<<"ElMatchedTriggerFilter.size() = "<<evt.DESYtriggerElMatchedTriggerFilter->size()<<endl;
511 os<<"-----------------------------------------"<<endl;
512 os<<"triggered.size() = "<<evt.triggered->size()<<endl;
513 os<<"prescaled.size() = "<<evt.prescaled->size()<<endl;
514 os<<"DESYTriggerMap.size() = "<<evt.DESYtriggerNameMap->size()<<endl;
515 os<<"-----------------------------------------"<<endl;
516 os<<endl;
517
518 //tools::printVector<string>("jetP4", *(evt.jetP4));
519 //tools::printVector<string>("DESYtriggerMuMatchedTriggerFilter", *(evt.DESYtriggerMuMatchedTriggerFilter));
520 //tools::printVector<string>("DESYtriggerElMatchedTriggerFilter", *(evt.DESYtriggerElMatchedTriggerFilter));
521 //tools::printVector<string>("DESYtriggerNameMap", *(evt.DESYtriggerNameMap));
522 //tools::printVector<string>("prescaled", *(evt.prescaled));
523 //tools::printVector<string>("triggered", *(evt.triggered));
524 }
525
526
527 //================================================================================= printMuonInfo
528 void printMuonInfo(eventInfo& evt, ostream& os){
529 os<<endl;
530 os<<"///////////////////////////////////////////////////////"<<endl;
531 os<<"////////////////////// Muon Info //////////////////////"<<endl;
532 os<<"///////////////////////////////////////////////////////"<<endl;
533 os<<endl;
534 tools::printVector<LorentzM>("muonP4", *(evt.muonP4), os);
535 tools::printVector<int>("muonCharge", *(evt.muonCharge), os);
536 //------------------------ print invariant mass from all possible muon pairs
537 if(evt.muonP4->size()>=2){
538 os<<"-------------> Invariant masses of all possible muon pairs:"<<endl;
539 for(Int_t i=0,N=evt.muonP4->size(); i<N-1; ++i){
540 for(Int_t j=i+1; j<N; ++j){
541 Double_t minv = (evt.muonP4->at(i) + evt.muonP4->at(j)).mass();
542 os<<"Minv("<<i<<" + " <<j<<") = "<<minv<<endl;
543 }
544 }
545 os<<endl;
546 }
547 //------------------------ print matching info
548 for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
549 os<<"----------------------------------------------------------------------------> Matches for muon "<<i<<":"<<
550 endl;
551 tools::printVector<TString>("matched TriggerFilter", evt.triggerFilterMu.at(i), os);
552 }
553 os<<endl;
554 }
555
556
557 //================================================================================= printElectronInfo
558 void printElectronInfo(eventInfo& evt, ostream& os){
559 os<<endl;
560 os<<"///////////////////////////////////////////////////////"<<endl;
561 os<<"///////////////////// Electron Info ///////////////////"<<endl;
562 os<<"///////////////////////////////////////////////////////"<<endl;
563 os<<endl;
564 tools::printVector<LorentzM>("electronP4", *(evt.electronP4), os);
565 tools::printVector<int>("electronCharge", *(evt.electronCharge), os);
566 //------------------------ print invariant mass from all possible electron pairs
567 if(evt.electronP4->size()>=2){
568 os<<"-------------> Invariant masses of all possible electron pairs:"<<endl;
569 for(Int_t i=0,N=evt.electronP4->size(); i<N-1; ++i){
570 for(Int_t j=i+1; j<N; ++j){
571 Double_t minv = (evt.electronP4->at(i) + evt.electronP4->at(j)).mass();
572 os<<"Minv("<<i<<" + " <<j<<") = "<<minv<<endl;
573 }
574 }
575 os<<endl;
576 }
577 //------------------------ print matching info
578 for(Int_t i=0,N=evt.electronP4->size(); i<N; ++i){
579 os<<"----------------------------------------------------------------------------> Matches for electron "<<i<<":"<<
580 endl;
581 tools::printVector<TString>("matched TriggerFilter", evt.triggerFilterEl.at(i), os);
582 }
583 os<<endl;
584 }
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600 //Thats another version of implementing the T&P method
601 //This version uses ALL Objects within an event to extract the trigger efficiency.
602 //However, this method leads to a more complex treatment of the counting and the error propagation.
603 //The advantage is higher statistict. But the statistis is anayhow large enough...
604 // //================================================================================= investigateTriggerFilterMu
605 // // eff = 2*NTT / (2*NTT + NTP)
606 // void investigateTriggerFilterMu(eventInfo& evt, map<TString, plotMaker> &plot, TString filterName){
607 // //create vector which stores, if a muon is matched to the trigger filter with the name "filterName"
608 // vector<bool> matched;
609 // vector<int> probeIndex;
610 // for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
611 // matched.push_back(false);
612 // for(Int_t j=0,M=evt.triggerFilterMu.at(i).size(); j<M; ++j){
613 // if(filterName == evt.triggerFilterMu.at(i).at(j)){
614 // matched.back() = true;
615 // break;
616 // }
617 // }
618 // if(!matched.at(i)) probeIndex.push_back(i);
619 // }
620 //
621 // if(probeIndex.size()>2){
622 // cout<<"WARNING: Something is awkward with the number of Probe leptons: " << probeIndex.size() << endl;
623 // cout<<"The number of probe leptons is expected to be 1 or 2, because here a T&P method is applied!" << endl;
624 // }
625 // else if(probeIndex.size()==2) return; //make only plots if at least one Tag exists in an event
626 //
627 // //Fill the Tag histograms
628 // //here, we have TWO Tag leptons -> use BOTH for the efficiency
629 // if(probeIndex.size()==0){
630 // for(Int_t i=0,N=evt.muonP4->size(); i<N; ++i){
631 // plot[filterName].addh1d("pt_tag", "pt_tag", "p_{T}(tag) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(i).Pt(), evt.Weight);
632 // plot[filterName].addh1d("eta_tag", "eta_tag", "#eta(tag)", "events", binEta, lowEta, upEta, evt.muonP4->at(i).eta(), evt.Weight);
633 // }
634 // //the factor of '2' is due to the fact, that there are two matched tag leptons
635 // plot[filterName].addh1d("PUVertices_tag", "PUVertices_tag", "number of primary vertices (tag)", "events", binPU, lowPU, upPU, evt.PUInter, 2*evt.Weight);
636 // }
637 //
638 // //Fill the Probe histograms
639 // //here, we have ONE Tag lepton -> the probe is NOT matched to the trigger filter
640 // if(probeIndex.size()==1){
641 // plot[filterName].addh1d("pt_probe", "pt_probe", "p_{T}(probe) [GeV/c]", "events", binPt, lowPt, upPt, evt.muonP4->at(probeIndex.back()).pt(), evt.Weight);
642 // plot[filterName].addh1d("eta_probe", "eta_probe", "#eta(probe)", "events", binEta, lowEta, upEta, evt.muonP4->at(probeIndex.back()).eta(), evt.Weight);
643 // plot[filterName].addh1d("PUVertices_probe", "PUVertices_probe", "number of primary vertices (probe)", "events", binPU, lowPU, upPU, evt.PUInter, evt.Weight);
644 // }
645 // }
646 //
647 //
648 // //================================================================================= makeEfficiencyPlots
649 // //eff = tag/(tag+probe)
650 // void makeEfficiencyPlots(map<TString, plotMaker> &plot){
651 // cout<<"----------> Create efficiency histograms"<<endl;
652 //
653 // vector<string> allTriggerFiltersToCheck = triggerFilterToCheckMu;
654 // for(Int_t i=0,N=triggerFilterToCheckEl.size(); i<N; ++i){
655 // allTriggerFiltersToCheck.push_back(triggerFilterToCheckEl.at(i));
656 // }
657 // vector<string> histNames;
658 // histNames.push_back("pt_tag");
659 // histNames.push_back("eta_tag");
660 // histNames.push_back("PUVertices_tag");
661 //
662 // for(Int_t i=0,N=allTriggerFiltersToCheck.size(); i<N; ++i){
663 // TString triggerFilter = allTriggerFiltersToCheck.at(i);
664 // for(Int_t j=0,M=histNames.size(); j<M; ++j){
665 // TString tagName = histNames.at(j);
666 // TString probeName = histNames.at(j);
667 // TString effName = histNames.at(j);
668 // probeName.ReplaceAll("tag", "probe");
669 // effName.ReplaceAll("tag", "eff");
670 // if(plot[triggerFilter].isNewh1d(tagName)) continue;
671 // if(plot[triggerFilter].isNewh1d(probeName)) continue;
672 //
673 // TH1D *h_tag = (TH1D*)(plot[triggerFilter].h1d[tagName])->Clone();
674 // TH1D *h_probe = (TH1D*)(plot[triggerFilter].h1d[probeName])->Clone();
675 // TH1D *den = (TH1D*)(h_tag)->Clone();
676 // den->Add((TH1D*)(h_probe)->Clone());
677 //
678 // TH1D *eff = (TH1D*)(h_tag)->Clone();
679 // eff->Divide(den);
680 // setEfficiencyErrors(eff, h_tag, h_probe);
681 //
682 // eff->SetName(effName);
683 // eff->SetTitle(effName);
684 //
685 // TString xTitle_eff = eff->GetXaxis()->GetTitle();
686 // xTitle_eff.ReplaceAll("(tag)","");
687 // eff->GetXaxis()->SetTitle(xTitle_eff);
688 // eff->GetYaxis()->SetTitle("efficiency");
689 // //eff->SetAxisRange(0.,1.,"Y");
690 // plot[triggerFilter].addh1d(eff);
691 // }
692 // }
693 // }
694 //
695 //
696 // //================================================================================= setEfficiencyErrors
697 // // f(eff) = t/(t+p)
698 // // df/dt = 1/(t+p) - t/(t+p)^2
699 // // df/dp = - t/(t+p)^2
700 // // eff_sig = sqrt(t_sig^2*(df/dt)^2 + p_sig^2*(df/dp)^2)
701 // void setEfficiencyErrors(TH1D* &eff, TH1D* &h_tag, TH1D* &h_probe){
702 // if(h_tag->GetNbinsX() != h_probe->GetNbinsX()){
703 // cout<<"WARNING: the number of bins of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
704 // cout<<"Number of bins (Tag | Probe): " << h_tag->GetNbinsX() << " | " << h_probe->GetNbinsX() << endl;
705 // }
706 // if(h_tag->GetXaxis()->GetXmax() != h_probe->GetXaxis()->GetXmax()){
707 // cout<<"WARNING: the upper x-range of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
708 // cout<<"Number of bins (Tag | Probe): " << h_tag->GetXaxis()->GetXmax() << " | " << h_probe->GetXaxis()->GetXmax() << endl;
709 // }
710 // if(h_tag->GetXaxis()->GetXmin() != h_probe->GetXaxis()->GetXmin()){
711 // cout<<"WARNING: the lower x-range of the tag and probe histograms differ! The efficiency calculation might be depricated!"<<endl;
712 // cout<<"Number of bins (Tag | Probe): " << h_tag->GetXaxis()->GetXmin() << " | " << h_probe->GetXaxis()->GetXmin() << endl;
713 // }
714 // //loop over single bins of histograms an set the efficiency error
715 // for(Int_t i=1,N=h_tag->GetNbinsX(); i<=N; ++i){
716 // Double_t t = h_tag ->GetBinContent(i);
717 // Double_t p = h_probe->GetBinContent(i);
718 // Double_t t_sig = h_tag ->GetBinError(i);
719 // Double_t p_sig = h_probe->GetBinError(i);
720 //
721 // if(eff->GetBinContent(i)==0.) continue;
722 //
723 // Double_t dedt = 1/(t+p)-t/((t+p)*(t+p));
724 // Double_t dedp = -t/((t+p)*(t+p));
725 //
726 // Double_t eff_sig = sqrt(t_sig*t_sig*dedt*dedt + p_sig*p_sig*dedp*dedp);
727 // eff->SetBinError(i, eff_sig);
728 // }
729 // }
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755