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