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