PWGCF/Correlations/DPhi/AliAnalysisTaskPhiCorrelations.cxx

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