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 |
|
|
19 |
|
BaseMod(name,title), |
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), |
28 |
|
fAllEventHeader(0), |
29 |
|
fReload(kFALSE), |
30 |
|
fCurEnt(-2), |
31 |
< |
fNEventsSkimmed(0) |
31 |
> |
fNEventsSkimmed(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 |
|
} |
55 |
|
// depending on entry in RunInfo. |
56 |
|
|
57 |
|
if (fReload) { |
39 |
– |
|
58 |
|
// reset to be (re-)loaded variables |
59 |
|
fReload = 0; |
60 |
|
fAllHeadTree = 0; |
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", |
73 |
< |
"Can not find tree '%s' in file '%s'", |
73 |
> |
"Cannot find tree '%s' in file '%s'", |
74 |
|
fAllHeadTreeName.Data(),file->GetName()); |
75 |
|
return; |
76 |
|
} |
59 |
– |
} |
77 |
|
|
78 |
< |
// get HLT trigger name branch |
79 |
< |
if (fAllHeadTree->GetBranch(fAllHeadBrName)) { |
80 |
< |
fAllHeadTree->SetBranchAddress(fAllHeadBrName, &fAllEventHeader); |
81 |
< |
} else { |
82 |
< |
SendError(kWarning, "BeginRun", |
83 |
< |
"Can not 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 |
|
|
96 |
|
|
97 |
|
const EventHeader *curev = GetEventHeader(); |
98 |
|
if (!curev) { |
99 |
< |
SendError(kWarning, "CopyAllEventHeaders", "Can not obtain current event"); |
99 |
> |
SendError(kWarning, "CopyAllEventHeaders", "Cannot obtain current event"); |
100 |
|
return; |
101 |
|
} |
102 |
|
|
108 |
|
return; |
109 |
|
} |
110 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
111 |
< |
while(fAllEventHeader->Skimmed() && fCurEnt<nemax) { |
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(fAllEventHeader->Skimmed() && testEnt<nemax) |
139 |
> |
while (testEnt<=nemax && fAllEventHeader->Skimmed()) |
140 |
|
fAllHeadTree->GetEntry(testEnt++); |
141 |
< |
if (testEnt==nemax) { // need to add remaining skimmed events |
141 |
> |
if (testEnt==nemax+1) { // need to add remaining skimmed events |
142 |
|
fAllHeadTree->GetEntry(fCurEnt++); |
143 |
< |
while(fAllEventHeader->Skimmed() && fCurEnt<nemax) { |
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) { |
152 |
> |
if (fCurEnt != nemax+1) { |
153 |
|
SendError(kAbortEvent, "CopyAllEventHeaders", |
154 |
|
"End of all events tree unexpectedly not reached (%d!=%d)", fCurEnt, nemax); |
155 |
|
return; |
171 |
|
//-------------------------------------------------------------------------------------------------- |
172 |
|
void AnaFwkMod::Process() |
173 |
|
{ |
174 |
< |
// Do event counting and print out timing information. |
174 |
> |
// Do event skipping and counting and print out timing information. |
175 |
|
|
176 |
+ |
// counting processed events |
177 |
+ |
IncNEventsProcessed(); |
178 |
|
|
179 |
|
// get skimmed event headers |
180 |
|
CopyAllEventHeaders(); |
181 |
|
fNEventsSkimmed += fAllHeaders.GetEntries(); |
182 |
|
|
183 |
< |
// counting events |
184 |
< |
IncNEventsProcessed(); |
183 |
> |
// check if events should be skipped |
184 |
> |
if (fNEventsSkipped<fSkipNEvents) { |
185 |
> |
++fNEventsSkipped; |
186 |
> |
MDB(kAnalysis, 3) { |
187 |
> |
Info("Process", "Skipping (aborting) %d out of %lld first events.", |
188 |
> |
fNEventsSkipped, fSkipNEvents); |
189 |
> |
} |
190 |
> |
AbortEvent(); |
191 |
> |
return; |
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() % 1000 == 0) |
199 |
> |
if (GetNEventsProcessed() % (fPrintScale) == 0) |
200 |
|
doPrint = 1; |
201 |
|
} else { |
202 |
|
MDB(kAnalysis, 3) { |
203 |
< |
if (GetNEventsProcessed() % 10000 == 0) |
203 |
> |
if (GetNEventsProcessed() % (fPrintScale*10) == 0) |
204 |
|
doPrint = 1; |
205 |
+ |
nProcessed = fPrintScale*10; |
206 |
|
} else { |
207 |
|
MDB(kAnalysis, 2) { |
208 |
< |
if (GetNEventsProcessed() % 50000 == 0) |
208 |
> |
if (GetNEventsProcessed() % (fPrintScale*100) == 0) |
209 |
|
doPrint = 1; |
210 |
+ |
nProcessed = fPrintScale*100; |
211 |
|
} else { |
212 |
|
MDB(kAnalysis, 1) { |
213 |
< |
if (GetNEventsProcessed() % 250000 == 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", "Events %d -> %.2gs real, %.2gs cpu (%.2g real, %.2g cpu per event)", |
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(), 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 |
|
//-------------------------------------------------------------------------------------------------- |
333 |
|
hDAllEvents->SetEntries(fNEventsSkimmed+GetNEventsProcessed()); |
334 |
|
AddOutput(hDAllEvents); |
335 |
|
|
336 |
+ |
TH1D *hDSkippedEvents = new TH1D("hDSkippedEvents","Number of skipped events",1,-0.5,0.5); |
337 |
+ |
hDSkippedEvents->Fill(0.0,fNEventsSkipped); |
338 |
+ |
hDSkippedEvents->SetEntries(fNEventsSkipped); |
339 |
+ |
AddOutput(hDSkippedEvents); |
340 |
+ |
|
341 |
|
fSWtotal->Stop(); |
342 |
|
fSWevent->Stop(); |
343 |
|
|
344 |
|
MDB(kAnalysis, 1) |
345 |
< |
Info("SlaveTerminate", "Events %d -> %.2gs real, %.2gs cpu (%.2gs real, %.2gs cpu per event)", |
345 |
> |
Info("SlaveTerminate", |
346 |
> |
"Events %d -> %.2gs real, %.2gs cpu (%.2gs real, %.2gs cpu per event)", |
347 |
|
GetNEventsProcessed(), fSWtotal->RealTime(), fSWtotal->CpuTime(), |
348 |
< |
fSWtotal->RealTime()/GetNEventsProcessed(), fSWtotal->CpuTime()/GetNEventsProcessed()); |
348 |
> |
fSWtotal->RealTime()/GetNEventsProcessed(), |
349 |
> |
fSWtotal->CpuTime()/GetNEventsProcessed()); |
350 |
|
|
351 |
|
delete fSWtotal; |
352 |
|
delete fSWevent; |