3 |
|
#include "MitAna/TreeMod/interface/AnaFwkMod.h" |
4 |
|
#include "MitAna/DataUtil/interface/Debug.h" |
5 |
|
#include "MitAna/DataTree/interface/Names.h" |
6 |
+ |
#include "MitAna/DataTree/interface/PileupInfo.h" |
7 |
|
#include <TFile.h> |
8 |
|
#include <TH1D.h> |
9 |
+ |
#include <TH3D.h> |
10 |
|
#include <TStopwatch.h> |
11 |
|
#include <TTree.h> |
12 |
|
|
20 |
|
fAllHeadTreeName(Names::gkAllEvtTreeName), |
21 |
|
fAllHeadBrName(Names::gkAllEvtHeaderBrn), |
22 |
|
fSkipNEvents(0), |
23 |
+ |
fPrintScale(1), |
24 |
|
fSWtotal(0), |
25 |
|
fSWevent(0), |
26 |
|
fAllHeaders(0,Names::gkSkimmedHeaders), |
29 |
|
fReload(kFALSE), |
30 |
|
fCurEnt(-2), |
31 |
|
fNEventsSkimmed(0), |
32 |
< |
fNEventsSkipped(0) |
32 |
> |
fNEventsSkipped(0), |
33 |
> |
fPileupInfoName("PileupInfo"), |
34 |
> |
fDoPUInfo(kFALSE), |
35 |
> |
hNPU(0), |
36 |
> |
hNPU50ns(0), |
37 |
> |
hNPUTrue(0), |
38 |
> |
fMCEventInfo(0), |
39 |
> |
fMCEventInfoName(Names::gkMCEvtInfoBrn), |
40 |
> |
hDTotalMCWeight(0), |
41 |
> |
hNPURunABObs(0), |
42 |
> |
hNPURunCObs(0), |
43 |
> |
hNPURunDObs(0), |
44 |
> |
hNPURunABTrue(0), |
45 |
> |
hNPURunCTrue(0), |
46 |
> |
hNPURunDTrue(0) |
47 |
|
{ |
48 |
|
// Constructor. |
49 |
|
} |
66 |
|
if (!file) |
67 |
|
return; |
68 |
|
|
69 |
< |
// get HLT tree |
69 |
> |
// get all event header tree |
70 |
|
fAllHeadTree = dynamic_cast<TTree*>(file->Get(fAllHeadTreeName)); |
71 |
|
if (!fAllHeadTree) { |
72 |
|
SendError(kWarning, "BeginRun", |
74 |
|
fAllHeadTreeName.Data(),file->GetName()); |
75 |
|
return; |
76 |
|
} |
60 |
– |
} |
77 |
|
|
78 |
< |
// get HLT trigger name branch |
79 |
< |
if (fAllHeadTree->GetBranch(fAllHeadBrName)) { |
80 |
< |
fAllHeadTree->SetBranchAddress(fAllHeadBrName, &fAllEventHeader); |
81 |
< |
} else { |
82 |
< |
SendError(kWarning, "BeginRun", |
83 |
< |
"Cannot find branch '%s' in tree '%s'", |
84 |
< |
fAllHeadBrName.Data(), fAllHeadTreeName.Data()); |
85 |
< |
return; |
78 |
> |
// get all event header branch |
79 |
> |
if (fAllHeadTree->GetBranch(fAllHeadBrName)) { |
80 |
> |
fAllHeadTree->SetBranchAddress(fAllHeadBrName, &fAllEventHeader); |
81 |
> |
} else { |
82 |
> |
SendError(kWarning, "BeginRun", |
83 |
> |
"Cannot find branch '%s' in tree '%s'", |
84 |
> |
fAllHeadBrName.Data(), fAllHeadTreeName.Data()); |
85 |
> |
return; |
86 |
> |
} |
87 |
|
} |
88 |
|
} |
89 |
|
|
108 |
|
return; |
109 |
|
} |
110 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
111 |
< |
while(fCurEnt<=nemax && fAllEventHeader->Skimmed()) { |
111 |
> |
while (fCurEnt <= nemax && fAllEventHeader->Skimmed()) { |
112 |
|
EventHeader *eh = fAllHeaders.AddNew(); |
113 |
< |
eh->SetRunNum(fAllEventHeader->RunNum()); |
114 |
< |
eh->SetEvtNum(fAllEventHeader->EvtNum()); |
115 |
< |
eh->SetLumiSec(fAllEventHeader->LumiSec()); |
113 |
> |
eh->SetRunNum (fAllEventHeader->RunNum()); |
114 |
> |
eh->SetEvtNum (fAllEventHeader->EvtNum()); |
115 |
> |
eh->SetLumiSec (fAllEventHeader->LumiSec()); |
116 |
|
eh->SetRunEntry(fAllEventHeader->RunEntry()); |
117 |
< |
eh->SetSkimmed(fAllEventHeader->Skimmed()); |
117 |
> |
eh->SetSkimmed (fAllEventHeader->Skimmed()); |
118 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
119 |
|
} |
120 |
< |
if ((fAllEventHeader->RunNum()!=curev->RunNum()) || |
121 |
< |
(fAllEventHeader->EvtNum()!=curev->EvtNum()) || |
122 |
< |
(fAllEventHeader->LumiSec()!=curev->LumiSec()) || |
123 |
< |
(fAllEventHeader->RunEntry()!=curev->RunEntry())) { |
120 |
> |
if ((fAllEventHeader->RunNum() != curev->RunNum()) || |
121 |
> |
(fAllEventHeader->EvtNum() != curev->EvtNum()) || |
122 |
> |
(fAllEventHeader->LumiSec() != curev->LumiSec()) || |
123 |
> |
(fAllEventHeader->RunEntry() != curev->RunEntry())) { |
124 |
|
SendError(kWarning, "CopyAllEventHeaders", |
125 |
|
"Event header information for entry %d inconsistent: " |
126 |
|
"%d==%d, %d==%d, %d==%d, %d==%d", |
136 |
|
if (fCurEnt<nemax) { |
137 |
|
Int_t testEnt = fCurEnt; |
138 |
|
fAllHeadTree->GetEntry(testEnt++); |
139 |
< |
while(testEnt<=nemax && fAllEventHeader->Skimmed()) |
139 |
> |
while (testEnt<=nemax && fAllEventHeader->Skimmed()) |
140 |
|
fAllHeadTree->GetEntry(testEnt++); |
141 |
|
if (testEnt==nemax+1) { // need to add remaining skimmed events |
142 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
143 |
|
while(fCurEnt<=nemax) { |
144 |
|
EventHeader *eh = fAllHeaders.AddNew(); |
145 |
< |
eh->SetRunNum(fAllEventHeader->RunNum()); |
146 |
< |
eh->SetEvtNum(fAllEventHeader->EvtNum()); |
147 |
< |
eh->SetLumiSec(fAllEventHeader->LumiSec()); |
145 |
> |
eh->SetRunNum (fAllEventHeader->RunNum()); |
146 |
> |
eh->SetEvtNum (fAllEventHeader->EvtNum()); |
147 |
> |
eh->SetLumiSec (fAllEventHeader->LumiSec()); |
148 |
|
eh->SetRunEntry(fAllEventHeader->RunEntry()); |
149 |
< |
eh->SetSkimmed(fAllEventHeader->Skimmed()); |
149 |
> |
eh->SetSkimmed (fAllEventHeader->Skimmed()); |
150 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
151 |
|
} |
152 |
|
if (fCurEnt != nemax+1) { |
184 |
|
if (fNEventsSkipped<fSkipNEvents) { |
185 |
|
++fNEventsSkipped; |
186 |
|
MDB(kAnalysis, 3) { |
187 |
< |
Info("Process", "Skipping (aborting) %d out of %d first events.", |
187 |
> |
Info("Process", "Skipping (aborting) %d out of %lld first events.", |
188 |
|
fNEventsSkipped, fSkipNEvents); |
189 |
|
} |
190 |
|
AbortEvent(); |
192 |
|
} |
193 |
|
|
194 |
|
// check if printout should be done |
195 |
< |
Bool_t doPrint = 0; |
195 |
> |
Bool_t doPrint = 0; |
196 |
> |
UInt_t nProcessed = fPrintScale; |
197 |
|
|
198 |
|
MDB(kAnalysis, 4) { |
199 |
< |
if (GetNEventsProcessed() % 250 == 0) |
199 |
> |
if (GetNEventsProcessed() % (fPrintScale) == 0) |
200 |
|
doPrint = 1; |
201 |
|
} else { |
202 |
|
MDB(kAnalysis, 3) { |
203 |
< |
if (GetNEventsProcessed() % 2500 == 0) |
203 |
> |
if (GetNEventsProcessed() % (fPrintScale*10) == 0) |
204 |
|
doPrint = 1; |
205 |
+ |
nProcessed = fPrintScale*10; |
206 |
|
} else { |
207 |
|
MDB(kAnalysis, 2) { |
208 |
< |
if (GetNEventsProcessed() % 10000 == 0) |
208 |
> |
if (GetNEventsProcessed() % (fPrintScale*100) == 0) |
209 |
|
doPrint = 1; |
210 |
+ |
nProcessed = fPrintScale*100; |
211 |
|
} else { |
212 |
|
MDB(kAnalysis, 1) { |
213 |
< |
if (GetNEventsProcessed() % 50000 == 0) |
213 |
> |
if (GetNEventsProcessed() % (fPrintScale*1000) == 0) |
214 |
|
doPrint = 1; |
215 |
+ |
nProcessed = fPrintScale*1000; |
216 |
|
} |
217 |
|
} |
218 |
|
} |
219 |
|
} |
220 |
|
|
221 |
+ |
if (GetEventHeader()->IsMC()) { |
222 |
+ |
LoadBranch(fPileupInfoName); |
223 |
+ |
Double_t npu[4] = {0.,0.,0.,0.}; |
224 |
+ |
for (UInt_t i=0; i<fPileupInfo->GetEntries(); ++i) { |
225 |
+ |
const PileupInfo *puinfo = fPileupInfo->At(i); |
226 |
+ |
if (puinfo->GetBunchCrossing()==0) npu[0]= puinfo->GetPU_NumInteractions(); |
227 |
+ |
else if (puinfo->GetBunchCrossing()==-1) npu[1] = puinfo->GetPU_NumInteractions(); |
228 |
+ |
else if (puinfo->GetBunchCrossing()==1) npu[2] = puinfo->GetPU_NumInteractions(); |
229 |
+ |
if (puinfo->GetBunchCrossing()==0) npu[3]= puinfo->GetPU_NumMean(); |
230 |
+ |
} |
231 |
+ |
|
232 |
+ |
LoadBranch(fMCEventInfoName); |
233 |
+ |
double mcweight = fMCEventInfo->Weight(); |
234 |
+ |
|
235 |
+ |
hNPU->Fill(npu[0],mcweight); |
236 |
+ |
hNPU50ns->Fill(npu[0],npu[1],npu[2],mcweight); |
237 |
+ |
hNPUTrue->Fill(npu[3],mcweight); |
238 |
+ |
|
239 |
+ |
UInt_t run = GetEventHeader()->RunNum(); |
240 |
+ |
|
241 |
+ |
//Josh: fill pileup histograms for 2012 run-dependent Monte Carlo which currently has one of three run numbers |
242 |
+ |
//More generic solution to be implemented depending on future run/lumi dependent Monte Carlo production strategy |
243 |
+ |
if (run==194533) { |
244 |
+ |
hNPURunABObs->Fill(npu[0],mcweight); |
245 |
+ |
hNPURunABTrue->Fill(npu[3],mcweight); |
246 |
+ |
} |
247 |
+ |
else if (run==200519) { |
248 |
+ |
hNPURunCObs->Fill(npu[0],mcweight); |
249 |
+ |
hNPURunCTrue->Fill(npu[3],mcweight); |
250 |
+ |
} |
251 |
+ |
else if (run==206859) { |
252 |
+ |
hNPURunDObs->Fill(npu[0],mcweight); |
253 |
+ |
hNPURunDTrue->Fill(npu[3],mcweight); |
254 |
+ |
} |
255 |
+ |
|
256 |
+ |
hDTotalMCWeight->Fill(0., mcweight); |
257 |
+ |
|
258 |
+ |
} |
259 |
+ |
|
260 |
+ |
|
261 |
|
if (doPrint) { |
262 |
|
fSWevent->Stop(); |
263 |
|
Info("Process", |
264 |
|
"Events %d -> %.2gs real, %.2gs cpu (%.2g real, %.2g cpu per event)", |
265 |
|
GetNEventsProcessed(), fSWevent->RealTime(), fSWevent->CpuTime(), |
266 |
< |
fSWevent->RealTime()/GetNEventsProcessed(), |
267 |
< |
fSWevent->CpuTime()/GetNEventsProcessed()); |
266 |
> |
fSWevent->RealTime()/nProcessed, |
267 |
> |
fSWevent->CpuTime()/nProcessed); |
268 |
|
fSWevent->Start(); |
269 |
< |
} |
269 |
> |
} |
270 |
> |
|
271 |
> |
|
272 |
|
} |
273 |
|
|
274 |
|
//-------------------------------------------------------------------------------------------------- |
284 |
|
"Could not publish all event headers with name %s.", fAllHeaders.GetName()); |
285 |
|
return; |
286 |
|
} |
287 |
+ |
|
288 |
+ |
ReqBranch(fPileupInfoName, fPileupInfo); |
289 |
+ |
ReqBranch(fMCEventInfoName, fMCEventInfo); |
290 |
+ |
|
291 |
+ |
hNPU = new TH1D("hNPU", "hNPU", 201, -0.5, 200.5); |
292 |
+ |
AddOutput(hNPU); |
293 |
+ |
|
294 |
+ |
hNPU50ns = new TH3D("hNPU50ns", "hNPU50ns", 201, -0.5, 200.5, 201, -0.5, 200.5, 201, -0.5, 200.5); |
295 |
+ |
AddOutput(hNPU50ns); |
296 |
+ |
|
297 |
+ |
hNPUTrue = new TH1D("hNPUTrue", "hNPUTrue", 2000, 0.0, 200.0); |
298 |
+ |
AddOutput(hNPUTrue); |
299 |
+ |
|
300 |
+ |
hDTotalMCWeight = new TH1D("hDTotalMCWeight","hDTotalMCWeight",1,-0.5,0.5); |
301 |
+ |
AddOutput(hDTotalMCWeight); |
302 |
+ |
|
303 |
+ |
hNPURunABObs = new TH1D("hNPURunABObs","hNPURunABObs",201, -0.5, 200.5); |
304 |
+ |
AddOutput(hNPURunABObs); |
305 |
+ |
|
306 |
+ |
hNPURunCObs = new TH1D("hNPURunCObs","hNPURunCObs",201, -0.5, 200.5); |
307 |
+ |
AddOutput(hNPURunCObs); |
308 |
+ |
|
309 |
+ |
hNPURunDObs = new TH1D("hNPURunDObs","hNPURunDObs",201, -0.5, 200.5); |
310 |
+ |
AddOutput(hNPURunDObs); |
311 |
+ |
|
312 |
+ |
hNPURunABTrue = new TH1D("hNPURunABTrue","hNPURunABTrue",2000, 0.0, 200.0); |
313 |
+ |
AddOutput(hNPURunABTrue); |
314 |
+ |
|
315 |
+ |
hNPURunCTrue = new TH1D("hNPURunCTrue","hNPURunCTrue",2000, 0.0, 200.0); |
316 |
+ |
AddOutput(hNPURunCTrue); |
317 |
+ |
|
318 |
+ |
hNPURunDTrue = new TH1D("hNPURunDTrue","hNPURunDTrue",2000, 0.0, 200.0); |
319 |
+ |
AddOutput(hNPURunDTrue); |
320 |
+ |
|
321 |
|
} |
322 |
|
|
323 |
|
//-------------------------------------------------------------------------------------------------- |