2d3f91d1099eee94bb05e19a6dfc7188547ce3a0
[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskPileup.cxx
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 /* $Id$ */
17
18 //-------------------------------------------------------------------------
19 //             Implementation of class AliAnalysisTaskPileup
20 //
21 // The class re-weights the number of analysed events with the correction
22 // factor to account for the event pileup.
23 //
24 // The final idea will be to include all of the methods for the pileup
25 // calculation.
26 //
27 // For the moment, the implemented method is based on a work 
28 // by L. Aphecetche.
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
41 //
42 // CAVEATS:
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 //-------------------------------------------------------------------------
49
50 #include <Riostream.h>
51
52 // ROOT includes
53 #include "TH1.h"
54 #include "TH2.h"
55 #include "TCanvas.h"
56 #include "TROOT.h"
57 #include "TString.h"
58 #include "TObjArray.h"
59 #include "TObjString.h"
60 #include "TMath.h"
61
62 // STEER includes
63 #include "AliESDEvent.h"
64 #include "AliESDInputHandler.h"
65 #include "AliAODEvent.h"
66 #include "AliAODInputHandler.h"
67 #include "AliTriggerScalersESD.h"
68 #include "AliCentrality.h"
69
70 // ANALYSIS includes
71 #include "AliAnalysisTaskSE.h"
72 #include "AliAnalysisManager.h"
73 #include "AliAnalysisDataSlot.h"
74 #include "AliAnalysisDataContainer.h"
75
76 #include "AliCounterCollection.h"
77
78 #include "AliAnalysisTaskPileup.h"
79
80 #if defined(READOCDB)
81 #include "AliTriggerRunScalers.h"
82 #include "AliTriggerConfiguration.h"
83 #include "AliTriggerClass.h"
84 #include "AliCDBManager.h"
85 #include "AliCDBEntry.h"
86 #include "AliCDBPath.h"
87 #endif
88
89
90 ClassImp(AliAnalysisTaskPileup)
91
92
93 //________________________________________________________________________
94 AliAnalysisTaskPileup::AliAnalysisTaskPileup() :
95 AliAnalysisTaskSE(),
96   fEventCounters(0x0),
97   fHistoEventsList(0x0),
98   fTriggerClasses(0x0),
99   fTriggerClassIndex(0x0),
100   fIsInitCDB(0),
101   fCentralityClasses(0x0)
102 #if defined(READOCDB)
103   , fTriggerRunScalers(0x0),
104   fStorageList(0x0)
105 #endif
106
107 {
108   /// Default Constructor
109 }
110
111
112 //________________________________________________________________________
113 AliAnalysisTaskPileup::AliAnalysisTaskPileup(const char *name) :
114   AliAnalysisTaskSE(name),
115   fEventCounters(0x0),
116   fHistoEventsList(0x0),
117   fTriggerClasses(0x0),
118   fTriggerClassIndex(0x0),
119   fIsInitCDB(0),
120   fCentralityClasses(0x0)
121 #if defined(READOCDB)
122   , fTriggerRunScalers(0x0),
123   fStorageList(0x0)
124 #endif
125
126 {
127   /// Constructor
128
129   DefineOutput(1,AliCounterCollection::Class());
130   DefineOutput(2,TObjArray::Class());
131 }
132
133 //________________________________________________________________________
134 AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
135 {
136   /// Destructor
137
138   // For proof: do not delete output containers
139   if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
140     delete fEventCounters;
141   }
142
143   delete fHistoEventsList;
144   delete fTriggerClasses;
145   delete fTriggerClassIndex;
146
147 #if defined(READOCDB)
148   delete fTriggerRunScalers;
149   delete fStorageList;
150 #endif
151
152 }
153
154
155 //___________________________________________________________________________
156 void AliAnalysisTaskPileup::NotifyRun()
157 {
158   /// Notify run
159
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);
165       ientry++;
166       TObjString* dbStr = (TObjString*)fStorageList->At(ientry);
167       TString calibName = calibStr->GetString();
168       if ( ! calibName.CompareTo("default") ) {
169         AliCDBManager::Instance()->SetDefaultStorage(dbStr->GetName());
170       }
171       else {
172         AliCDBManager::Instance()->SetSpecificStorage(calibStr->GetName(), dbStr->GetName());
173       }
174     }
175   }
176
177   // Default storage was not correclty set: nothing done
178   if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) return;
179
180   AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber());
181
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();
188
189   if ( fTriggerClasses ) delete fTriggerClasses;
190
191   fTriggerClasses = new TObjArray(nclasses);
192   fTriggerClasses->SetOwner();
193
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) {
198       currActive++;
199       Int_t currPos = currActive;
200
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") )
204         currPos = 0;
205       else if ( ! fTriggerClasses->At(0) )
206         currPos++;
207
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);
213     } // is class active
214   } // loop on trigger classes
215
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 = dynamic_cast<AliTriggerRunScalers*> (entry->GetObject());
220   entry->SetOwner(0);
221   if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow");
222
223   fIsInitCDB = kTRUE;
224
225 #endif
226
227 }
228
229
230 //___________________________________________________________________________
231 void AliAnalysisTaskPileup::UserCreateOutputObjects()
232 {
233   /// Create histograms and counters
234
235   fTriggerClassIndex = new TArrayI(50);
236   fTriggerClassIndex->Reset(-1);
237
238   fCentralityClasses = new TAxis(20, 0., 100.);
239
240   TString centralityClassesStr = "", currClass = "";
241   for ( Int_t ibin=1; ibin<=fCentralityClasses->GetNbins(); ibin++ ){
242     if ( ! centralityClassesStr.IsNull() )
243       centralityClassesStr.Append("/");
244     currClass = Form("%.0f-%.0f",fCentralityClasses->GetBinLowEdge(ibin),fCentralityClasses->GetBinUpEdge(ibin));
245     centralityClassesStr += currClass;
246     fCentralityClasses->SetBinLabel(ibin, currClass.Data());
247   }
248   
249   // The framework has problems if the name of the object
250   // and the one of container differ
251   // To be complaint, get the name from container and set it
252   TString containerName = GetOutputSlot(1)->GetContainer()->GetName();
253
254   // initialize event counters
255   fEventCounters = new AliCounterCollection(containerName.Data());
256   fEventCounters->AddRubric("event", "any/correctedL0");
257   fEventCounters->AddRubric("trigger", 1000000);
258   fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
259   fEventCounters->AddRubric("selected", "yes/no");
260   fEventCounters->AddRubric("centrality", centralityClassesStr.Data());
261   fEventCounters->Init(kTRUE);
262
263   // Post data at least once per task to ensure data synchronisation (required for merging)
264   PostData(1, fEventCounters);
265 }
266
267 //________________________________________________________________________
268 void AliAnalysisTaskPileup::UserExec(Option_t *)
269 {
270   /// Called for each event
271
272   AliAODEvent* aodEvent = 0x0;
273
274   AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
275   if ( ! esdEvent ) {
276     aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
277   }
278
279   if ( ! aodEvent && ! esdEvent ) {
280     AliError ("AOD or ESD event not found. Nothing done!");
281     return;
282   }
283  
284   // check physics selection
285   Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
286   TString selected = ( isPhysicsSelected ) ? "yes" : "no";
287
288   TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();
289
290   if ( ! fIsInitCDB ) {
291     delete fTriggerClasses;
292     fTriggerClasses = firedTrigClasses.Tokenize(" ");
293     fTriggerClasses->SetOwner();
294   }
295
296 #if defined(READOCDB)
297
298   //if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events
299
300   if ( fDebug >= 2 ) printf("\nAliAnalysisTaskPileup: Event %lli\n", Entry());
301
302   // If scalers is not correctly loaded, the task becomes a mere event counter
303   Int_t nPoints = ( fTriggerRunScalers ) ? fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast() : 0;
304
305   AliTriggerScalersRecordESD* trigScalerRecords1 = 0x0;
306   AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
307   if ( nPoints > 1 ) {
308     // Add protection for MC (no scalers there!)
309
310     AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber());
311     Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp);
312     if ( position < 0 ) {
313       AliWarning("Position out of range: put to 1");
314       position = 1;
315     } 
316     if ( position == 0 ) position++; // Don't trust the first one
317     else if ( position + 1 >= nPoints ) position--;
318     if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: position %i\n", position);
319     trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);
320
321     // Sometimes scalers are filled very close to each others
322     // in this case skip to the next entry
323     for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) {
324       trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos);
325       Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() );
326       if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Pos %i  TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime);
327
328       if ( deltaTime > 1 )
329         break;
330     } // loop on position
331   } // nPoins > 0
332
333 #endif
334
335   ULong64_t trigMask = 0;
336   Double_t correctFactor = 1.;
337
338   Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
339   Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8);
340
341   Double_t centralityClass = InputEvent()->GetCentrality()->GetCentralityPercentile("V0M");
342   // If percentile is exactly 100. the bin chosen is in overflow
343   // fix it by setting 99.999 instead of 100.
344   if ( centralityClass >= 100. ) centralityClass = 99.999;
345   Int_t centralityBin = fCentralityClasses->FindBin(centralityClass);
346
347   TString vtxSelKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
348   Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )};
349
350   //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT
351
352   TString trigName = "";
353   TString eventType = "";
354
355   // Get entries counts the number of entries != 0
356   // However the position 0 can be empty in MC (it is reserved to CBEAMB)
357   // In this case with GetEntries the last entry is lost
358   // Use GetEntriesFast instead
359   Int_t nTriggerClasses = fTriggerClasses->GetEntriesFast();
360   Int_t classIndex = -1;
361 #if defined(READOCDB)
362   Double_t deltaScalersBeam = 0., deltaScalers = 0.;
363 #endif
364   Bool_t isFiredOnce = kFALSE;
365   for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) {
366
367     Double_t correctFactorL0 = 1.;
368
369     Bool_t isClassFired = kFALSE;
370
371     if ( itrig < nTriggerClasses ) {
372
373       if ( fIsInitCDB ) {
374         // Check if current mask contains trigger
375         classIndex = (*fTriggerClassIndex)[itrig];
376         if ( classIndex < 0 ) continue; // Protection for MC (where BEAMB not present
377         trigMask = ( 1ull << classIndex );
378         isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
379       }
380       else
381         isClassFired = kTRUE;
382
383       trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
384
385 #if defined(READOCDB)
386       if ( trigScalerRecords2 && ( isClassFired || itrig == 0 ) ) {
387         // Get scalers
388         const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1);
389         const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1);
390         deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB();
391  
392         if ( itrig == 0 )
393           deltaScalersBeam = deltaScalers;
394         else if ( isClassFired ) {
395           correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam);
396           if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Scalers: %s %.0f  %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0);
397         }
398       }
399 #endif
400     } // if ( itrig < nTriggerClasses )
401     else {
402       classIndex = -1;
403       trigName = "any";
404       isClassFired = isFiredOnce;
405     }
406
407     if ( ! isClassFired ) continue;
408     isFiredOnce = kTRUE;
409
410     //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT
411     //if ( classIndex > 1 ) printf("Index: trigger %i  scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT
412
413     if ( fDebug >= 2 ) printf("AliAnalysisTaskPileup: Fired trig %s\n", trigName.Data());
414
415     for ( Int_t ievType=0; ievType<2; ievType++ ){
416       switch ( ievType ) {
417       case kHeventsCorrectL0:
418         correctFactor = correctFactorL0;
419         eventType = "correctedL0";
420         break;
421       default:
422         correctFactor = 1.;
423         eventType = "any";
424         break;
425       }
426
427       for ( Int_t isel=0; isel<3; isel++ ) {
428         if ( ! fillSel[isel] ) continue;
429         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);
430       } // loop on vertex selection
431     } // loop on event type
432   } // loop on trigger classes
433
434   // Post final data. It will be written to a file with option "RECREATE"
435   PostData(1, fEventCounters);
436 }
437
438
439 //________________________________________________________________________
440 void AliAnalysisTaskPileup::Terminate(Option_t *)
441 {
442   //
443   /// Save final histograms
444   /// and draw result to the screen
445   //
446
447   // Fill the container only at the very last step
448   // i.e. in local, when done interactively
449   if ( gROOT->IsBatch() )
450     return;
451
452   fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
453   if ( ! fEventCounters ) return;
454
455   if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(0);
456   fHistoEventsList->SetOwner();
457
458   TH2D* histo = 0x0;
459   const Int_t kNevtTypes = 2;
460   TString evtSel[kNevtTypes] = {"event:any", "event:correctedL0"};
461   TString evtName[kNevtTypes] = {"", "CorrectL0"};
462   TString evtTitle[kNevtTypes] = {"Events", "L0 corrected events"};
463   const Int_t kNphysSel = 2;
464   TString physSel[kNphysSel] = {"selected:any","selected:yes"};
465   TString physName[kNphysSel] = {"", "PhysSel"};
466   TString physTitle[kNphysSel] = {"", "w/ physics selection"};
467
468   TString typeName[2] = {"", "Centrality"};
469
470   Int_t currHisto = -1;
471   TString currName = "";
472   for ( Int_t itype=0; itype<2; itype++ ) {
473     for ( Int_t isel=0; isel<kNphysSel; isel++ ) {
474       for ( Int_t iev=0; iev<kNevtTypes; iev++ ) {
475         currName = Form("%s/%s", evtSel[iev].Data(), physSel[isel].Data());
476         if ( itype == 0 )
477           histo = fEventCounters->Get("trigger","vtxSelection",currName.Data());
478         else {
479           currName.Append("/vtxSelection:any");
480           histo = fEventCounters->Get("trigger","centrality",currName.Data());
481         }
482         if ( ! histo ) continue;
483         currHisto++;
484         currName = Form("hEvents%s%s%s", typeName[itype].Data(), evtName[iev].Data(), physName[isel].Data());
485         histo->SetName(currName.Data());
486         currName = Form("%s %s", evtTitle[iev].Data(), physTitle[isel].Data());
487         histo->SetTitle(currName.Data());
488         fHistoEventsList->AddAtAndExpand(histo, currHisto);
489       } // loop on event type
490       TH2D* num = (TH2D*)fHistoEventsList->At(currHisto);
491       TH2D* den = (TH2D*)fHistoEventsList->At(currHisto-1);
492       if ( ! num || ! den ) continue;
493       currName = Form("hPileupL0%s%s_%s", typeName[itype].Data(), physName[isel].Data(), GetName());
494       TH2D* histoPileupL0 = (TH2D*)num->Clone(currName.Data());
495       histoPileupL0->Divide(den);
496       currName.ReplaceAll("hPileupL0","canPileupL0");
497       TCanvas *can = new TCanvas(currName.Data(),"Pileup",10,10,310,310);
498       can->SetFillColor(10); can->SetHighLightColor(10);
499       can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
500       histoPileupL0->DrawCopy("text");
501       delete histoPileupL0;
502     } // loop on physics selection
503   } // loop on histo type
504
505   PostData(2, fHistoEventsList);
506 }
507
508
509 //________________________________________________________________________
510 Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB)
511 {
512   /// Get the correction factor for L0 calculation
513
514   if ( nCBEAMB == 0. )
515     return 1.;
516
517   Double_t ratio = nCINT1B / nCBEAMB;
518
519   if ( ratio >= 1. || ratio == 0. )
520     return 1.;
521
522   Double_t mu = -TMath::Log(1-ratio);
523
524   return mu / ( 1. - TMath::Exp(-mu) );
525
526 }
527
528 //________________________________________________________________________
529 void AliAnalysisTaskPileup::SetDefaultStorage(TString dbString)
530 {
531   /// Set default storage
532   SetSpecificStorage("default", dbString);
533 }
534
535
536 //________________________________________________________________________
537 void AliAnalysisTaskPileup::SetSpecificStorage(TString calibType, TString dbString)
538 {
539   /// Set specific storage
540 #if defined(READOCDB)
541   if ( ! fStorageList ) {
542     fStorageList = new TObjArray(5);
543     fStorageList->SetOwner();
544   }
545   fStorageList->AddLast(new TObjString(calibType));
546   fStorageList->AddLast(new TObjString(dbString));
547 #else
548   calibType = "";
549   dbString  = "";
550   AliWarning(Form("Class was not compiled to run on OCDB. Command will not have effect"));
551 #endif
552 }
553