Fix to correctly set mother type with Phojet simulations. Account for centrality...
[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
c705eb52 106 fTriggerClassIndex->Reset(-1);
107
0fa651e5 108 DefineOutput(1,AliCounterCollection::Class());
109 DefineOutput(2,TObjArray::Class());
110}
111
112//________________________________________________________________________
113AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
114{
115 /// Destructor
c705eb52 116
117 // For proof: do not delete output containers
118 if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
119 delete fEventCounters;
120 }
121
0fa651e5 122 delete fHistoEventsList;
123 delete fTriggerClasses;
124 delete fTriggerClassIndex;
125
126#if defined(READOCDB)
127 delete fTriggerRunScalers;
128 delete fDefaultStorage;
129#endif
130
131}
132
133
134//___________________________________________________________________________
135void AliAnalysisTaskPileup::NotifyRun()
136{
137 /// Notify run
138
139#if defined(READOCDB)
140 if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) {
141 AliCDBManager::Instance()->SetDefaultStorage(fDefaultStorage->Data());
142 }
143
144 AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber());
145
146 AliCDBEntry *entry = 0x0;
147 entry = AliCDBManager::Instance()->Get("GRP/CTP/Config");
148 if ( entry ) {
149 AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject();
150 const TObjArray& classesArray = trigConf->GetClasses();
151 Int_t nclasses = classesArray.GetEntriesFast();
152
153 if ( fTriggerClasses ) delete fTriggerClasses;
154
155 fTriggerClasses = new TObjArray(nclasses);
156 fTriggerClasses->SetOwner();
157
158 Int_t currActive = -1;
159 for( Int_t iclass=0; iclass < nclasses; iclass++ ) {
160 AliTriggerClass* trclass = (AliTriggerClass*)classesArray.At(iclass);
161 if (trclass && trclass->GetMask()>0) {
162 currActive++;
163 Int_t currPos = currActive;
164
165 // Store the CBEAMB class at the first position
166 TString trigName = trclass->GetName();
c705eb52 167 if ( ( trigName.Contains("CBEAMB") || trigName.Contains("CTRUE") ) && ! trigName.Contains("WU") )
0fa651e5 168 currPos = 0;
169 else if ( ! fTriggerClasses->At(0) )
170 currPos++;
171
172 Int_t trindex = TMath::Nint(TMath::Log2(trclass->GetMask()));
173 TObjString* objString = new TObjString(trigName.Data());
174 fTriggerClasses->AddAtAndExpand(objString, currPos);
175 (*fTriggerClassIndex)[currPos] = trindex;
176 AliDebug(3, Form("Current class %s index %i position %i", trigName.Data(), trindex, currPos));
177 } // is class active
178 } // loop on trigger classes
179 } // if entry
180
181 entry = AliCDBManager::Instance()->Get("GRP/CTP/Scalers");
182 if (entry) {
183 AliInfo("Found an AliTriggerRunScalers in GRP/CTP/Scalers, reading it");
184 fTriggerRunScalers = dynamic_cast<AliTriggerRunScalers*> (entry->GetObject());
185 entry->SetOwner(0);
186 if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow");
187 }
188#endif
189
190}
191
192
193//___________________________________________________________________________
194void AliAnalysisTaskPileup::UserCreateOutputObjects()
195{
196 /// Create histograms and counters
197
198 // The framework has problems if the name of the object
199 // and the one of container differ
200 // To be complaint, get the name from container and set it
201 TString containerName = GetOutputSlot(1)->GetContainer()->GetName();
202
203 // initialize event counters
204 fEventCounters = new AliCounterCollection(containerName.Data());
205 fEventCounters->AddRubric("event", "any/correctedL0");
206 fEventCounters->AddRubric("trigger", 1000000);
c705eb52 207 fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
208 fEventCounters->AddRubric("selected", "yes/no");
0fa651e5 209 fEventCounters->Init(kTRUE);
210
211 // Post data at least once per task to ensure data synchronisation (required for merging)
212 PostData(1, fEventCounters);
213}
214
215//________________________________________________________________________
216void AliAnalysisTaskPileup::UserExec(Option_t *)
217{
218 /// Called for each event
219
220 AliAODEvent* aodEvent = 0x0;
221
222 AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
223 if ( ! esdEvent ) {
224 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
225 }
226
227 if ( ! aodEvent && ! esdEvent ) {
228 AliError ("AOD or ESD event not found. Nothing done!");
229 return;
230 }
231
232 // check physics selection
c705eb52 233 Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
234 TString selected = ( isPhysicsSelected ) ? "yes" : "no";
0fa651e5 235
236#if defined(READOCDB)
237
238 TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
239
240 if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events
241
242 AliDebug(2, Form("\nEvent %lli\n", Entry()));
243
244 Int_t nPoints = fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast();
245
c705eb52 246 // Add protection for MC (no scalers there!)
247 AliTriggerScalersRecordESD* trigScalerRecords1 = 0x0;
0fa651e5 248 AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
c705eb52 249 if ( nPoints > 1 ) {
250
251 AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber());
252 Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp);
253 if ( position < 0 ) {
254 AliWarning("Position out of range: put to 1");
255 position = 1;
256 }
257 if ( position == 0 ) position++; // Don't trust the first one
258 else if ( position + 1 >= nPoints ) position--;
259 AliDebug(2, Form("position %i\n", position));
260 trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);
261
262 // Sometimes scalers are filled very close to each others
263 // in this case skip to the next entry
264 for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) {
265 trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos);
266 Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() );
267 AliDebug(2, Form("Pos %i TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime));
268
269 if ( deltaTime > 1 )
270 break;
271 } // loop on position
272 } // nPoins > 0
0fa651e5 273
274#endif
275
276 ULong64_t trigMask = 0;
277 Double_t correctFactor = 1.;
278
279 Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
280 Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8);
281
c705eb52 282 TString vtxSelKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
0fa651e5 283 Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )};
284
285 //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT
286
287 TString trigName = "";
288 TString eventType = "";
289
290 Int_t nTriggerClasses = fTriggerClasses->GetEntries();
291 Int_t classIndex = -1;
292 Double_t deltaScalersBeam = 0., deltaScalers = 0.;
293 Bool_t isFiredOnce = kFALSE;
294 for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) {
295
296 Double_t correctFactorL0 = 1.;
297
298 Bool_t isClassFired = kFALSE;
299
300 if ( itrig < nTriggerClasses ) {
301
302 // Check if current mask contains trigger
0fa651e5 303 classIndex = (*fTriggerClassIndex)[itrig];
c705eb52 304 if ( classIndex < 0 ) continue; // Protection for MC (where BEAMB not present
305 trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
0fa651e5 306 trigMask = ( 1ull << classIndex );
307 isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
308
c705eb52 309#if defined(READOCDB)
310 if ( trigScalerRecords2 && ( isClassFired || itrig == 0 ) ) {
0fa651e5 311 // Get scalers
312 const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1);
313 const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1);
314 deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB();
c705eb52 315
0fa651e5 316 if ( itrig == 0 )
317 deltaScalersBeam = deltaScalers;
318 else if ( isClassFired ) {
319 correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam);
320 AliDebug(2, Form("Scalers: %s %.0f %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0));
321 }
322 }
c705eb52 323#endif
324 } // if ( itrig < nTriggerClasses )
0fa651e5 325 else {
0fa651e5 326 classIndex = -1;
c705eb52 327 trigName = "any";
0fa651e5 328 isClassFired = isFiredOnce;
329 }
330
331 if ( ! isClassFired ) continue;
332 isFiredOnce = kTRUE;
333
334 //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT
335 //if ( classIndex > 1 ) printf("Index: trigger %i scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT
336
337 AliDebug(2, Form("Fired trig %s\n", trigName.Data()));
338
339 for ( Int_t ievType=0; ievType<2; ievType++ ){
340 switch ( ievType ) {
341#if defined(READOCDB)
342 case kHeventsCorrectL0:
343 correctFactor = correctFactorL0;
344 eventType = "correctedL0";
345 break;
346#endif
347 default:
348 correctFactor = 1.;
349 eventType = "any";
350 }
351
352 for ( Int_t isel=0; isel<3; isel++ ) {
353 if ( ! fillSel[isel] ) continue;
c705eb52 354 fEventCounters->Count(Form("event:%s/trigger:%s/vtxSelection:%s/selected:%s",eventType.Data(),trigName.Data(), vtxSelKey[isel].Data(), selected.Data()),correctFactor);
0fa651e5 355 } // loop on vertex selection
356 } // loop on event type
357 } // loop on trigger classes
358
359 // Post final data. It will be written to a file with option "RECREATE"
360 PostData(1, fEventCounters);
361}
362
363
364//________________________________________________________________________
365void AliAnalysisTaskPileup::Terminate(Option_t *)
366{
367 //
368 /// Save final histograms
369 /// and draw result to the screen
370 //
371
c705eb52 372 // Fill the container only at the very last step
373 // i.e. in local, when done interactively
374 if ( gROOT->IsBatch() )
375 return;
376
0fa651e5 377 fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
378 if ( ! fEventCounters ) return;
379
c705eb52 380 if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(0);
0fa651e5 381 fHistoEventsList->SetOwner();
382
383 TH2D* histo = 0x0;
c705eb52 384 const Int_t kNevtTypes = 2;
385 TString evtSel[kNevtTypes] = {"event:any", "event:correctedL0"};
386 TString evtName[kNevtTypes] = {"", "CorrectL0"};
387 TString evtTitle[kNevtTypes] = {"Events", "L0 corrected events"};
388 const Int_t kNphysSel = 2;
389 TString physSel[kNphysSel] = {"selected:any","selected:yes"};
390 TString physName[kNphysSel] = {"", "PhysSel"};
391 TString physTitle[kNphysSel] = {"", "w/ physics selection"};
392
393 Int_t currHisto = -1;
394 TString currName = "";
395 for ( Int_t isel=0; isel<kNphysSel; isel++ ) {
396 for ( Int_t iev=0; iev<kNevtTypes; iev++ ) {
397 currName = Form("%s/%s", evtSel[iev].Data(), physSel[isel].Data());
398 histo = fEventCounters->Get("trigger","vtxSelection",currName.Data());
399 if ( ! histo ) continue;
400 currHisto++;
401 currName = Form("hEvents%s%s", evtName[iev].Data(), physName[isel].Data());
402 histo->SetName(currName.Data());
403 currName = Form("%s %s", evtTitle[iev].Data(), physTitle[isel].Data());
404 histo->SetTitle(currName.Data());
405 fHistoEventsList->AddAtAndExpand(histo, currHisto);
406 } // loop on event type
407 TH2D* num = (TH2D*)fHistoEventsList->At(currHisto);
408 TH2D* den = (TH2D*)fHistoEventsList->At(currHisto-1);
409 if ( ! num || ! den ) continue;
410 currName = Form("hPileupL0%s_%s", physName[isel].Data(), GetName());
411 TH2D* histoPileupL0 = (TH2D*)num->Clone(currName.Data());
412 histoPileupL0->Divide(den);
413 currName = Form("c%i_%s", isel, GetName());
414 TCanvas *can = new TCanvas(currName.Data(),"Pileup",10,10,310,310);
415 can->SetFillColor(10); can->SetHighLightColor(10);
416 can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
417 histoPileupL0->DrawCopy("text");
418 delete histoPileupL0;
419 } // loop on physics selection
0fa651e5 420
421 PostData(2, fHistoEventsList);
0fa651e5 422}
423
424
425//________________________________________________________________________
426Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB)
427{
428 /// Get the correction factor for L0 calculation
429
430 if ( nCBEAMB == 0. )
431 return 1.;
432
433 Double_t ratio = nCINT1B / nCBEAMB;
434
435 if ( ratio >= 1. || ratio == 0. )
436 return 1.;
437
438 Double_t mu = -TMath::Log(1-ratio);
439
440 return mu / ( 1. - TMath::Exp(-mu) );
441
442}
443