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