[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

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

//________________________________________________________________________
AliAnalysisTaskPileup::~AliAnalysisTaskPileup()
{
  /// Destructor
  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("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("selection", "any/hasVtxContrib/nonPileupSPD");
  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());

#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();

  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));
  AliTriggerScalersRecordESD* trigScalerRecords1 = (AliTriggerScalersRecordESD*)fTriggerRunScalers->GetScalersRecordsESD()->At(position);
  AliTriggerScalersRecordESD* trigScalerRecords2 = 0x0;

  // 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;
  }

#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 selKey[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
      trigName = ((TObjString*)fTriggerClasses->At(itrig))->GetString();
      classIndex = (*fTriggerClassIndex)[itrig];
      trigMask = ( 1ull << classIndex );
      isClassFired = ( trigMask & InputEvent()->GetTriggerMask() );

      if ( 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));
	}
      }
    }
    else {
      trigName = "any";
      classIndex = -1;
      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/selection:%s",eventType.Data(),trigName.Data(), selKey[isel].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
  //

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

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

  TH2D* histo = 0x0;
  histo = fEventCounters->Draw("trigger","selection","event:any");
  if ( histo ) {
    histo->SetName("hEvents");
    histo->SetTitle("Events per trigger");
    fHistoEventsList->AddAtAndExpand(histo, kHevents);
  }

  histo = fEventCounters->Draw("trigger","selection","event:correctedL0");
  if ( histo ) {
    histo->SetName("hEventsCorrectL0");
    histo->SetTitle("L0 corrected events");
    fHistoEventsList->AddAtAndExpand(histo, kHeventsCorrectL0);
  }

  PostData(2, fHistoEventsList);

  if ( gROOT->IsBatch() )
    return;

  TH2D* histoPileupL0 = (TH2D*)fHistoEventsList->At(kHeventsCorrectL0)->Clone("hPileupL0");
  histoPileupL0->Divide((TH2D*)fHistoEventsList->At(kHevents));
  
  TCanvas *can = new TCanvas("can1_Pileup","Pileup",10,10,310,310);
  can->SetFillColor(10); can->SetHighLightColor(10);
  can->SetLeftMargin(0.15); can->SetBottomMargin(0.15);
  histoPileupL0->DrawCopy("text");
}


//________________________________________________________________________
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
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