PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx

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