ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitAna/TreeMod/src/AnaFwkMod.cc
(Generate patch)

Comparing UserCode/MitAna/TreeMod/src/AnaFwkMod.cc (file contents):
Revision 1.4 by loizides, Mon Mar 9 17:07:05 2009 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  
# Line 14 | Line 17 | ClassImp(mithep::AnaFwkMod)
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,"AllEventHeaders"),
26 >  fAllHeaders(0,Names::gkSkimmedHeaders),
27    fAllHeadTree(0),
28    fAllEventHeader(0),
29    fReload(kFALSE),
30 <  fCurEnt(-2)
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   }
# Line 32 | Line 55 | void AnaFwkMod::BeginRun()
55    // depending on entry in RunInfo.
56  
57    if (fReload) {
35
58      // reset to be (re-)loaded variables
59      fReload         = 0;
60      fAllHeadTree    = 0;
# Line 44 | Line 66 | void AnaFwkMod::BeginRun()
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 +                "Cannot find tree '%s' in file '%s'",
74 +                fAllHeadTreeName.Data(),file->GetName());
75        return;
76      }
52  }
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 with name %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  
# Line 71 | Line 96 | void AnaFwkMod::CopyAllEventHeaders()
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  
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(fAllEventHeader->IsRemoved()) {
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->SetIsRemoved(fAllEventHeader->RunEntry());
117 <      eh->SetRunEntry(fAllEventHeader->IsRemoved());
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())) {
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 inconsistent: %d==%d, %d==%d, %d==%d, %d==%d",
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  
# Line 113 | Line 171 | Bool_t AnaFwkMod::Notify()
171   //--------------------------------------------------------------------------------------------------
172   void AnaFwkMod::Process()
173   {
174 <  // Do event counting and print out timing information.
117 <
118 <  
174 >  // Do event skipping and counting and print out timing information.
175  
176 <  // counting events
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 (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   //--------------------------------------------------------------------------------------------------
# Line 165 | Line 284 | void AnaFwkMod::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 175 | Line 328 | void AnaFwkMod::SlaveTerminate()
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)",
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;
353   }
190

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines