ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/TreeMod/src/AnaFwkMod.cc
Revision: 1.23
Committed: Mon Aug 26 22:43:59 2013 UTC (11 years, 8 months ago) by bendavid
Content type: text/plain
Branch: MAIN
Changes since 1.22: +12 -2 lines
Log Message:
Add histogram to track total MC weight, needed for proper normalization of sherpa samples

File Contents

# Content
1 // $Id: AnaFwkMod.cc,v 1.22 2012/10/03 18:13:59 ceballos Exp $
2
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
13 using namespace mithep;
14
15 ClassImp(mithep::AnaFwkMod)
16
17 //--------------------------------------------------------------------------------------------------
18 AnaFwkMod::AnaFwkMod(const char *name, const char *title) :
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),
27 fAllHeadTree(0),
28 fAllEventHeader(0),
29 fReload(kFALSE),
30 fCurEnt(-2),
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 {
42 // Constructor.
43 }
44
45 //--------------------------------------------------------------------------------------------------
46 void AnaFwkMod::BeginRun()
47 {
48 // Get HLT tree and set branches if new file was opened. Read next entry in HLT key
49 // depending on entry in RunInfo.
50
51 if (fReload) {
52 // reset to be (re-)loaded variables
53 fReload = 0;
54 fAllHeadTree = 0;
55 fAllEventHeader = 0;
56 fCurEnt = 0;
57
58 // get current file
59 TFile *file = GetCurrentFile();
60 if (!file)
61 return;
62
63 // get all event header tree
64 fAllHeadTree = dynamic_cast<TTree*>(file->Get(fAllHeadTreeName));
65 if (!fAllHeadTree) {
66 SendError(kWarning, "BeginRun",
67 "Cannot find tree '%s' in file '%s'",
68 fAllHeadTreeName.Data(),file->GetName());
69 return;
70 }
71
72 // get all event header branch
73 if (fAllHeadTree->GetBranch(fAllHeadBrName)) {
74 fAllHeadTree->SetBranchAddress(fAllHeadBrName, &fAllEventHeader);
75 } else {
76 SendError(kWarning, "BeginRun",
77 "Cannot find branch '%s' in tree '%s'",
78 fAllHeadBrName.Data(), fAllHeadTreeName.Data());
79 return;
80 }
81 }
82 }
83
84 //--------------------------------------------------------------------------------------------------
85 void AnaFwkMod::CopyAllEventHeaders()
86 {
87 // Deal with the headers from all events if needed.
88
89 fAllHeaders.Reset();
90
91 const EventHeader *curev = GetEventHeader();
92 if (!curev) {
93 SendError(kWarning, "CopyAllEventHeaders", "Cannot obtain current event");
94 return;
95 }
96
97 if (fAllHeadTree) {
98 const Int_t nemax = fAllHeadTree->GetEntries();
99 if (fCurEnt == nemax) {
100 SendError(kAbortEvent, "CopyAllEventHeaders",
101 "End of all events tree reached (%d=%d)", fCurEnt, nemax);
102 return;
103 }
104 fAllHeadTree->GetEntry(fCurEnt++);
105 while (fCurEnt <= nemax && fAllEventHeader->Skimmed()) {
106 EventHeader *eh = fAllHeaders.AddNew();
107 eh->SetRunNum (fAllEventHeader->RunNum());
108 eh->SetEvtNum (fAllEventHeader->EvtNum());
109 eh->SetLumiSec (fAllEventHeader->LumiSec());
110 eh->SetRunEntry(fAllEventHeader->RunEntry());
111 eh->SetSkimmed (fAllEventHeader->Skimmed());
112 fAllHeadTree->GetEntry(fCurEnt++);
113 }
114 if ((fAllEventHeader->RunNum() != curev->RunNum()) ||
115 (fAllEventHeader->EvtNum() != curev->EvtNum()) ||
116 (fAllEventHeader->LumiSec() != curev->LumiSec()) ||
117 (fAllEventHeader->RunEntry() != curev->RunEntry())) {
118 SendError(kWarning, "CopyAllEventHeaders",
119 "Event header information for entry %d inconsistent: "
120 "%d==%d, %d==%d, %d==%d, %d==%d",
121 fCurEnt,
122 fAllEventHeader->RunNum(), curev->RunNum(),
123 fAllEventHeader->EvtNum(), curev->EvtNum(),
124 fAllEventHeader->LumiSec(), curev->LumiSec(),
125 fAllEventHeader->RunEntry(), curev->RunEntry());
126 return;
127 }
128
129 // read-ahead to check if more events are coming
130 if (fCurEnt<nemax) {
131 Int_t testEnt = fCurEnt;
132 fAllHeadTree->GetEntry(testEnt++);
133 while (testEnt<=nemax && fAllEventHeader->Skimmed())
134 fAllHeadTree->GetEntry(testEnt++);
135 if (testEnt==nemax+1) { // need to add remaining skimmed events
136 fAllHeadTree->GetEntry(fCurEnt++);
137 while(fCurEnt<=nemax) {
138 EventHeader *eh = fAllHeaders.AddNew();
139 eh->SetRunNum (fAllEventHeader->RunNum());
140 eh->SetEvtNum (fAllEventHeader->EvtNum());
141 eh->SetLumiSec (fAllEventHeader->LumiSec());
142 eh->SetRunEntry(fAllEventHeader->RunEntry());
143 eh->SetSkimmed (fAllEventHeader->Skimmed());
144 fAllHeadTree->GetEntry(fCurEnt++);
145 }
146 if (fCurEnt != nemax+1) {
147 SendError(kAbortEvent, "CopyAllEventHeaders",
148 "End of all events tree unexpectedly not reached (%d!=%d)", fCurEnt, nemax);
149 return;
150 }
151 }
152 }
153 }
154 }
155
156 //--------------------------------------------------------------------------------------------------
157 Bool_t AnaFwkMod::Notify()
158 {
159 // Make sure to get the new "AllEvents" tree when the file changes.
160
161 fReload = kTRUE;
162 return kTRUE;
163 }
164
165 //--------------------------------------------------------------------------------------------------
166 void AnaFwkMod::Process()
167 {
168 // Do event skipping and counting and print out timing information.
169
170 // counting processed events
171 IncNEventsProcessed();
172
173 // get skimmed event headers
174 CopyAllEventHeaders();
175 fNEventsSkimmed += fAllHeaders.GetEntries();
176
177 // check if events should be skipped
178 if (fNEventsSkipped<fSkipNEvents) {
179 ++fNEventsSkipped;
180 MDB(kAnalysis, 3) {
181 Info("Process", "Skipping (aborting) %d out of %lld first events.",
182 fNEventsSkipped, fSkipNEvents);
183 }
184 AbortEvent();
185 return;
186 }
187
188 // check if printout should be done
189 Bool_t doPrint = 0;
190 UInt_t nProcessed = fPrintScale;
191
192 MDB(kAnalysis, 4) {
193 if (GetNEventsProcessed() % (fPrintScale) == 0)
194 doPrint = 1;
195 } else {
196 MDB(kAnalysis, 3) {
197 if (GetNEventsProcessed() % (fPrintScale*10) == 0)
198 doPrint = 1;
199 nProcessed = fPrintScale*10;
200 } else {
201 MDB(kAnalysis, 2) {
202 if (GetNEventsProcessed() % (fPrintScale*100) == 0)
203 doPrint = 1;
204 nProcessed = fPrintScale*100;
205 } else {
206 MDB(kAnalysis, 1) {
207 if (GetNEventsProcessed() % (fPrintScale*1000) == 0)
208 doPrint = 1;
209 nProcessed = fPrintScale*1000;
210 }
211 }
212 }
213 }
214
215 if (GetEventHeader()->IsMC()) {
216 LoadBranch(fPileupInfoName);
217 Double_t npu[4] = {0.,0.,0.,0.};
218 for (UInt_t i=0; i<fPileupInfo->GetEntries(); ++i) {
219 const PileupInfo *puinfo = fPileupInfo->At(i);
220 if (puinfo->GetBunchCrossing()==0) npu[0]= puinfo->GetPU_NumInteractions();
221 else if (puinfo->GetBunchCrossing()==-1) npu[1] = puinfo->GetPU_NumInteractions();
222 else if (puinfo->GetBunchCrossing()==1) npu[2] = puinfo->GetPU_NumInteractions();
223 if (puinfo->GetBunchCrossing()==0) npu[3]= puinfo->GetPU_NumMean();
224 }
225 hNPU->Fill(npu[0]);
226 hNPU50ns->Fill(npu[0],npu[1],npu[2]);
227 hNPUTrue->Fill(npu[3]);
228
229 LoadBranch(fMCEventInfoName);
230 hDTotalMCWeight->Fill(0., fMCEventInfo->Weight());
231
232 }
233
234
235 if (doPrint) {
236 fSWevent->Stop();
237 Info("Process",
238 "Events %d -> %.2gs real, %.2gs cpu (%.2g real, %.2g cpu per event)",
239 GetNEventsProcessed(), fSWevent->RealTime(), fSWevent->CpuTime(),
240 fSWevent->RealTime()/nProcessed,
241 fSWevent->CpuTime()/nProcessed);
242 fSWevent->Start();
243 }
244
245
246 }
247
248 //--------------------------------------------------------------------------------------------------
249 void AnaFwkMod::SlaveBegin()
250 {
251 // Book our histogram and start the stop watches.
252
253 fSWtotal = new TStopwatch;
254 fSWevent = new TStopwatch;
255
256 if (!PublishObj(&fAllHeaders)) {
257 SendError(kAbortAnalysis, "SlaveBegin",
258 "Could not publish all event headers with name %s.", fAllHeaders.GetName());
259 return;
260 }
261
262 ReqBranch(fPileupInfoName, fPileupInfo);
263
264 hNPU = new TH1D("hNPU", "hNPU", 201, -0.5, 200.5);
265 AddOutput(hNPU);
266
267 hNPU50ns = new TH3D("hNPU50ns", "hNPU50ns", 201, -0.5, 200.5, 201, -0.5, 200.5, 201, -0.5, 200.5);
268 AddOutput(hNPU50ns);
269
270 hNPUTrue = new TH1D("hNPUTrue", "hNPUTrue", 2000, 0.0, 200.0);
271 AddOutput(hNPUTrue);
272
273 hDTotalMCWeight = new TH1D("hDTotalMCWeight","hDTotalMCWeight",1,-0.5,0.5);
274 AddOutput(hDTotalMCWeight);
275
276 }
277
278 //--------------------------------------------------------------------------------------------------
279 void AnaFwkMod::SlaveTerminate()
280 {
281 // Fill event histogram and printout timing information.
282
283 RetractObj(fAllHeaders.GetName());
284
285 SaveNEventsProcessed();
286 TH1D *hDAllEvents = new TH1D("hDAllEvents","Sum of processed and skimmed events",1,-0.5,0.5);
287 hDAllEvents->Fill(0.0,fNEventsSkimmed+GetNEventsProcessed());
288 hDAllEvents->SetEntries(fNEventsSkimmed+GetNEventsProcessed());
289 AddOutput(hDAllEvents);
290
291 TH1D *hDSkippedEvents = new TH1D("hDSkippedEvents","Number of skipped events",1,-0.5,0.5);
292 hDSkippedEvents->Fill(0.0,fNEventsSkipped);
293 hDSkippedEvents->SetEntries(fNEventsSkipped);
294 AddOutput(hDSkippedEvents);
295
296 fSWtotal->Stop();
297 fSWevent->Stop();
298
299 MDB(kAnalysis, 1)
300 Info("SlaveTerminate",
301 "Events %d -> %.2gs real, %.2gs cpu (%.2gs real, %.2gs cpu per event)",
302 GetNEventsProcessed(), fSWtotal->RealTime(), fSWtotal->CpuTime(),
303 fSWtotal->RealTime()/GetNEventsProcessed(),
304 fSWtotal->CpuTime()/GetNEventsProcessed());
305
306 delete fSWtotal;
307 delete fSWevent;
308 }