PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id:$ */
  17
  18 #include <TROOT.h>
  19 #include <TChain.h>
  20 #include <TFile.h>
  21 #include <TList.h>
  22 #include <TMath.h>
  23 #include <TTree.h>
  24 #include <TH2F.h>
  25 #include <TH3F.h>
  26 #include <TRandom.h>
  27
  28 #include "AliAnalysisTaskPhiCorrelations.h"
  29 #include "AliAnalyseLeadingTrackUE.h"
  30 #include "AliUEHistograms.h"
  31 #include "AliUEHist.h"
  32
  33 #include "AliAnalysisManager.h"
  34 #include "AliAODHandler.h"
  35 #include "AliAODInputHandler.h"
  36 #include "AliAODMCParticle.h"
  37 #include "AliInputEventHandler.h"
  38 #include "AliLog.h"
  39 #include "AliMCEventHandler.h"
  40 #include "AliVParticle.h"
  41 #include "AliCFContainer.h"
  42
  43 #include "AliESDEvent.h"
  44 #include "AliESDInputHandler.h"
  45 #include "AliMultiplicity.h"
  46 #include "AliCentrality.h"
  47 #include "AliStack.h"
  48 #include "AliAODMCHeader.h"
  49 #include "AliGenCocktailEventHeader.h"
  50 #include "AliGenEventHeader.h"
  51
  52 #include "AliEventPoolManager.h"
  53
  54 #include "AliESDZDC.h"
  55 #include "AliESDtrackCuts.h"
  56
  57 #include "AliHelperPID.h"
  58
  59 ////////////////////////////////////////////////////////////////////////
  60 //
  61 // Analysis class for azimuthal correlation studies
  62 // Based on the UE task from Sara Vallero and Jan Fiete
  63 //
  64 // This class needs input AODs.
  65 // The output is a list of analysis-specific containers.
  66 //
  67 // The AOD can be either connected to the InputEventHandler
  68 // for a chain of AOD files
  69 // or
  70 // to the OutputEventHandler
  71 // for a chain of ESD files,
  72 // in this case the class should be in the train after the jet-finder
  73 //
  74 //    Authors:
  75 //    Jan Fiete Grosse-Oetringhaus
  76 //
  77 ////////////////////////////////////////////////////////////////////////
  78
  79
  80 ClassImp( AliAnalysisTaskPhiCorrelations )
  81 ClassImp( AliDPhiBasicParticle )
  82
  83 //____________________________________________________________________
  84 AliAnalysisTaskPhiCorrelations:: AliAnalysisTaskPhiCorrelations(const char* name):
  85 AliAnalysisTask(name,""),
  86 // general configuration
  87 fDebug(0),
  88 fMode(0),
  89 fReduceMemoryFootprint(kFALSE),
  90 fFillMixed(kTRUE),
  91 fMixingTracks(50000),
  92 fCompareCentralities(kFALSE),
  93 fTwoTrackEfficiencyStudy(kFALSE),
  94 fTwoTrackEfficiencyCut(0),
  95 fUseVtxAxis(kFALSE),
  96 fCourseCentralityBinning(kFALSE),
  97 fSkipTrigger(kFALSE),
  98 fInjectedSignals(kFALSE),
  99 // pointers to UE classes
 100 fHelperPID(0x0),
 101 fAnalyseUE(0x0),
 102 fHistos(0x0),
 103 fHistosMixed(0),
 104 fEfficiencyCorrectionTriggers(0),
 105 fEfficiencyCorrectionAssociated(0),
 106 // handlers and events
 107 fAOD(0x0),
 108 fESD(0x0),
 109 fArrayMC(0x0),
 110 fInputHandler(0x0),
 111 fMcEvent(0x0),
 112 fMcHandler(0x0),
 113 fPoolMgr(0x0),
 114 // histogram settings
 115 fListOfHistos(0x0),
 116 // event QA
 117 fnTracksVertex(1),  // QA tracks pointing to principal vertex (= 3 default)
 118 fZVertex(7.),
 119 fCentralityMethod("V0M"),
 120 // track cuts
 121 fTrackEtaCut(0.8),
 122 fOnlyOneEtaSide(0),
 123 fPtMin(0.5),
 124 fFilterBit(0xFF),
 125 fTrackStatus(0),
 126 fSelectBit(AliVEvent::kMB|AliVEvent::kUserDefined),
 127 fUseChargeHadrons(kFALSE),
 128 fParticleSpeciesTrigger(-1),
 129 fParticleSpeciesAssociated(-1),
 130 fCheckMotherPDG(kTRUE),
 131 fSelectCharge(0),
 132 fTriggerSelectCharge(0),
 133 fAssociatedSelectCharge(0),
 134 fTriggerRestrictEta(-1),
 135 fEtaOrdering(kFALSE),
 136 fCutConversions(kFALSE),
 137 fCutResonances(kFALSE),
 138 fFillOnlyStep0(kFALSE),
 139 fSkipStep6(kFALSE),
 140 fRejectCentralityOutliers(kFALSE),
 141 fRemoveWeakDecays(kFALSE),
 142 fRemoveDuplicates(kFALSE),
 143 fSkipFastCluster(kFALSE),
 144 fWeightPerEvent(kFALSE),
 145 fCustomBinning(),
 146 fPtOrder(kTRUE),
 147 fFillpT(kFALSE)
 148 {
 149   // Default constructor
 150   // Define input and output slots here
 151   // Input slot #0 works with a TChain
 152   DefineInput(0, TChain::Class());
 153   // Output slot #0 writes into a TList container
 154   DefineOutput(0, TList::Class());
 155 }
 156
 157 AliAnalysisTaskPhiCorrelations::~AliAnalysisTaskPhiCorrelations()
 158 {
 159   // destructor
 160
 161   if (fListOfHistos  && !AliAnalysisManager::GetAnalysisManager()->IsProofMode())
 162     delete fListOfHistos;
 163 }
 164
 165 //____________________________________________________________________
 166 void AliAnalysisTaskPhiCorrelations::ConnectInputData(Option_t* /*option*/)
 167 {
 168
 169   // Connect the input data
 170   if (fDebug > 1) AliInfo("ConnectInputData() ");
 171
 172   // Since AODs can either be connected to the InputEventHandler
 173   // or to the OutputEventHandler ( the AOD is created by a previus task in the train )
 174   // we need to get the pointer to the AODEvent correctly.
 175
 176   // Delta AODs are also accepted.
 177
 178   TObject* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 179
 180   if( handler && handler->InheritsFrom("AliAODInputHandler") )
 181   { // input AOD
 182     fAOD = ((AliAODInputHandler*)handler)->GetEvent();
 183     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODInputHandler");
 184   }
 185   else
 186   {  //output AOD
 187     handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 188     if (handler && handler->InheritsFrom("AliAODHandler") )
 189     {
 190       fAOD = ((AliAODHandler*)handler)->GetAOD();
 191       if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliAODHandler");
 192     }
 193     else
 194     {  // no AOD
 195       AliWarning("I can't get any AOD Event Handler");
 196     }
 197   }
 198
 199   if (handler && handler->InheritsFrom("AliESDInputHandler") )
 200   { // input ESD
 201     // pointer received per event in ::Exec
 202     if (fDebug > 1) AliInfo(" ==== Tracks and Jets from AliESDInputHandler");
 203   }
 204
 205   // Initialize common pointers
 206   Initialize();
 207 }
 208
 209 //____________________________________________________________________
 210 void  AliAnalysisTaskPhiCorrelations::CreateOutputObjects()
 211 {
 212   // Create the output container
 213
 214   if (fDebug > 1) AliInfo("CreateOutputObjects()");
 215
 216   // Initialize class with main algorithms, event and track selection.
 217   fAnalyseUE = new AliAnalyseLeadingTrackUE();
 218   fAnalyseUE->SetParticleSelectionCriteria(fFilterBit, fUseChargeHadrons, fTrackEtaCut, fPtMin);
 219   fAnalyseUE->SetTrackStatus(fTrackStatus);
 220   fAnalyseUE->SetCheckMotherPDG(fCheckMotherPDG);
 221   fAnalyseUE->SetDebug(fDebug);
 222   fAnalyseUE->DefineESDCuts(fFilterBit);
 223   fAnalyseUE->SetEventSelection(fSelectBit);
 224   fAnalyseUE->SetHelperPID(fHelperPID);
 225   if ((fParticleSpeciesTrigger != -1 || fParticleSpeciesAssociated != -1) && !fHelperPID)
 226     AliFatal("HelperPID object should be set in the steering macro");
 227
 228   // Initialize output list of containers
 229   if (fListOfHistos != NULL){
 230         delete fListOfHistos;
 231         fListOfHistos = NULL;
 232         }
 233   if (!fListOfHistos){
 234         fListOfHistos = new TList();
 235         fListOfHistos->SetOwner(kTRUE);
 236         }
 237
 238   // Initialize class to handle histograms
 239   TString histType = "4R";
 240   if (fUseVtxAxis == 1)
 241     histType = "5R";
 242   else if (fUseVtxAxis == 2)
 243     histType = "6R";
 244   if (fCourseCentralityBinning)
 245     histType += "C";
 246   fHistos = new AliUEHistograms("AliUEHistogramsSame", histType, fCustomBinning);
 247   fHistosMixed = new AliUEHistograms("AliUEHistogramsMixed", histType, fCustomBinning);
 248
 249   fHistos->SetSelectCharge(fSelectCharge);
 250   fHistosMixed->SetSelectCharge(fSelectCharge);
 251
 252   fHistos->SetSelectTriggerCharge(fTriggerSelectCharge);
 253   fHistosMixed->SetSelectTriggerCharge(fTriggerSelectCharge);
 254
 255   fHistos->SetSelectAssociatedCharge(fAssociatedSelectCharge);
 256   fHistosMixed->SetSelectAssociatedCharge(fAssociatedSelectCharge);
 257
 258   fHistos->SetTriggerRestrictEta(fTriggerRestrictEta);
 259   fHistosMixed->SetTriggerRestrictEta(fTriggerRestrictEta);
 260
 261   fHistos->SetOnlyOneEtaSide(fOnlyOneEtaSide);
 262   fHistosMixed->SetOnlyOneEtaSide(fOnlyOneEtaSide);
 263
 264   fHistos->SetEtaOrdering(fEtaOrdering);
 265   fHistosMixed->SetEtaOrdering(fEtaOrdering);
 266
 267   fHistos->SetPairCuts(fCutConversions, fCutResonances);
 268   fHistosMixed->SetPairCuts(fCutConversions, fCutResonances);
 269
 270   fHistos->SetTrackEtaCut(fTrackEtaCut);
 271   fHistosMixed->SetTrackEtaCut(fTrackEtaCut);
 272
 273   fHistos->SetWeightPerEvent(fWeightPerEvent);
 274   fHistosMixed->SetWeightPerEvent(fWeightPerEvent);
 275
 276   fHistos->SetPtOrder(fPtOrder);
 277   fHistosMixed->SetPtOrder(fPtOrder);
 278
 279   if (fEfficiencyCorrectionTriggers)
 280    {
 281     fHistos->SetEfficiencyCorrectionTriggers(fEfficiencyCorrectionTriggers);
 282     fHistosMixed->SetEfficiencyCorrectionTriggers((THnF*) fEfficiencyCorrectionTriggers->Clone());
 283    }
 284   if (fEfficiencyCorrectionAssociated)
 285   {
 286     fHistos->SetEfficiencyCorrectionAssociated(fEfficiencyCorrectionAssociated);
 287     fHistosMixed->SetEfficiencyCorrectionAssociated((THnF*) fEfficiencyCorrectionAssociated->Clone());
 288   }
 289
 290   // add histograms to list
 291   fListOfHistos->Add(fHistos);
 292   fListOfHistos->Add(fHistosMixed);
 293   // add HelperPID to list
 294   if (fHelperPID)
 295     fListOfHistos->Add(fHelperPID);
 296
 297   fListOfHistos->Add(new TH2F("trackletsVsV0Cent", ";L1 clusters;v0 centrality", 100, -0.5, 9999.5, 101, 0, 101));
 298   fListOfHistos->Add(new TH2F("processIDs", ";#Delta#phi;process id", 100, -0.5 * TMath::Pi(), 1.5 * TMath::Pi(), kPNoProcess + 1, -0.5, kPNoProcess + 0.5));
 299   fListOfHistos->Add(new TH1F("eventStat", ";;events", 4, -0.5, 3.5));
 300   fListOfHistos->Add(new TH2F("mixedDist", ";centrality;tracks;events", 101, 0, 101, 200, 0, fMixingTracks * 1.5));
 301   fListOfHistos->Add(new TH1F("pids", ";pdg;tracks", 2001, -1000.5, 1000.5));
 302   fListOfHistos->Add(new TH2F("referenceMultiplicity", ";centrality;tracks;events", 101, 0, 101, 200, 0, 200));
 303
 304   PostData(0,fListOfHistos);
 305
 306   // Add task configuration to output list
 307   AddSettingsTree();
 308
 309   // event mixing
 310   Int_t poolsize   = 1000;  // Maximum number of events, ignored in the present implemention of AliEventPoolManager
 311
 312   Int_t nCentralityBins  = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(1);
 313   Double_t* centralityBins = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(1, 0)->GetXbins()->GetArray();
 314
 315   const Int_t kNZvtxBins  = 10+(1+10)*4;
 316   // bins for further buffers are shifted by 100 cm
 317   Double_t vertexBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10,
 318                                        90,  92,  94,  96,  98, 100, 102, 104, 106, 108, 110,
 319                                       190, 192, 194, 196, 198, 200, 202, 204, 206, 208, 210,
 320                                       290, 292, 294, 296, 298, 300, 302, 304, 306, 308, 310,
 321                                       390, 392, 394, 396, 398, 400, 402, 404, 406, 408, 410 };
 322
 323   Int_t nZvtxBins  = kNZvtxBins;
 324   Double_t* zvtxbin = vertexBins;
 325
 326   if (fMode == 0 && fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
 327   {
 328     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
 329     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
 330   }
 331
 332   fPoolMgr = new AliEventPoolManager(poolsize, fMixingTracks, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
 333   fPoolMgr->SetTargetValues(fMixingTracks, 0.1, 5);
 334 }
 335
 336 //____________________________________________________________________
 337 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
 338 {
 339   // receive ESD pointer if we are not running AOD analysis
 340   if (!fAOD)
 341   {
 342     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 343     if (handler && handler->InheritsFrom("AliESDInputHandler"))
 344       fESD = ((AliESDInputHandler*)handler)->GetEvent();
 345   }
 346
 347   if (fMode)
 348   {
 349     // correction mode
 350
 351     if (fMcHandler)
 352       fMcEvent = fMcHandler->MCEvent();
 353
 354     if (fAOD)
 355     {
 356       // array of MC particles
 357       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 358       if (!fArrayMC)
 359         AliFatal("No array of MC particles found !!!");
 360     }
 361
 362     AnalyseCorrectionMode();
 363   }
 364   else AnalyseDataMode();
 365 }
 366
 367 /******************** ANALYSIS METHODS *****************************/
 368
 369 //____________________________________________________________________
 370 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
 371 {
 372   //Write settings to output list
 373   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
 374   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
 375   settingsTree->Branch("fZVertex", &fZVertex,"ZVertex/D");
 376   //settingsTree->Branch("fCentralityMethod", fCentralityMethod.Data(),"CentralityMethod/C");
 377   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
 378   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
 379   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
 380   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
 381   settingsTree->Branch("fTrackStatus", &fTrackStatus,"TrackStatus/I");
 382   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
 383   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
 384   settingsTree->Branch("fParticleSpeciesTrigger", &fParticleSpeciesTrigger,"ParticleSpeciesTrigger/I");
 385   settingsTree->Branch("fParticleSpeciesAssociated", &fParticleSpeciesAssociated,"ParticleSpeciesAssociated/I");
 386   settingsTree->Branch("fCheckMotherPDG", &fCheckMotherPDG,"CheckMotherPDG/I");
 387   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
 388   settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
 389   settingsTree->Branch("fAssociatedSelectCharge", &fAssociatedSelectCharge,"fAssociatedSelectCharge/I");
 390   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
 391   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
 392   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
 393   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
 394   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
 395   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
 396   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
 397   settingsTree->Branch("fInjectedSignals", &fInjectedSignals,"SkipTrigger/O");
 398   settingsTree->Branch("fRejectCentralityOutliers", &fRejectCentralityOutliers,"RejectCentralityOutliers/O");
 399   settingsTree->Branch("fRemoveWeakDecays", &fRemoveWeakDecays,"RemoveWeakDecays/O");
 400   settingsTree->Branch("fRemoveDuplicates", &fRemoveDuplicates,"RemoveDuplicates/O");
 401   settingsTree->Branch("fSkipFastCluster", &fSkipFastCluster,"SkipFastCluster/O");
 402   settingsTree->Branch("fWeightPerEvent", &fWeightPerEvent,"WeightPerEvent/O");
 403   settingsTree->Branch("fPtOrder", &fPtOrder,"PtOrder/O");
 404   //fCustomBinning
 405
 406   settingsTree->Fill();
 407   fListOfHistos->Add(settingsTree);
 408 }
 409
 410 //____________________________________________________________________
 411 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
 412 {
 413   // Run the analysis on MC to get the correction maps
 414   //
 415
 416   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
 417
 418   Double_t centrality = 0;
 419
 420   if (fCentralityMethod.Length() > 0)
 421   {
 422     AliCentrality *centralityObj = 0;
 423     if (fAOD)
 424       centralityObj = fAOD->GetHeader()->GetCentralityP();
 425     else if (fESD)
 426       centralityObj = fESD->GetCentrality();
 427
 428     if (centralityObj)
 429     {
 430       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
 431       AliInfo(Form("Centrality is %f", centrality));
 432     }
 433     else
 434     {
 435       Printf("WARNING: Centrality object is 0");
 436       centrality = -1;
 437      }
 438   }
 439
 440   // Support for ESD and AOD based analysis
 441   AliVEvent* inputEvent = fAOD;
 442   if (!inputEvent)
 443     inputEvent = fESD;
 444
 445   Float_t bSign = 0;
 446
 447   if (inputEvent)
 448   {
 449     fHistos->SetRunNumber(inputEvent->GetRunNumber());
 450     bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
 451   }
 452
 453   TObject* mc = fArrayMC;
 454   if (!mc)
 455     mc = fMcEvent;
 456
 457   // count all events
 458   fHistos->FillEvent(centrality, -1);
 459
 460   if (centrality < 0)
 461     return;
 462
 463   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
 464   TObject* vertexSupplier = fMcEvent;
 465   if (fAOD) // AOD
 466     vertexSupplier = fAOD->FindListObject(AliAODMCHeader::StdBranchName());
 467
 468   if (!fAnalyseUE->VertexSelection(vertexSupplier, 0, fZVertex))
 469     return;
 470
 471   Float_t zVtx = 0;
 472   if (fAOD)
 473     zVtx = ((AliAODMCHeader*) vertexSupplier)->GetVtxZ();
 474   else
 475     zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
 476   Float_t weight = 1;
 477   if (fFillpT)
 478     weight = -1;
 479
 480   // For productions with injected signals, figure out above which label to skip particles/tracks
 481   Int_t skipParticlesAbove = 0;
 482   if (fInjectedSignals)
 483   {
 484     AliGenEventHeader* eventHeader = 0;
 485     Int_t headers = 0;
 486
 487     if (fMcEvent)
 488     {
 489       // ESD
 490       AliHeader* header = (AliHeader*) fMcEvent->Header();
 491       if (!header)
 492         AliFatal("fInjectedSignals set but no MC header found");
 493
 494       AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
 495       if (!cocktailHeader)
 496       {
 497         header->Dump();
 498         AliFatal("fInjectedSignals set but no MC cocktail header found");
 499       }
 500
 501       headers = cocktailHeader->GetHeaders()->GetEntries();
 502       eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
 503     }
 504     else
 505     {
 506       // AOD
 507       AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 508       if (!header)
 509         AliFatal("fInjectedSignals set but no MC header found");
 510
 511       headers = header->GetNCocktailHeaders();
 512       eventHeader = header->GetCocktailHeader(0);
 513     }
 514
 515     if (!eventHeader)
 516     {
 517       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 518       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 519       AliError("First event header not found. Skipping this event.");
 520       fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
 521       return;
 522     }
 523
 524     skipParticlesAbove = eventHeader->NProduced();
 525     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 526   }
 527
 528   // Get MC primaries
 529   // triggers
 530   TObjArray* tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
 531   CleanUp(tmpList, mc, skipParticlesAbove);
 532   TObjArray* tracksMC = CloneAndReduceTrackList(tmpList);
 533   delete tmpList;
 534
 535   // associated
 536   TObjArray* tracksCorrelateMC = tracksMC;
 537   if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
 538   {
 539     // TODO for MC this uses the PDG of the mother of the particle
 540     tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
 541     CleanUp(tmpList, mc, skipParticlesAbove);
 542     tracksCorrelateMC = CloneAndReduceTrackList(tmpList);
 543     delete tmpList;
 544   }
 545
 546   /*
 547   if (fAOD)
 548   {
 549     for (Int_t i=0; i<fArrayMC->GetEntriesFast(); i++)
 550       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(((AliAODMCParticle*) fArrayMC->At(i))->PdgCode());
 551   }
 552   else
 553   {
 554     for (Int_t i=0; i<fMcEvent->GetNumberOfTracks(); i++)
 555       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(fMcEvent->GetTrack(i)->PdgCode());
 556   }
 557   */
 558
 559   if (fFillOnlyStep0)
 560     zVtx = 0;
 561
 562   // (MC-true all particles)
 563   // STEP 0
 564   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, tracksCorrelateMC, weight);
 565
 566   // mixed event
 567   if (fFillMixed)
 568   {
 569     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
 570     if (pool->IsReady())
 571       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 572         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 573     pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateMC));
 574   }
 575
 576 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
 577
 578   // Trigger selection ************************************************
 579   if (!fFillOnlyStep0 && (fSkipTrigger || fAnalyseUE->TriggerSelection(fInputHandler)))
 580   {
 581     // (MC-true all particles)
 582     // STEP 1
 583     if (!fReduceMemoryFootprint)
 584       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, tracksCorrelateMC, weight);
 585     else
 586       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
 587
 588     if (!inputEvent) {
 589       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
 590       return;
 591     }
 592
 593     // Vertex selection *************************************************
 594     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
 595     {
 596       // fill here for tracking efficiency
 597       // loop over particle species
 598
 599       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
 600       {
 601         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
 602         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kFALSE);
 603         TObjArray* allRecoTracksMatched  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kFALSE);
 604         TObjArray* primRecoTracksMatchedPID = 0;
 605         TObjArray* allRecoTracksMatchedPID  = 0;
 606
 607         if (fHelperPID)
 608         {
 609           primRecoTracksMatchedPID = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE, kTRUE);
 610           allRecoTracksMatchedPID  = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE, kTRUE);
 611         }
 612
 613         CleanUp(primMCParticles, mc, skipParticlesAbove);
 614         CleanUp(primRecoTracksMatched, mc, skipParticlesAbove);
 615         CleanUp(allRecoTracksMatched, mc, skipParticlesAbove);
 616         CleanUp(primRecoTracksMatchedPID, mc, skipParticlesAbove);
 617         CleanUp(allRecoTracksMatchedPID, mc, skipParticlesAbove);
 618
 619         // select charges
 620         if (fTriggerSelectCharge != 0)
 621         {
 622           SelectCharge(primMCParticles);
 623           SelectCharge(primRecoTracksMatched);
 624           SelectCharge(allRecoTracksMatched);
 625           SelectCharge(primRecoTracksMatchedPID);
 626           SelectCharge(allRecoTracksMatchedPID);
 627         }
 628
 629         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, primRecoTracksMatchedPID, allRecoTracksMatchedPID, 0, particleSpecies, centrality, zVtx);
 630
 631 //      Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
 632
 633         delete primMCParticles;
 634         delete primRecoTracksMatched;
 635         delete allRecoTracksMatched;
 636       }
 637
 638       TObjArray* fakeParticles = fAnalyseUE->GetFakeParticles(inputEvent, mc, kFALSE, -1, kTRUE);
 639       CleanUp((TObjArray*) fakeParticles->At(0), mc, skipParticlesAbove);
 640       CleanUp((TObjArray*) fakeParticles->At(1), mc, skipParticlesAbove);
 641
 642       fHistos->FillTrackingEfficiency(0, 0, 0, 0, 0, (TObjArray*) fakeParticles->At(2), 0, centrality, zVtx);
 643       fHistos->FillFakePt(fakeParticles, centrality);
 644 //       Printf(">>>>> %d %d %d fakes", ((TObjArray*) fakeParticles->At(0))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(1))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(2))->GetEntriesFast());
 645       delete fakeParticles;
 646
 647       // (MC-true all particles)
 648       // STEP 2
 649       if (!fReduceMemoryFootprint)
 650         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, tracksCorrelateMC, weight);
 651       else
 652         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
 653
 654       // Get MC primaries that match reconstructed track
 655       // triggers
 656       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesTrigger, kTRUE);
 657       CleanUp(tmpList, mc, skipParticlesAbove);
 658       TObjArray* tracksRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
 659       delete tmpList;
 660
 661       // associated
 662       TObjArray* tracksCorrelateRecoMatchedPrim = tracksRecoMatchedPrim;
 663       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
 664       {
 665         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, fParticleSpeciesAssociated, kTRUE);
 666         CleanUp(tmpList, mc, skipParticlesAbove);
 667         tracksCorrelateRecoMatchedPrim = CloneAndReduceTrackList(tmpList);
 668         delete tmpList;
 669       }
 670
 671       // (RECO-matched (quantities from MC particle) primary particles)
 672       // STEP 4
 673       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, tracksCorrelateRecoMatchedPrim, weight);
 674
 675       // mixed event
 676       if (fFillMixed)
 677       {
 678         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 200);
 679         if (pool->IsReady())
 680           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 681             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 682         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedPrim));
 683       }
 684
 685       // Get MC primaries + secondaries that match reconstructed track
 686       // triggers
 687       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesTrigger, kTRUE);
 688       CleanUp(tmpList, mc, skipParticlesAbove);
 689       TObjArray* tracksRecoMatchedAll = CloneAndReduceTrackList(tmpList);
 690       delete tmpList;
 691
 692       // associated
 693       TObjArray* tracksCorrelateRecoMatchedAll = tracksRecoMatchedAll;
 694       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
 695       {
 696         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, fParticleSpeciesAssociated, kTRUE);
 697         CleanUp(tmpList, mc, skipParticlesAbove);
 698         tracksCorrelateRecoMatchedAll = CloneAndReduceTrackList(tmpList);
 699         delete tmpList;
 700       }
 701
 702       // (RECO-matched (quantities from MC particle) all particles)
 703       // STEP 5
 704       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, tracksCorrelateRecoMatchedAll, weight);
 705
 706       // mixed event
 707       if (fFillMixed)
 708       {
 709         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 300);
 710         if (pool->IsReady())
 711           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 712             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 713         pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelateRecoMatchedAll));
 714       }
 715
 716       // Get RECO tracks
 717       // triggers
 718       tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
 719       CleanUp(tmpList, mc, skipParticlesAbove);
 720       TObjArray* tracks = CloneAndReduceTrackList(tmpList);
 721       delete tmpList;
 722
 723       // associated
 724       TObjArray* tracksCorrelate = tracks;
 725       if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
 726       {
 727         tmpList = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
 728         CleanUp(tmpList, mc, skipParticlesAbove);
 729         tracksCorrelate = CloneAndReduceTrackList(tmpList);
 730         delete tmpList;
 731       }
 732
 733       // (RECO all tracks)
 734       // STEP 6
 735       if (!fSkipStep6)
 736         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight);
 737
 738       // two track cut, STEP 8
 739       if (fTwoTrackEfficiencyCut > 0)
 740         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
 741
 742       // apply correction efficiency, STEP 10
 743       if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
 744       {
 745           // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
 746         Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
 747
 748         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, tracksCorrelate, weight, kTRUE, twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
 749       }
 750
 751       // mixed event
 752       if (fFillMixed)
 753       {
 754         AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
 755         if (pool2->IsReady())
 756         {
 757           for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
 758           {
 759             // STEP 6
 760             if (!fSkipStep6)
 761               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
 762
 763             // two track cut, STEP 8
 764             if (fTwoTrackEfficiencyCut > 0)
 765               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
 766
 767             // apply correction efficiency, STEP 10
 768             if (fEfficiencyCorrectionTriggers || fEfficiencyCorrectionAssociated)
 769             {
 770               // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
 771               Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
 772
 773               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
 774             }
 775           }
 776         }
 777         pool2->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
 778       }
 779
 780       if (0 && !fReduceMemoryFootprint)
 781       {
 782         // make list of secondaries (matched with MC)
 783         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
 784         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntriesFast(); i++)
 785           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
 786             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
 787
 788         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
 789         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll, weight);
 790
 791         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
 792         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries, weight);
 793
 794         // plot delta phi vs process id of secondaries
 795         // trigger particles: primaries in 4 < pT < 10
 796         // associated particles: secondaries in 1 < pT < 10
 797
 798         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntriesFast(); i++)
 799         {
 800           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
 801
 802           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
 803             continue;
 804
 805           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntriesFast(); j++)
 806           {
 807             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
 808
 809             if (particle->Pt() < 1 || particle->Pt() > 10)
 810               continue;
 811
 812             if (particle->Pt() > triggerParticle->Pt())
 813               continue;
 814
 815             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
 816             if (deltaPhi > 1.5 * TMath::Pi())
 817               deltaPhi -= TMath::TwoPi();
 818             if (deltaPhi < -0.5 * TMath::Pi())
 819               deltaPhi += TMath::TwoPi();
 820
 821             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
 822
 823             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
 824           }
 825         }
 826
 827         delete tracksRecoMatchedSecondaries;
 828       }
 829
 830       if (tracksCorrelateRecoMatchedPrim != tracksRecoMatchedPrim)
 831         delete tracksCorrelateRecoMatchedPrim;
 832       delete tracksRecoMatchedPrim;
 833
 834       if (tracksCorrelateRecoMatchedAll != tracksRecoMatchedAll)
 835         delete tracksCorrelateRecoMatchedAll;
 836       delete tracksRecoMatchedAll;
 837
 838       if (tracksCorrelate != tracks)
 839         delete tracksCorrelate;
 840       delete tracks;
 841     }
 842   }
 843
 844   if (tracksMC != tracksCorrelateMC)
 845     delete tracksCorrelateMC;
 846   delete tracksMC;
 847 }
 848
 849 //____________________________________________________________________
 850 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
 851 {
 852
 853   // Run the analysis on DATA or MC to get raw distributions
 854
 855   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
 856
 857   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
 858
 859   if (!fInputHandler)
 860     return;
 861
 862   // skip not selected events here (the AOD is not updated for those)
 863   if (!fSkipTrigger && !(fInputHandler->IsEventSelected() & fSelectBit))
 864     return;
 865
 866   // skip fast cluster events here if requested
 867   if (fSkipFastCluster && (fInputHandler->IsEventSelected() & AliVEvent::kFastOnly))
 868     return;
 869
 870   // Support for ESD and AOD based analysis
 871   AliVEvent* inputEvent = fAOD;
 872   if (!inputEvent)
 873     inputEvent = fESD;
 874
 875   Double_t centrality = 0;
 876
 877   AliCentrality *centralityObj = 0;
 878   if (fCentralityMethod.Length() > 0)
 879   {
 880     if (fCentralityMethod == "ZNA_MANUAL")
 881     {
 882       Bool_t zna = kFALSE;
 883       for(Int_t j = 0; j < 4; ++j) {
 884         if (fESD->GetZDCData()->GetZDCTDCData(12,j) != 0) {
 885           zna = kTRUE;
 886         }
 887       }
 888
 889 //       Printf("%d %f", zna, fZNAtower[0]);
 890       if (zna)
 891       {
 892         // code from Chiara O (23.10.12)
 893         const Double_t *fZNAtower = fESD->GetZDCData()->GetZN2TowerEnergy();
 894         Float_t znacut[4] = {681., 563., 413., 191.};
 895
 896         if(fZNAtower[0]>znacut[0]) centrality = 1;
 897         else if(fZNAtower[0]>znacut[1]) centrality = 21;
 898         else if(fZNAtower[0]>znacut[2]) centrality = 41;
 899         else if(fZNAtower[0]>znacut[3]) centrality = 61;
 900         else centrality = 81;
 901       }
 902       else
 903         centrality = -1;
 904     }
 905     else if (fCentralityMethod == "TRACKS_MANUAL")
 906     {
 907       // for pp
 908       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
 909       centrality = tracks->GetEntriesFast();
 910       if (centrality > 40)
 911         centrality = 41;
 912 //       Printf("%d %f", tracks->GetEntriesFast(), centrality);
 913       delete tracks;
 914     }
 915     else
 916     {
 917       if (fAOD)
 918         centralityObj = fAOD->GetHeader()->GetCentralityP();
 919       else if (fESD)
 920         centralityObj = fESD->GetCentrality();
 921
 922       if (centralityObj)
 923         centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
 924         //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
 925       else
 926         centrality = -1;
 927
 928       if (fAOD)
 929       {
 930         // remove outliers
 931         if (centrality == 0)
 932         {
 933           if (fAOD->GetVZEROData())
 934           {
 935             Float_t multV0 = 0;
 936             for (Int_t i=0; i<64; i++)
 937               multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
 938             if (multV0 < 19500)
 939             {
 940               centrality = -1;
 941               AliInfo("Rejecting event due to too small V0 multiplicity");
 942             }
 943           }
 944         }
 945       }
 946     }
 947
 948     AliInfo(Form("Centrality is %f", centrality));
 949   }
 950
 951   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
 952
 953   fHistos->SetRunNumber(inputEvent->GetRunNumber());
 954
 955   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 956   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
 957
 958   // Trigger selection ************************************************
 959   if (!fSkipTrigger && !fAnalyseUE->TriggerSelection(fInputHandler)) return;
 960
 961   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 962   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
 963
 964   // Vertex selection *************************************************
 965   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
 966
 967   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 968   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
 969
 970   // optimization
 971   if (centrality < 0 && !fCompareCentralities)
 972     return;
 973
 974   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesTrigger, kTRUE);
 975   //Printf("Accepted %d tracks", tracks->GetEntries());
 976
 977   // check for outlier in centrality vs number of tracks (rough constants extracted from correlation histgram)
 978   Bool_t reject = kFALSE;
 979   if (fRejectCentralityOutliers)
 980   {
 981     if (centrality > 40 && centrality <= 50 && tracks->GetEntriesFast() > 1160)
 982       reject = kTRUE;
 983     if (centrality > 50 && centrality <= 60 && tracks->GetEntriesFast() > 650)
 984       reject = kTRUE;
 985     if (centrality > 60 && centrality <= 70 && tracks->GetEntriesFast() > 370)
 986       reject = kTRUE;
 987     if (centrality > 70 && centrality <= 80 && tracks->GetEntriesFast() > 220)
 988       reject = kTRUE;
 989     if (centrality > 80 && centrality <= 90 && tracks->GetEntriesFast() > 130)
 990       reject = kTRUE;
 991     if (centrality > 90 && tracks->GetEntriesFast() > 75)
 992       reject = kTRUE;
 993   }
 994
 995   if (reject)
 996   {
 997     AliInfo(Form("Rejecting event due to centrality vs tracks correlation: %f %d", centrality, tracks->GetEntriesFast()));
 998     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
 999     delete tracks;
1000     return;
1001   }
1002
1003   // correlate particles with...
1004   TObjArray* tracksCorrelate = 0;
1005   if (fParticleSpeciesAssociated != fParticleSpeciesTrigger)
1006     tracksCorrelate = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, fParticleSpeciesAssociated, kTRUE);
1007
1008   // reference multiplicity
1009   Int_t referenceMultiplicity = -1;
1010   if (fESD)
1011     referenceMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity(fESD);
1012   else if (fAOD)
1013     referenceMultiplicity = tracks->GetEntriesFast(); // TODO to be replaced by the estimator once available in the AOD
1014 //    referenceMultiplicity = fAOD->GetHeader()->GetRefMultiplicityComb05();
1015
1016   ((TH2F*) fListOfHistos->FindObject("referenceMultiplicity"))->Fill(centrality, referenceMultiplicity);
1017
1018   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
1019   Double_t zVtx = vertex->GetZ();
1020
1021   Float_t weight = 1;
1022   if (fFillpT)
1023     weight = -1;
1024
1025   // fill two different centralities (syst study)
1026   // the zvtx axis is used to distinguish the results of both centralities: configured centrality in zvtx = 0, SPD in zvtx = 2
1027   if (fCompareCentralities)
1028     zVtx = 0;
1029
1030   // Fill containers at STEP 6 (reconstructed)
1031   if (centrality >= 0)
1032   {
1033     if (!fSkipStep6)
1034       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
1035
1036     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
1037
1038     if (fTwoTrackEfficiencyCut > 0)
1039       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, tracksCorrelate, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1040   }
1041
1042   // fill second time with SPD centrality
1043   if (fCompareCentralities && centralityObj)
1044   {
1045     centrality = centralityObj->GetCentralityPercentile("CL1");
1046     if (centrality >= 0 && !fSkipStep6)
1047       fHistos->FillCorrelations(centrality, 2, AliUEHist::kCFStepReconstructed, tracks, tracksCorrelate, weight, kFALSE, kFALSE, 0, 0.02, kTRUE);
1048   }
1049
1050   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
1051   TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
1052   delete tracks;
1053
1054   if (fFillMixed)
1055   {
1056     // event mixing
1057
1058     // 1. First get an event pool corresponding in mult (cent) and
1059     //    zvertex to the current event. Once initialized, the pool
1060     //    should contain nMix (reduced) events. This routine does not
1061     //    pre-scan the chain. The first several events of every chain
1062     //    will be skipped until the needed pools are filled to the
1063     //    specified depth. If the pool categories are not too rare, this
1064     //    should not be a problem. If they are rare, you could lose
1065     //    statistics.
1066
1067     // 2. Collect the whole pool's content of tracks into one TObjArray
1068     //    (bgTracks), which is effectively a single background super-event.
1069
1070     // 3. The reduced and bgTracks arrays must both be passed into
1071     //    FillCorrelations(). Also nMix should be passed in, so a weight
1072     //    of 1./nMix can be applied.
1073
1074     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
1075
1076     if (!pool)
1077       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
1078
1079 //     pool->SetDebug(1);
1080
1081     if (pool->IsReady())
1082     {
1083
1084       Int_t nMix = pool->GetCurrentNEvents();
1085 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
1086
1087       ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
1088       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
1089       if (pool->IsReady())
1090         ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
1091
1092       // Fill mixed-event histos here
1093       for (Int_t jMix=0; jMix<nMix; jMix++)
1094       {
1095         TObjArray* bgTracks = pool->GetEvent(jMix);
1096
1097         if (!fSkipStep6)
1098           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kFALSE, 0, 0.02, kTRUE);
1099
1100         if (fTwoTrackEfficiencyCut > 0)
1101           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
1102       }
1103     }
1104
1105     // ownership is with the pool now
1106     if (tracksCorrelate)
1107     {
1108       pool->UpdatePool(CloneAndReduceTrackList(tracksCorrelate));
1109       delete tracksClone;
1110     }
1111     else
1112       pool->UpdatePool(tracksClone);
1113     //pool->PrintInfo();
1114   }
1115   else
1116   {
1117     delete tracksClone;
1118     if (tracksCorrelate)
1119       delete tracksCorrelate;
1120   }
1121 }
1122
1123 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
1124 {
1125   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
1126
1127   TObjArray* tracksClone = new TObjArray;
1128   tracksClone->SetOwner(kTRUE);
1129
1130   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
1131   {
1132     AliVParticle* particle = (AliVParticle*) tracks->At(i);
1133     tracksClone->Add(new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge()));
1134   }
1135
1136   return tracksClone;
1137 }
1138
1139 //____________________________________________________________________
1140 void  AliAnalysisTaskPhiCorrelations::Initialize()
1141 {
1142   // input handler
1143   fInputHandler = (AliInputEventHandler*)
1144          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1145   // MC handler
1146   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1147 }
1148
1149 //____________________________________________________________________
1150 void AliAnalysisTaskPhiCorrelations::RemoveDuplicates(TObjArray* tracks)
1151 {
1152   // remove particles with the same label
1153
1154   Int_t before = tracks->GetEntriesFast();
1155
1156   for (Int_t i=0; i<before; ++i)
1157   {
1158     AliVParticle* part = (AliVParticle*) tracks->At(i);
1159
1160     for (Int_t j=i+1; j<before; ++j)
1161     {
1162       AliVParticle* part2 = (AliVParticle*) tracks->At(j);
1163
1164       if (part->GetLabel() == part2->GetLabel())
1165       {
1166         Printf("Removing %d with label %d (duplicated in %d)", i, part->GetLabel(), j); part->Dump(); part2->Dump();
1167         TObject* object = tracks->RemoveAt(i);
1168         if (tracks->IsOwner())
1169           delete object;
1170         break;
1171       }
1172     }
1173   }
1174
1175   tracks->Compress();
1176
1177   if (before > tracks->GetEntriesFast())
1178     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
1179 }
1180
1181 void AliAnalysisTaskPhiCorrelations::CleanUp(TObjArray* tracks, TObject* mcObj, Int_t maxLabel)
1182 {
1183   // calls RemoveInjectedSignals, RemoveWeakDecays and RemoveDuplicates
1184
1185   if (!tracks)
1186     return;
1187
1188   if (fInjectedSignals)
1189     fAnalyseUE->RemoveInjectedSignals(tracks, mcObj, maxLabel);
1190   if (fRemoveWeakDecays)
1191     fAnalyseUE->RemoveWeakDecays(tracks, mcObj);
1192   if (fRemoveDuplicates)
1193     RemoveDuplicates(tracks);
1194 }
1195
1196 //____________________________________________________________________
1197 void AliAnalysisTaskPhiCorrelations::SelectCharge(TObjArray* tracks)
1198 {
1199   // remove particles with charge not selected (depending on fTriggerSelectCharge)
1200
1201   if (!tracks)
1202     return;
1203
1204   Int_t before = tracks->GetEntriesFast();
1205
1206   for (Int_t i=0; i<before; ++i)
1207   {
1208     AliVParticle* part = (AliVParticle*) tracks->At(i);
1209
1210     if (part->Charge() * fTriggerSelectCharge < -1)
1211     {
1212 //       Printf("Removing %d with charge %d", i, part->Charge());
1213       TObject* object = tracks->RemoveAt(i);
1214       if (tracks->IsOwner())
1215         delete object;
1216     }
1217   }
1218
1219   tracks->Compress();
1220
1221   if (before > tracks->GetEntriesFast())
1222     AliInfo(Form("Reduced from %d to %d", before, tracks->GetEntriesFast()));
1223 }