PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx

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