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