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