ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitAna/TreeMod/src/AnaFwkMod.cc

Comparing UserCode/MitAna/TreeMod/src/AnaFwkMod.cc (file contents):
Revision 1.1 by loizides, Wed Nov 19 15:30:26 2008 UTC vs.
Revision 1.25 by bendavid, Fri Sep 13 10:16:35 2013 UTC

#	Line 2 | Line 2
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
#	Line 12 | Line 17 | ClassImp(mithep::AnaFwkMod)
17		//--------------------------------------------------------------------------------------------------
18		AnaFwkMod::AnaFwkMod(const char *name, const char *title) :
19		BaseMod(name,title),
20	<	fNEventsProcessed(0),
21	<	hDEvents(0),
20	>	fAllHeadTreeName(Names::gkAllEvtTreeName),
21	>	fAllHeadBrName(Names::gkAllEvtHeaderBrn),
22	>	fSkipNEvents(0),
23	>	fPrintScale(1),
24		fSWtotal(0),
25	<	fSWevent(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	>	hNPURunABObs(0),
42	>	hNPURunCObs(0),
43	>	hNPURunDObs(0),
44	>	hNPURunABTrue(0),
45	>	hNPURunCTrue(0),
46	>	hNPURunDTrue(0)
47		{
48		// Constructor.
49		}
50
51		//--------------------------------------------------------------------------------------------------
52	+	void AnaFwkMod::BeginRun()
53	+	{
54	+	// Get HLT tree and set branches if new file was opened. Read next entry in HLT key
55	+	// depending on entry in RunInfo.
56	+
57	+	if (fReload) {
58	+	// reset to be (re-)loaded variables
59	+	fReload         = 0;
60	+	fAllHeadTree    = 0;
61	+	fAllEventHeader = 0;
62	+	fCurEnt         = 0;
63	+
64	+	// get current file
65	+	TFile *file = GetCurrentFile();
66	+	if (!file)
67	+	return;
68	+
69	+	// get all event header tree
70	+	fAllHeadTree = dynamic_cast<TTree*>(file->Get(fAllHeadTreeName));
71	+	if (!fAllHeadTree) {
72	+	SendError(kWarning, "BeginRun",
73	+	"Cannot find tree '%s' in file '%s'",
74	+	fAllHeadTreeName.Data(),file->GetName());
75	+	return;
76	+	}
77	+
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	+
90	+	//--------------------------------------------------------------------------------------------------
91	+	void AnaFwkMod::CopyAllEventHeaders()
92	+	{
93	+	// Deal with the headers from all events if needed.
94	+
95	+	fAllHeaders.Reset();
96	+
97	+	const EventHeader *curev = GetEventHeader();
98	+	if (!curev) {
99	+	SendError(kWarning, "CopyAllEventHeaders", "Cannot obtain current event");
100	+	return;
101	+	}
102	+
103	+	if (fAllHeadTree) {
104	+	const Int_t nemax = fAllHeadTree->GetEntries();
105	+	if (fCurEnt == nemax) {
106	+	SendError(kAbortEvent, "CopyAllEventHeaders",
107	+	"End of all events tree reached (%d=%d)", fCurEnt, nemax);
108	+	return;
109	+	}
110	+	fAllHeadTree->GetEntry(fCurEnt++);
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());
116	+	eh->SetRunEntry(fAllEventHeader->RunEntry());
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())) {
124	+	SendError(kWarning, "CopyAllEventHeaders",
125	+	"Event header information for entry %d inconsistent: "
126	+	"%d==%d, %d==%d, %d==%d, %d==%d",
127	+	fCurEnt,
128	+	fAllEventHeader->RunNum(),   curev->RunNum(),
129	+	fAllEventHeader->EvtNum(),   curev->EvtNum(),
130	+	fAllEventHeader->LumiSec(),  curev->LumiSec(),
131	+	fAllEventHeader->RunEntry(), curev->RunEntry());
132	+	return;
133	+	}
134	+
135	+	// read-ahead to check if more events are coming
136	+	if (fCurEnt<nemax) {
137	+	Int_t testEnt = fCurEnt;
138	+	fAllHeadTree->GetEntry(testEnt++);
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());
148	+	eh->SetRunEntry(fAllEventHeader->RunEntry());
149	+	eh->SetSkimmed (fAllEventHeader->Skimmed());
150	+	fAllHeadTree->GetEntry(fCurEnt++);
151	+	}
152	+	if (fCurEnt != nemax+1) {
153	+	SendError(kAbortEvent, "CopyAllEventHeaders",
154	+	"End of all events tree unexpectedly not reached (%d!=%d)", fCurEnt, nemax);
155	+	return;
156	+	}
157	+	}
158	+	}
159	+	}
160	+	}
161	+
162	+	//--------------------------------------------------------------------------------------------------
163	+	Bool_t AnaFwkMod::Notify()
164	+	{
165	+	// Make sure to get the new "AllEvents" tree when the file changes.
166	+
167	+	fReload = kTRUE;
168	+	return kTRUE;
169	+	}
170	+
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 events
177	<	fNEventsProcessed++;
176	>	// counting processed events
177	>	IncNEventsProcessed();
178	>
179	>	// get skimmed event headers
180	>	CopyAllEventHeaders();
181	>	fNEventsSkimmed += fAllHeaders.GetEntries();
182	>
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 (fNEventsProcessed % 1000  == 0)
199	>	if (GetNEventsProcessed() % (fPrintScale)  == 0)
200		doPrint = 1;
201		} else {
202		MDB(kAnalysis, 3) {
203	<	if (fNEventsProcessed % 10000  == 0)
203	>	if (GetNEventsProcessed() % (fPrintScale*10)  == 0)
204		doPrint = 1;
205	+	nProcessed = fPrintScale*10;
206		} else {
207		MDB(kAnalysis, 2) {
208	<	if (fNEventsProcessed % 50000  == 0)
208	>	if (GetNEventsProcessed() % (fPrintScale*100)  == 0)
209		doPrint = 1;
210	+	nProcessed = fPrintScale*100;
211		} else {
212		MDB(kAnalysis, 1) {
213	<	if (fNEventsProcessed % 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)",
264	<	fNEventsProcessed, fSWevent->RealTime(), fSWevent->CpuTime(),
265	<	fSWevent->RealTime()/fNEventsProcessed, fSWevent->CpuTime()/fNEventsProcessed);
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()/nProcessed,
267	>	fSWevent->CpuTime()/nProcessed);
268		fSWevent->Start();
269	<	}
269	>	}
270	>
271	>
272		}
273
274		//--------------------------------------------------------------------------------------------------
#	Line 65 | Line 276 | void AnaFwkMod::SlaveBegin()
276		{
277		// Book our histogram and start the stop watches.
278
68	–	hDEvents = new TH1D("hDEvents","hDEvents for overall event count",1,-0.5,0.5);
69	–	AddOutput(hDEvents);
70	–
279		fSWtotal = new TStopwatch;
280		fSWevent = new TStopwatch;
281	+
282	+	if (!PublishObj(&fAllHeaders)) {
283	+	SendError(kAbortAnalysis, "SlaveBegin",
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		//--------------------------------------------------------------------------------------------------
#	Line 77 | Line 325 | void AnaFwkMod::SlaveTerminate()
325		{
326		// Fill event histogram and printout timing information.
327
328	<	hDEvents->Fill(0.0,fNEventsProcessed);
329	<	hDEvents->SetEntries(fNEventsProcessed);
328	>	RetractObj(fAllHeaders.GetName());
329	>
330	>	SaveNEventsProcessed();
331	>	TH1D *hDAllEvents = new TH1D("hDAllEvents","Sum of processed and skimmed events",1,-0.5,0.5);
332	>	hDAllEvents->Fill(0.0,fNEventsSkimmed+GetNEventsProcessed());
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)",
346	<	fNEventsProcessed, fSWtotal->RealTime(), fSWtotal->CpuTime(),
347	<	fSWtotal->RealTime()/fNEventsProcessed, fSWtotal->CpuTime()/fNEventsProcessed);
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(),
349	>	fSWtotal->CpuTime()/GetNEventsProcessed());
350
351		delete fSWtotal;
352		delete fSWevent;

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines