Class for counting events and for pile-up correction factors (Diego)
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskPileup.cxx
CommitLineData
0fa651e5 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//-------------------------------------------------------------------------
17// Implementation of class AliAnalysisTaskPileup
18//
19// The class re-weights the number of analysed events with the correction
20// factor to account for the event pileup.
21//
22// The final idea will be to include all of the methods for the pileup
23// calculation.
24//
25// For the moment, the implemented method is based on a work
26// by L. Aphecetche.
27// The probability of having at least 1 collision is given by:
28// MB L0 trigger / max collision rate
29// The max collision rate is given by:
30// number of bunches X revolution frequency ~ CBEAMB L0 trigger
31// Assuming a poissonian distribution of the collision, with:
32// mu = mean number of collisions
33// the probability of having 0 collisions is:
34// exp{-mu} = 1 - L0b_MB / max collision rate
35// Since the MB trigger is measuring the probability of having at least
36// 1 collisions, the number should be corrected by the factor:
37// CF = mu / ( 1 - exp{-mu} )
38// The class weights the number of events with this correction factor
39//
40// CAVEATS:
41// - so far the class needs to access the OCDB
42// (hopefully it will be changed in the future)
43// Hence the following libraries need to be included (in addition to
44// the standard analysis ones):
45// libRAWDatabase.so libCDB.so libSTEER.so libPWG3base.so
46//-------------------------------------------------------------------------
47
48#include <Riostream.h>
49
50// ROOT includes
51#include "TH1.h"
52#include "TH2.h"
53#include "TCanvas.h"
54#include "TROOT.h"
55#include "TString.h"
56#include "TObjArray.h"
57#include "TObjString.h"
58#include "TMath.h"
59
60// STEER includes
61#include "AliESDEvent.h"
62#include "AliESDInputHandler.h"
63#include "AliAODEvent.h"
64#include "AliAODInputHandler.h"
65#include "AliTriggerScalersESD.h"
66#include "AliLog.h"
67
68// ANALYSIS includes
69#include "AliAnalysisTaskSE.h"
70#include "AliAnalysisManager.h"
71#include "AliAnalysisDataSlot.h"
72#include "AliAnalysisDataContainer.h"
73
74#include "AliCounterCollection.h"
75
76#include "AliAnalysisTaskPileup.h"
77
78#if defined(READOCDB)
79#include "AliTriggerRunScalers.h"
80#include "AliTriggerConfiguration.h"
81#include "AliTriggerClass.h"
82#include "AliCDBManager.h"
83#include "AliCDBEntry.h"
84#include "AliCDBPath.h"
85#endif
86
87
88ClassImp(AliAnalysisTaskPileup)
89
90
91//________________________________________________________________________
92AliAnalysisTaskPileup::AliAnalysisTaskPileup(const char *name) :
93 AliAnalysisTaskSE(name),
94 fEventCounters(0x0),
95 fHistoEventsList(0x0),
96 fTriggerClasses(0x0),
97 fTriggerClassIndex(new TArrayI(50))
98#if defined(READOCDB)
99 , fTriggerRunScalers(0x0),
100 fDefaultStorage(new TString())
101#endif
102
103{
104 /// Constructor
105
106 DefineOutput(1,AliCounterCollection::Class());
107 DefineOutput(2,TObjArray::Class());
108}
109
110//________________________________________________________________________
111AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
112{
113 /// Destructor
114 delete fEventCounters;
115 delete fHistoEventsList;
116 delete fTriggerClasses;
117 delete fTriggerClassIndex;
118
119#if defined(READOCDB)
120 delete fTriggerRunScalers;
121 delete fDefaultStorage;
122#endif
123
124}
125
126
127//___________________________________________________________________________
128void AliAnalysisTaskPileup::NotifyRun()
129{
130 /// Notify run
131
132#if defined(READOCDB)
133 if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) {
134 AliCDBManager::Instance()->SetDefaultStorage(fDefaultStorage->Data());
135 }
136
137 AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber());
138
139 AliCDBEntry *entry = 0x0;
140 entry = AliCDBManager::Instance()->Get("GRP/CTP/Config");
141 if ( entry ) {
142 AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject();
143 const TObjArray& classesArray = trigConf->GetClasses();
144 Int_t nclasses = classesArray.GetEntriesFast();
145
146 if ( fTriggerClasses ) delete fTriggerClasses;
147
148 fTriggerClasses = new TObjArray(nclasses);
149 fTriggerClasses->SetOwner();
150
151 Int_t currActive = -1;
152 for( Int_t iclass=0; iclass < nclasses; iclass++ ) {
153 AliTriggerClass* trclass = (AliTriggerClass*)classesArray.At(iclass);
154 if (trclass && trclass->GetMask()>0) {
155 currActive++;
156 Int_t currPos = currActive;
157
158 // Store the CBEAMB class at the first position
159 TString trigName = trclass->GetName();
160 if ( trigName.Contains("CBEAMB") && ! trigName.Contains("WU") )
161 currPos = 0;
162 else if ( ! fTriggerClasses->At(0) )
163 currPos++;
164
165 Int_t trindex = TMath::Nint(TMath::Log2(trclass->GetMask()));
166 TObjString* objString = new TObjString(trigName.Data());
167 fTriggerClasses->AddAtAndExpand(objString, currPos);
168 (*fTriggerClassIndex)[currPos] = trindex;
169 AliDebug(3, Form("Current class %s index %i position %i", trigName.Data(), trindex, currPos));
170 } // is class active
171 } // loop on trigger classes
172 } // if entry
173
174 entry = AliCDBManager::Instance()->Get("GRP/CTP/Scalers");
175 if (entry) {
176 AliInfo("Found an AliTriggerRunScalers in GRP/CTP/Scalers, reading it");
177 fTriggerRunScalers = dynamic_cast<AliTriggerRunScalers*> (entry->GetObject());
178 entry->SetOwner(0);
179 if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow");
180 }
181#endif
182
183}
184
185
186//___________________________________________________________________________
187void AliAnalysisTaskPileup::UserCreateOutputObjects()
188{
189 /// Create histograms and counters
190
191 // The framework has problems if the name of the object
192 // and the one of container differ
193 // To be complaint, get the name from container and set it
194 TString containerName = GetOutputSlot(1)->GetContainer()->GetName();
195
196 // initialize event counters
197 fEventCounters = new AliCounterCollection(containerName.Data());
198 fEventCounters->AddRubric("event", "any/correctedL0");
199 fEventCounters->AddRubric("trigger", 1000000);
200 fEventCounters->AddRubric("selection", "any/hasVtxContrib/nonPileupSPD");
201 fEventCounters->Init(kTRUE);
202
203 // Post data at least once per task to ensure data synchronisation (required for merging)
204 PostData(1, fEventCounters);
205}
206
207//________________________________________________________________________
208void AliAnalysisTaskPileup::UserExec(Option_t *)
209{
210 /// Called for each event
211
212 AliAODEvent* aodEvent = 0x0;
213
214 AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
215 if ( ! esdEvent ) {
216 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
217 }
218
219 if ( ! aodEvent && ! esdEvent ) {
220 AliError ("AOD or ESD event not found. Nothing done!");
221 return;
222 }
223
224 // check physics selection
225 //Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
226
227#if defined(READOCDB)
228
229 TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
230
231 if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events
232
233 AliDebug(2, Form("\nEvent %lli\n", Entry()));
234
235 Int_t nPoints = fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast();
236
237 AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber());
238 Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp);
239 if ( position < 0 ) {
240 AliWarning("Position out of range: put to 1");
241 position = 1;
242 }
243 if ( position == 0 ) position++; // Don't trust the first one
244 else if ( position + 1 >= nPoints ) position--;
245 AliDebug(2, Form("position %i\n", position));
246 AliTriggerScalersRecordESD* trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);
247 AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
248
249 // Sometimes scalers are filled very close to each others
250 // in this case skip to the next entry
251 for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) {
252 trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos);
253 Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() );
254 AliDebug(2, Form("Pos %i TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime));
255
256 if ( deltaTime > 1 )
257 break;
258 }
259
260#endif
261
262 ULong64_t trigMask = 0;
263 Double_t correctFactor = 1.;
264
265 Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
266 Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8);
267
268 TString selKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
269 Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )};
270
271 //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT
272
273 TString trigName = "";
274 TString eventType = "";
275
276 Int_t nTriggerClasses = fTriggerClasses->GetEntries();
277 Int_t classIndex = -1;
278 Double_t deltaScalersBeam = 0., deltaScalers = 0.;
279 Bool_t isFiredOnce = kFALSE;
280 for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) {
281
282 Double_t correctFactorL0 = 1.;
283
284 Bool_t isClassFired = kFALSE;
285
286 if ( itrig < nTriggerClasses ) {
287
288 // Check if current mask contains trigger
289 trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
290 classIndex = (*fTriggerClassIndex)[itrig];
291 trigMask = ( 1ull << classIndex );
292 isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
293
294 if ( isClassFired || itrig == 0 ) {
295 // Get scalers
296 const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1);
297 const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1);
298 deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB();
299
300 if ( itrig == 0 )
301 deltaScalersBeam = deltaScalers;
302 else if ( isClassFired ) {
303 correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam);
304 AliDebug(2, Form("Scalers: %s %.0f %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0));
305 }
306 }
307 }
308 else {
309 trigName = "any";
310 classIndex = -1;
311 isClassFired = isFiredOnce;
312 }
313
314 if ( ! isClassFired ) continue;
315 isFiredOnce = kTRUE;
316
317 //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT
318 //if ( classIndex > 1 ) printf("Index: trigger %i scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT
319
320 AliDebug(2, Form("Fired trig %s\n", trigName.Data()));
321
322 for ( Int_t ievType=0; ievType<2; ievType++ ){
323 switch ( ievType ) {
324#if defined(READOCDB)
325 case kHeventsCorrectL0:
326 correctFactor = correctFactorL0;
327 eventType = "correctedL0";
328 break;
329#endif
330 default:
331 correctFactor = 1.;
332 eventType = "any";
333 }
334
335 for ( Int_t isel=0; isel<3; isel++ ) {
336 if ( ! fillSel[isel] ) continue;
337 fEventCounters->Count(Form("event:%s/trigger:%s/selection:%s",eventType.Data(),trigName.Data(), selKey[isel].Data()),correctFactor);
338 } // loop on vertex selection
339 } // loop on event type
340 } // loop on trigger classes
341
342 // Post final data. It will be written to a file with option "RECREATE"
343 PostData(1, fEventCounters);
344}
345
346
347//________________________________________________________________________
348void AliAnalysisTaskPileup::Terminate(Option_t *)
349{
350 //
351 /// Save final histograms
352 /// and draw result to the screen
353 //
354
355 fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
356 if ( ! fEventCounters ) return;
357
358 if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(2);
359 fHistoEventsList->SetOwner();
360
361 TH2D* histo = 0x0;
362 histo = fEventCounters->Draw("trigger","selection","event:any");
363 if ( histo ) {
364 histo->SetName("hEvents");
365 histo->SetTitle("Events per trigger");
366 fHistoEventsList->AddAtAndExpand(histo, kHevents);
367 }
368
369 histo = fEventCounters->Draw("trigger","selection","event:correctedL0");
370 if ( histo ) {
371 histo->SetName("hEventsCorrectL0");
372 histo->SetTitle("L0 corrected events");
373 fHistoEventsList->AddAtAndExpand(histo, kHeventsCorrectL0);
374 }
375
376 PostData(2, fHistoEventsList);
377
378 if ( gROOT->IsBatch() )
379 return;
380
381 TH2D* histoPileupL0 = (TH2D*)fHistoEventsList->At(kHeventsCorrectL0)->Clone("hPileupL0");
382 histoPileupL0->Divide((TH2D*)fHistoEventsList->At(kHevents));
383
384 TCanvas *can = new TCanvas("can1_Pileup","Pileup",10,10,310,310);
385 can->SetFillColor(10); can->SetHighLightColor(10);
386 can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
387 histoPileupL0->DrawCopy("text");
388}
389
390
391//________________________________________________________________________
392Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB)
393{
394 /// Get the correction factor for L0 calculation
395
396 if ( nCBEAMB == 0. )
397 return 1.;
398
399 Double_t ratio = nCINT1B / nCBEAMB;
400
401 if ( ratio >= 1. || ratio == 0. )
402 return 1.;
403
404 Double_t mu = -TMath::Log(1-ratio);
405
406 return mu / ( 1. - TMath::Exp(-mu) );
407
408}
409