1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-------------------------------------------------------------------------
19 // Implementation of class AliAnalysisTaskPileup
21 // The class re-weights the number of analysed events with the correction
22 // factor to account for the event pileup.
24 // The final idea will be to include all of the methods for the pileup
27 // For the moment, the implemented method is based on a work
29 // The probability of having at least 1 collision is given by:
30 // MB L0 trigger / max collision rate
31 // The max collision rate is given by:
32 // number of bunches X revolution frequency ~ CBEAMB L0 trigger
33 // Assuming a poissonian distribution of the collision, with:
34 // mu = mean number of collisions
35 // the probability of having 0 collisions is:
36 // exp{-mu} = 1 - L0b_MB / max collision rate
37 // Since the MB trigger is measuring the probability of having at least
38 // 1 collisions, the number should be corrected by the factor:
39 // CF = mu / ( 1 - exp{-mu} )
40 // The class weights the number of events with this correction factor
43 // - so far the class needs to access the OCDB
44 // (hopefully it will be changed in the future)
45 // Hence the following libraries need to be included (in addition to
46 // the standard analysis ones):
47 // libRAWDatabase.so libCDB.so libSTEER.so libPWG3base.so
48 //-------------------------------------------------------------------------
50 #include <Riostream.h>
58 #include "TObjArray.h"
59 #include "TObjString.h"
63 #include "AliESDEvent.h"
64 #include "AliESDInputHandler.h"
65 #include "AliAODEvent.h"
66 #include "AliAODInputHandler.h"
67 #include "AliTriggerScalersESD.h"
68 #include "AliCentrality.h"
71 #include "AliAnalysisTaskSE.h"
72 #include "AliAnalysisManager.h"
73 #include "AliAnalysisDataSlot.h"
74 #include "AliAnalysisDataContainer.h"
76 #include "AliCounterCollection.h"
78 #include "AliAnalysisTaskPileup.h"
81 #include "AliTriggerRunScalers.h"
82 #include "AliTriggerConfiguration.h"
83 #include "AliTriggerClass.h"
84 #include "AliCDBManager.h"
85 #include "AliCDBEntry.h"
86 #include "AliCDBPath.h"
90 ClassImp(AliAnalysisTaskPileup)
93 //________________________________________________________________________
94 AliAnalysisTaskPileup::AliAnalysisTaskPileup() :
97 fHistoEventsList(0x0),
99 fTriggerClassIndex(0x0),
101 fCentralityClasses(0x0)
102 #if defined(READOCDB)
103 , fTriggerRunScalers(0x0),
108 /// Default Constructor
112 //________________________________________________________________________
113 AliAnalysisTaskPileup::AliAnalysisTaskPileup(const char *name) :
114 AliAnalysisTaskSE(name),
116 fHistoEventsList(0x0),
117 fTriggerClasses(0x0),
118 fTriggerClassIndex(0x0),
120 fCentralityClasses(0x0)
121 #if defined(READOCDB)
122 , fTriggerRunScalers(0x0),
129 DefineOutput(1,AliCounterCollection::Class());
130 DefineOutput(2,TObjArray::Class());
133 //________________________________________________________________________
134 AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
138 // For proof: do not delete output containers
139 if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
140 delete fEventCounters;
143 delete fHistoEventsList;
144 delete fTriggerClasses;
145 delete fTriggerClassIndex;
147 #if defined(READOCDB)
148 //delete fTriggerRunScalers; // Not owner -> Owned by OCDB
155 //___________________________________________________________________________
156 void AliAnalysisTaskPileup::NotifyRun()
160 #if defined(READOCDB)
161 fStorageList->Compress();
162 if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) {
163 for ( Int_t ientry=0; ientry<fStorageList->GetEntries(); ientry++ ) {
164 TObjString* calibStr = (TObjString*)fStorageList->At(ientry);
166 TObjString* dbStr = (TObjString*)fStorageList->At(ientry);
167 TString calibName = calibStr->GetString();
168 if ( ! calibName.CompareTo("default") ) {
169 AliCDBManager::Instance()->SetDefaultStorage(dbStr->GetName());
172 AliCDBManager::Instance()->SetSpecificStorage(calibStr->GetName(), dbStr->GetName());
177 // Default storage was not correclty set: nothing done
178 if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) return;
180 AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber());
182 AliCDBEntry *entry = 0x0;
183 entry = AliCDBManager::Instance()->Get("GRP/CTP/Config");
184 if ( ! entry ) return;
185 AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject();
186 const TObjArray& classesArray = trigConf->GetClasses();
187 Int_t nclasses = classesArray.GetEntriesFast();
189 if ( fTriggerClasses ) delete fTriggerClasses;
191 fTriggerClasses = new TObjArray(nclasses);
192 fTriggerClasses->SetOwner();
194 Int_t currActive = -1;
195 for( Int_t iclass=0; iclass < nclasses; iclass++ ) {
196 AliTriggerClass* trclass = (AliTriggerClass*)classesArray.At(iclass);
197 if (trclass && trclass->GetMask()>0) {
199 Int_t currPos = currActive;
201 // Store the CBEAMB class at the first position
202 TString trigName = trclass->GetName();
203 if ( ( trigName.Contains("CBEAMB") || trigName.Contains("CTRUE") ) && ! trigName.Contains("WU") )
205 else if ( ! fTriggerClasses->At(0) )
208 Int_t trindex = TMath::Nint(TMath::Log2(trclass->GetMask()));
209 TObjString* objString = new TObjString(trigName.Data());
210 fTriggerClasses->AddAtAndExpand(objString, currPos);
211 (*fTriggerClassIndex)[currPos] = trindex;
212 if ( fDebug >= 3 ) printf("AliAnalysisTaskPileup: Current class %s index %i position %i\n", trigName.Data(), trindex, currPos);
214 } // loop on trigger classes
216 entry = AliCDBManager::Instance()->Get("GRP/CTP/Scalers");
217 if ( ! entry ) return;
218 AliInfo("Found an AliTriggerRunScalers in GRP/CTP/Scalers, reading it");
219 fTriggerRunScalers = static_cast<AliTriggerRunScalers*> (entry->GetObject());
220 if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow");
229 //___________________________________________________________________________
230 void AliAnalysisTaskPileup::UserCreateOutputObjects()
232 /// Create histograms and counters
234 fTriggerClassIndex = new TArrayI(50);
235 fTriggerClassIndex->Reset(-1);
237 fCentralityClasses = new TAxis(20, 0., 100.);
239 TString centralityClassesStr = "", currClass = "";
240 for ( Int_t ibin=1; ibin<=fCentralityClasses->GetNbins(); ibin++ ){
241 if ( ! centralityClassesStr.IsNull() )
242 centralityClassesStr.Append("/");
243 currClass = Form("%.0f-%.0f",fCentralityClasses->GetBinLowEdge(ibin),fCentralityClasses->GetBinUpEdge(ibin));
244 centralityClassesStr += currClass;
245 fCentralityClasses->SetBinLabel(ibin, currClass.Data());
248 // The framework has problems if the name of the object
249 // and the one of container differ
250 // To be complaint, get the name from container and set it
251 TString containerName = GetOutputSlot(1)->GetContainer()->GetName();
253 // initialize event counters
254 fEventCounters = new AliCounterCollection(containerName.Data());
255 fEventCounters->AddRubric("event", "any/correctedL0");
256 fEventCounters->AddRubric("trigger", 1000000);
257 fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
258 fEventCounters->AddRubric("selected", "yes/no");
259 fEventCounters->AddRubric("centrality", centralityClassesStr.Data());
260 fEventCounters->Init(kTRUE);
262 // Post data at least once per task to ensure data synchronisation (required for merging)
263 PostData(1, fEventCounters);
266 //________________________________________________________________________
267 void AliAnalysisTaskPileup::UserExec(Option_t *)
269 /// Called for each event
271 AliAODEvent* aodEvent = 0x0;
273 AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
275 aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
278 if ( ! aodEvent && ! esdEvent ) {
279 AliError ("AOD or ESD event not found. Nothing done!");
283 // check physics selection
284 Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
285 TString selected = ( isPhysicsSelected ) ? "yes" : "no";
287 TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
289 if ( ! fIsInitCDB ) {
290 delete fTriggerClasses;
291 fTriggerClasses = firedTrigClasses.Tokenize(" ");
292 fTriggerClasses->SetOwner();
295 #if defined(READOCDB)
297 //if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events
299 if ( fDebug >= 2 ) printf("\nAliAnalysisTaskPileup: Event %lli\n", Entry());
301 // If scalers is not correctly loaded, the task becomes a mere event counter
302 Int_t nPoints = ( fTriggerRunScalers ) ? fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast() : 0;
304 AliTriggerScalersRecordESD* trigScalerRecords1 = 0x0;
305 AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
307 // Add protection for MC (no scalers there!)
309 AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber());
310 Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp);
311 if ( position < 0 ) {
312 AliWarning("Position out of range: put to 1");
315 if ( position == 0 ) position++; // Don't trust the first one
316 else if ( position + 1 >= nPoints ) position--;
317 if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: position %i\n", position);
318 trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);
320 // Sometimes scalers are filled very close to each others
321 // in this case skip to the next entry
322 for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) {
323 trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos);
324 Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() );
325 if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Pos %i TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime);
329 } // loop on position
334 ULong64_t trigMask = 0;
335 Double_t correctFactor = 1.;
337 Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
338 Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8);
340 Double_t centralityClass = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
341 // If percentile is exactly 100. the bin chosen is in overflow
342 // fix it by setting 99.999 instead of 100.
343 if ( centralityClass >= 100. ) centralityClass = 99.999;
344 Int_t centralityBin = fCentralityClasses->FindBin(centralityClass);
346 TString vtxSelKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
347 Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )};
349 //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT
351 TString trigName = "";
352 TString eventType = "";
354 // Get entries counts the number of entries != 0
355 // However the position 0 can be empty in MC (it is reserved to CBEAMB)
356 // In this case with GetEntries the last entry is lost
357 // Use GetEntriesFast instead
358 Int_t nTriggerClasses = fTriggerClasses->GetEntriesFast();
359 Int_t classIndex = -1;
360 #if defined(READOCDB)
361 Double_t deltaScalersBeam = 0., deltaScalers = 0.;
363 //Bool_t isFiredOnce = kFALSE;
364 for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) {
366 Double_t correctFactorL0 = 1.;
368 Bool_t isClassFired = kFALSE;
370 if ( itrig < nTriggerClasses ) {
373 // Check if current mask contains trigger
374 classIndex = (*fTriggerClassIndex)[itrig];
375 if ( classIndex < 0 ) continue; // Protection for MC (where BEAMB not present
376 trigMask = ( 1ull << classIndex );
377 isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
380 isClassFired = kTRUE;
382 trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
384 #if defined(READOCDB)
385 if ( trigScalerRecords2 && ( isClassFired || itrig == 0 ) ) {
387 const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1);
388 const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1);
389 deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB();
392 deltaScalersBeam = deltaScalers;
393 else if ( isClassFired ) {
394 correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam);
395 if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Scalers: %s %.0f %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0);
399 } // if ( itrig < nTriggerClasses )
403 isClassFired = kTRUE; // isFiredOnce;
406 if ( ! isClassFired ) continue;
407 //isFiredOnce = kTRUE;
409 //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT
410 //if ( classIndex > 1 ) printf("Index: trigger %i scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT
412 if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Fired trig %s\n", trigName.Data());
414 for ( Int_t ievType=0; ievType<2; ievType++ ){
416 case kHeventsCorrectL0:
417 correctFactor = correctFactorL0;
418 eventType = "correctedL0";
426 for ( Int_t isel=0; isel<3; isel++ ) {
427 if ( ! fillSel[isel] ) continue;
428 fEventCounters->Count(Form("event:%s/trigger:%s/vtxSelection:%s/selected:%s/centrality:%s",eventType.Data(),trigName.Data(), vtxSelKey[isel].Data(), selected.Data(),fCentralityClasses->GetBinLabel(centralityBin)),correctFactor);
429 } // loop on vertex selection
430 } // loop on event type
431 } // loop on trigger classes
433 // Post final data. It will be written to a file with option "RECREATE"
434 PostData(1, fEventCounters);
438 //________________________________________________________________________
439 void AliAnalysisTaskPileup::Terminate(Option_t *)
442 /// Save final histograms
443 /// and draw result to the screen
446 // Fill the container only at the very last step
447 // i.e. in local, when done interactively
448 if ( gROOT->IsBatch() )
451 fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
452 if ( ! fEventCounters ) return;
454 if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(0);
455 fHistoEventsList->SetOwner();
458 const Int_t kNevtTypes = 2;
459 TString evtSel[kNevtTypes] = {"event:any", "event:correctedL0"};
460 TString evtName[kNevtTypes] = {"", "CorrectL0"};
461 TString evtTitle[kNevtTypes] = {"Events", "L0 corrected events"};
462 const Int_t kNphysSel = 2;
463 TString physSel[kNphysSel] = {"selected:any","selected:yes"};
464 TString physName[kNphysSel] = {"", "PhysSel"};
465 TString physTitle[kNphysSel] = {"", "w/ physics selection"};
467 TString typeName[2] = {"", "Centrality"};
469 Int_t currHisto = -1;
470 TString currName = "";
471 for ( Int_t itype=0; itype<2; itype++ ) {
472 for ( Int_t isel=0; isel<kNphysSel; isel++ ) {
473 for ( Int_t iev=0; iev<kNevtTypes; iev++ ) {
474 currName = Form("%s/%s", evtSel[iev].Data(), physSel[isel].Data());
476 histo = fEventCounters->Get("trigger","vtxSelection",currName.Data());
478 currName.Append("/vtxSelection:any");
479 histo = fEventCounters->Get("trigger","centrality",currName.Data());
481 if ( ! histo ) continue;
483 currName = Form("hEvents%s%s%s", typeName[itype].Data(), evtName[iev].Data(), physName[isel].Data());
484 histo->SetName(currName.Data());
485 currName = Form("%s %s", evtTitle[iev].Data(), physTitle[isel].Data());
486 histo->SetTitle(currName.Data());
487 fHistoEventsList->AddAtAndExpand(histo, currHisto);
488 } // loop on event type
489 TH2D* num = (TH2D*)fHistoEventsList->At(currHisto);
490 TH2D* den = (TH2D*)fHistoEventsList->At(currHisto-1);
491 if ( ! num || ! den ) continue;
492 currName = Form("hPileupL0%s%s_%s", typeName[itype].Data(), physName[isel].Data(), GetName());
493 TH2D* histoPileupL0 = (TH2D*)num->Clone(currName.Data());
494 histoPileupL0->Divide(den);
495 currName.ReplaceAll("hPileupL0","canPileupL0");
496 TCanvas *can = new TCanvas(currName.Data(),"Pileup",10,10,310,310);
497 can->SetFillColor(10); can->SetHighLightColor(10);
498 can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
499 histoPileupL0->DrawCopy("text");
500 delete histoPileupL0;
501 } // loop on physics selection
502 } // loop on histo type
504 PostData(2, fHistoEventsList);
508 //________________________________________________________________________
509 Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB)
511 /// Get the correction factor for L0 calculation
516 Double_t ratio = nCINT1B / nCBEAMB;
518 if ( ratio >= 1. || ratio == 0. )
521 Double_t mu = -TMath::Log(1-ratio);
523 return mu / ( 1. - TMath::Exp(-mu) );
527 //________________________________________________________________________
528 void AliAnalysisTaskPileup::SetDefaultStorage(TString dbString)
530 /// Set default storage
531 SetSpecificStorage("default", dbString);
535 //________________________________________________________________________
536 void AliAnalysisTaskPileup::SetSpecificStorage(TString calibType, TString dbString)
538 /// Set specific storage
539 #if defined(READOCDB)
540 if ( ! fStorageList ) {
541 fStorageList = new TObjArray(5);
542 fStorageList->SetOwner();
544 fStorageList->AddLast(new TObjString(calibType));
545 fStorageList->AddLast(new TObjString(dbString));
549 AliWarning(Form("Class was not compiled to run on OCDB. Command will not have effect"));