[u/mrichter/AliRoot.git] / PWG3 / muondep / AliAnalysisTaskPileup.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

//-------------------------------------------------------------------------
//             Implementation of class AliAnalysisTaskPileup
//
// The class re-weights the number of analysed events with the correction
// factor to account for the event pileup.
//
// The final idea will be to include all of the methods for the pileup
// calculation.
//
// For the moment, the implemented method is based on a work 
// by L. Aphecetche.
// The probability of having at least 1 collision is given by:
// MB L0 trigger / max collision rate
// The max collision rate is given by:
//    number of bunches X revolution frequency ~ CBEAMB L0 trigger
// Assuming a poissonian distribution of the collision, with:
//    mu = mean number of collisions
// the probability of having 0 collisions is:
//    exp{-mu} = 1 - L0b_MB / max collision rate
// Since the MB trigger is measuring the probability of having at least
// 1 collisions, the number should be corrected by the factor:
//    CF = mu / ( 1 - exp{-mu} )
// The class weights the number of events with this correction factor
//
// CAVEATS:
// - so far the class needs to access the OCDB
//   (hopefully it will be changed in the future)
//   Hence the following libraries need to be included (in addition to
//   the standard analysis ones):
//   libRAWDatabase.so libCDB.so libSTEER.so libPWG3base.so
//-------------------------------------------------------------------------

#include <Riostream.h>

// ROOT includes
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TString.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "TMath.h"

// STEER includes
#include "AliESDEvent.h"
#include "AliESDInputHandler.h"
#include "AliAODEvent.h"
#include "AliAODInputHandler.h"
#include "AliTriggerScalersESD.h"
#include "AliLog.h"

// ANALYSIS includes
#include "AliAnalysisTaskSE.h"
#include "AliAnalysisManager.h"
#include "AliAnalysisDataSlot.h"
#include "AliAnalysisDataContainer.h"

#include "AliCounterCollection.h"

#include "AliAnalysisTaskPileup.h"

#if defined(READOCDB)
#include "AliTriggerRunScalers.h"
#include "AliTriggerConfiguration.h"
#include "AliTriggerClass.h"
#include "AliCDBManager.h"
#include "AliCDBEntry.h"
#include "AliCDBPath.h"
#endif


ClassImp(AliAnalysisTaskPileup)


//________________________________________________________________________
AliAnalysisTaskPileup::AliAnalysisTaskPileup(const char *name) :
  AliAnalysisTaskSE(name),
  fEventCounters(0x0),
  fHistoEventsList(0x0),
  fTriggerClasses(0x0),
  fTriggerClassIndex(new TArrayI(50))
#if defined(READOCDB)
  , fTriggerRunScalers(0x0),
  fDefaultStorage(new TString())
#endif

{
  /// Constructor

  fTriggerClassIndex->Reset(-1);

  DefineOutput(1,AliCounterCollection::Class());
  DefineOutput(2,TObjArray::Class());
}

//________________________________________________________________________
AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
{
  /// Destructor

  // For proof: do not delete output containers
  if ( ! AliAnalysisManager::GetAnalysisManager()->IsProofMode() ) {
    delete fEventCounters;
  }

  delete fHistoEventsList;
  delete fTriggerClasses;
  delete fTriggerClassIndex;

#if defined(READOCDB)
  delete fTriggerRunScalers;
  delete fDefaultStorage;
#endif

}


//___________________________________________________________________________
void AliAnalysisTaskPileup::NotifyRun()
{
  /// Notify run

#if defined(READOCDB)
  if ( ! AliCDBManager::Instance()->GetDefaultStorage() ) {
    AliCDBManager::Instance()->SetDefaultStorage(fDefaultStorage->Data());
  }

  AliCDBManager::Instance()->SetRun(InputEvent()->GetRunNumber());

  AliCDBEntry *entry = 0x0;
  entry = AliCDBManager::Instance()->Get("GRP/CTP/Config");
  if ( entry ) {
    AliTriggerConfiguration* trigConf = (AliTriggerConfiguration*)entry->GetObject();
    const TObjArray& classesArray = trigConf->GetClasses();
    Int_t nclasses = classesArray.GetEntriesFast();

    if ( fTriggerClasses ) delete fTriggerClasses;

    fTriggerClasses = new TObjArray(nclasses);
    fTriggerClasses->SetOwner();

    Int_t currActive = -1;
    for( Int_t iclass=0; iclass < nclasses; iclass++ ) {
      AliTriggerClass* trclass = (AliTriggerClass*)classesArray.At(iclass);
      if (trclass && trclass->GetMask()>0) {
	currActive++;
	Int_t currPos = currActive;

	// Store the CBEAMB class at the first position
	TString trigName = trclass->GetName();
	if ( ( trigName.Contains("CBEAMB") ||  trigName.Contains("CTRUE") ) && ! trigName.Contains("WU") )
	  currPos = 0;
	else if ( ! fTriggerClasses->At(0) )
	  currPos++;

	Int_t trindex = TMath::Nint(TMath::Log2(trclass->GetMask()));
	TObjString* objString = new TObjString(trigName.Data());
	fTriggerClasses->AddAtAndExpand(objString, currPos);
	(*fTriggerClassIndex)[currPos] = trindex;
	AliDebug(3, Form("Current class %s  index %i  position %i", trigName.Data(), trindex, currPos));
      } // is class active
    } // loop on trigger classes
  } // if entry

  entry = AliCDBManager::Instance()->Get("GRP/CTP/Scalers");
  if (entry) {    
       AliInfo("Found an AliTriggerRunScalers in GRP/CTP/Scalers, reading it");
       fTriggerRunScalers = dynamic_cast<AliTriggerRunScalers*> (entry->GetObject());
       entry->SetOwner(0);
       if (fTriggerRunScalers->CorrectScalersOverflow() == 0) AliInfo("32bit Trigger counters corrected for overflow");
  }
#endif

}


//___________________________________________________________________________
void AliAnalysisTaskPileup::UserCreateOutputObjects()
{
  /// Create histograms and counters
  
  // The framework has problems if the name of the object
  // and the one of container differ
  // To be complaint, get the name from container and set it
  TString containerName = GetOutputSlot(1)->GetContainer()->GetName();

  // initialize event counters
  fEventCounters = new AliCounterCollection(containerName.Data());
  fEventCounters->AddRubric("event", "any/correctedL0");
  fEventCounters->AddRubric("trigger", 1000000);
  fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
  fEventCounters->AddRubric("selected", "yes/no");
  fEventCounters->Init(kTRUE);

  // Post data at least once per task to ensure data synchronisation (required for merging)
  PostData(1, fEventCounters);
}

//________________________________________________________________________
void AliAnalysisTaskPileup::UserExec(Option_t *)
{
  /// Called for each event

  AliAODEvent* aodEvent = 0x0;

  AliESDEvent* esdEvent = dynamic_cast<AliESDEvent*> (InputEvent());
  if ( ! esdEvent ) {
    aodEvent = dynamic_cast<AliAODEvent*> (InputEvent());
  }

  if ( ! aodEvent && ! esdEvent ) {
    AliError ("AOD or ESD event not found. Nothing done!");
    return;
  }
 
  // check physics selection
  Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
  TString selected = ( isPhysicsSelected ) ? "yes" : "no";

#if defined(READOCDB)

  TString firedTrigClasses = ( esdEvent ) ? esdEvent->GetFiredTriggerClasses() : aodEvent->GetFiredTriggerClasses();

  if ( firedTrigClasses.IsNull() ) return; // Reject non-physics events

  AliDebug(2, Form("\nEvent %lli\n", Entry()));

  Int_t nPoints = fTriggerRunScalers->GetScalersRecordsESD()->GetEntriesFast();

  // Add protection for MC (no scalers there!)
  AliTriggerScalersRecordESD* trigScalerRecords1 = 0x0;
  AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;
  if ( nPoints > 1 ) {

    AliTimeStamp timeStamp(InputEvent()->GetOrbitNumber(), InputEvent()->GetPeriodNumber(), InputEvent()->GetBunchCrossNumber());
    Int_t position = fTriggerRunScalers->FindNearestScalersRecord(&timeStamp);
    if ( position < 0 ) {
      AliWarning("Position out of range: put to 1");
      position = 1;
    } 
    if ( position == 0 ) position++; // Don't trust the first one
    else if ( position + 1 >= nPoints ) position--;
    AliDebug(2, Form("position %i\n", position));
    trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);

    // Sometimes scalers are filled very close to each others
    // in this case skip to the next entry
    for ( Int_t ipos=position+1; ipos<nPoints; ipos++ ) {
      trigScalerRecords2 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(ipos);
      Double_t deltaTime = (Double_t)( trigScalerRecords2->GetTimeStamp()->GetSeconds() - trigScalerRecords1->GetTimeStamp()->GetSeconds() );
      AliDebug(2, Form("Pos %i  TimeStamp %u - %u = %.0f\n", ipos, trigScalerRecords2->GetTimeStamp()->GetSeconds(), trigScalerRecords1->GetTimeStamp()->GetSeconds(), deltaTime));

      if ( deltaTime > 1 )
	break;
    } // loop on position
  } // nPoins > 0

#endif

  ULong64_t trigMask = 0;
  Double_t correctFactor = 1.;

  Int_t nVtxContrib = ( esdEvent ) ? esdEvent->GetPrimaryVertex()->GetNContributors() : aodEvent->GetPrimaryVertex()->GetNContributors();
  Bool_t isPileupSPD = ( esdEvent ) ? esdEvent->IsPileupFromSPD(3,0.8) : aodEvent->IsPileupFromSPD(3,0.8);

  TString vtxSelKey[3] = {"any","hasVtxContrib","nonPileupSPD"};
  Bool_t fillSel[3] = {kTRUE, ( nVtxContrib > 0 ), ( ( nVtxContrib > 0 ) && ( ! isPileupSPD ) )};

  //const AliTriggerScalersRecordESD* trigScalerRecords = esdEvent->GetHeader()->GetTriggerScalersRecord(); // REMEMBER TO CUT

  TString trigName = "";
  TString eventType = "";

  Int_t nTriggerClasses = fTriggerClasses->GetEntries();
  Int_t classIndex = -1;
  Double_t deltaScalersBeam = 0., deltaScalers = 0.;
  Bool_t isFiredOnce = kFALSE;
  for (Int_t itrig=0; itrig<nTriggerClasses+1; itrig++) {

    Double_t correctFactorL0 = 1.;

    Bool_t isClassFired = kFALSE;

    if ( itrig < nTriggerClasses ) {

      // Check if current mask contains trigger
      classIndex = (*fTriggerClassIndex)[itrig];
      if ( classIndex < 0 ) continue; // Protection for MC (where BEAMB not present
      trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
      trigMask = ( 1ull << classIndex );
      isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );

#if defined(READOCDB)
      if ( trigScalerRecords2 && ( isClassFired || itrig == 0 ) ) {
	// Get scalers
	const AliTriggerScalersESD* scaler1 = trigScalerRecords1->GetTriggerScalersForClass(classIndex+1);
	const AliTriggerScalersESD* scaler2 = trigScalerRecords2->GetTriggerScalersForClass(classIndex+1);
	deltaScalers = scaler2->GetLOCB() - scaler1->GetLOCB();
 
	if ( itrig == 0 )
	  deltaScalersBeam = deltaScalers;
	else if ( isClassFired ) {
	  correctFactorL0 = GetL0Correction(deltaScalers, deltaScalersBeam);
	  AliDebug(2, Form("Scalers: %s %.0f  %s %.0f -> CF %f\n", fTriggerClasses->At(itrig)->GetName(), deltaScalers, fTriggerClasses->At(0)->GetName(), deltaScalersBeam, correctFactorL0));
	}
      }
#endif
    } // if ( itrig < nTriggerClasses )
    else {
      classIndex = -1;
      trigName = "any";
      isClassFired = isFiredOnce;
    }

    if ( ! isClassFired ) continue;
    isFiredOnce = kTRUE;

    //const AliTriggerScalersESD* trigScaler = trigScalerRecords->GetTriggerScalersForClass(classIndex+1); // REMEMBER TO CUT
    //if ( classIndex > 1 ) printf("Index: trigger %i  scaler %i\n", classIndex+1, trigScaler->GetClassIndex()); // REMEMBER TO CUT

    AliDebug(2, Form("Fired trig %s\n", trigName.Data()));

    for ( Int_t ievType=0; ievType<2; ievType++ ){
      switch ( ievType ) {
#if defined(READOCDB)
      case kHeventsCorrectL0:
	correctFactor = correctFactorL0;
	eventType = "correctedL0";
	break;
#endif
      default:
	correctFactor = 1.;
	eventType = "any";
      }

      for ( Int_t isel=0; isel<3; isel++ ) {
	if ( ! fillSel[isel] ) continue;
	fEventCounters->Count(Form("event:%s/trigger:%s/vtxSelection:%s/selected:%s",eventType.Data(),trigName.Data(), vtxSelKey[isel].Data(), selected.Data()),correctFactor);
      } // loop on vertex selection
    } // loop on event type
  } // loop on trigger classes

  // Post final data. It will be written to a file with option "RECREATE"
  PostData(1, fEventCounters);
}


//________________________________________________________________________
void AliAnalysisTaskPileup::Terminate(Option_t *)
{
  //
  /// Save final histograms
  /// and draw result to the screen
  //

  // Fill the container only at the very last step
  // i.e. in local, when done interactively
  if ( gROOT->IsBatch() )
    return;

  fEventCounters = dynamic_cast<AliCounterCollection*>(GetOutputData(1));
  if ( ! fEventCounters ) return;

  if ( ! fHistoEventsList ) fHistoEventsList = new TObjArray(0);
  fHistoEventsList->SetOwner();

  TH2D* histo = 0x0;
  const Int_t kNevtTypes = 2;
  TString evtSel[kNevtTypes] = {"event:any", "event:correctedL0"};
  TString evtName[kNevtTypes] = {"", "CorrectL0"};
  TString evtTitle[kNevtTypes] = {"Events", "L0 corrected events"};
  const Int_t kNphysSel = 2;
  TString physSel[kNphysSel] = {"selected:any","selected:yes"};
  TString physName[kNphysSel] = {"", "PhysSel"};
  TString physTitle[kNphysSel] = {"", "w/ physics selection"};

  Int_t currHisto = -1;
  TString currName = "";
  for ( Int_t isel=0; isel<kNphysSel; isel++ ) {
    for ( Int_t iev=0; iev<kNevtTypes; iev++ ) {
      currName = Form("%s/%s", evtSel[iev].Data(), physSel[isel].Data());
      histo = fEventCounters->Get("trigger","vtxSelection",currName.Data());
      if ( ! histo ) continue;
      currHisto++;
      currName = Form("hEvents%s%s", evtName[iev].Data(), physName[isel].Data());
      histo->SetName(currName.Data());
      currName = Form("%s %s", evtTitle[iev].Data(), physTitle[isel].Data());
      histo->SetTitle(currName.Data());
      fHistoEventsList->AddAtAndExpand(histo, currHisto);
    } // loop on event type
    TH2D* num = (TH2D*)fHistoEventsList->At(currHisto);
    TH2D* den = (TH2D*)fHistoEventsList->At(currHisto-1);
    if ( ! num || ! den ) continue;
    currName = Form("hPileupL0%s_%s", physName[isel].Data(), GetName());
    TH2D* histoPileupL0 = (TH2D*)num->Clone(currName.Data());
    histoPileupL0->Divide(den);
    currName = Form("c%i_%s", isel, GetName());
    TCanvas *can = new TCanvas(currName.Data(),"Pileup",10,10,310,310);
    can->SetFillColor(10); can->SetHighLightColor(10);
    can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
    histoPileupL0->DrawCopy("text");
    delete histoPileupL0;
  } // loop on physics selection

  PostData(2, fHistoEventsList);
}


//________________________________________________________________________
Double_t AliAnalysisTaskPileup::GetL0Correction(Double_t nCINT1B, Double_t nCBEAMB)
{
  /// Get the correction factor for L0 calculation

  if ( nCBEAMB == 0. )
    return 1.;

  Double_t ratio = nCINT1B / nCBEAMB;

  if ( ratio >= 1. || ratio == 0. )
    return 1.;

  Double_t mu = -TMath::Log(1-ratio);

  return mu / ( 1. - TMath::Exp(-mu) );

}
Commit	Line	Data
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
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
c705eb52	106	fTriggerClassIndex->Reset(-1);
c705eb52	107
0fa651e5	108	DefineOutput(1,AliCounterCollection::Class());
	109	DefineOutput(2,TObjArray::Class());
	110	}
	111
	112	//________________________________________________________________________
	113	AliAnalysisTaskPileup::~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	//___________________________________________________________________________
	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();
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	//___________________________________________________________________________
	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);
c705eb52	207	fEventCounters->AddRubric("vtxSelection", "any/hasVtxContrib/nonPileupSPD");
c705eb52	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	//________________________________________________________________________
	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
c705eb52	233	Bool_t isPhysicsSelected = (fInputHandler && fInputHandler->IsEventSelected());
c705eb52	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!)
c705eb52	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
c705eb52	305	trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
0fa651e5	306	trigMask = ( 1ull << classIndex );
	307	isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );
	308
c705eb52	309	#if defined(READOCDB)
c705eb52	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
c705eb52	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	//________________________________________________________________________
	365	void 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
0fa651e5	421	PostData(2, fHistoEventsList);
0fa651e5	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