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   DefineOutput(1,AliCounterCollection::Class());
 107   DefineOutput(2,TObjArray::Class());
 108 }
 109
 110 //________________________________________________________________________
 111 AliAnalysisTaskPileup::~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 //___________________________________________________________________________
 128 void 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 //___________________________________________________________________________
 187 void 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 //________________________________________________________________________
 208 void 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 //________________________________________________________________________
 348 void 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 //________________________________________________________________________
 392 Double_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