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 //-------------------------------------------------------------------------
  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
  88 ClassImp(AliAnalysisTaskPileup)
  89
  90
  91 //________________________________________________________________________
  92 AliAnalysisTaskPileup::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   fTriggerClassIndex->Reset(-1);
 107
 108   DefineOutput(1,AliCounterCollection::Class());
 109   DefineOutput(2,TObjArray::Class());
 110 }
 111
 112 //________________________________________________________________________
 113 AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
 114 {
 115   /// Destructor
 116
 117   // For proof: do not delete output containers
 118   if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
 119     delete fEventCounters;
 120   }
 121
 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 //___________________________________________________________________________
 135 void 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();
 167         if ( ( trigName.Contains("CBEAMB") ||  trigName.Contains("CTRUE") ) && ! trigName.Contains("WU") )
 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 //___________________________________________________________________________
 194 void 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);
 207   fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
 208   fEventCounters->AddRubric("selected", "yes/no");
 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 //________________________________________________________________________
 216 void 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
 233   Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
 234   TString selected = ( isPhysicsSelected ) ? "yes" : "no";
 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
 246   // Add protection for MC (no scalers there!)
 247   AliTriggerScalersRecordESD* trigScalerRecords1 = 0x0;
 248   AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
 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
 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
 282   TString vtxSelKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
 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
 303       classIndex = (*fTriggerClassIndex)[itrig];
 304       if ( classIndex < 0 ) continue; // Protection for MC (where BEAMB not present
 305       trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
 306       trigMask = ( 1ull << classIndex );
 307       isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
 308
 309 #if defined(READOCDB)
 310       if ( trigScalerRecords2 && ( isClassFired || itrig == 0 ) ) {
 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();
 315
 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       }
 323 #endif
 324     } // if ( itrig < nTriggerClasses )
 325     else {
 326       classIndex = -1;
 327       trigName = "any";
 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;
 354         fEventCounters->Count(Form("event:%s/trigger:%s/vtxSelection:%s/selected:%s",eventType.Data(),trigName.Data(), vtxSelKey[isel].Data(), selected.Data()),correctFactor);
 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 //________________________________________________________________________
 365 void AliAnalysisTaskPileup::Terminate(Option_t *)
 366 {
 367   //
 368   /// Save final histograms
 369   /// and draw result to the screen
 370   //
 371
 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
 377   fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
 378   if ( ! fEventCounters ) return;
 379
 380   if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(0);
 381   fHistoEventsList->SetOwner();
 382
 383   TH2D* histo = 0x0;
 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
 420
 421   PostData(2, fHistoEventsList);
 422 }
 423
 424
 425 //________________________________________________________________________
 426 Double_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