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   Int_t nZvtxBins  = 7+(1+7)*4;
 278   // bins for further buffers are shifted by 100 cm
 279   Double_t vertexBins[] = { -7, -5, -3, -1, 1, 3, 5, 7,
 280                             93, 95, 97, 99, 101, 103, 105, 107,
 281                             193, 195, 197, 199, 201, 203, 205, 207,
 282                             293, 295, 297, 299, 301, 303, 305, 307,
 283                             393, 395, 397, 399, 401, 403, 405, 407 };
 284   Double_t* zvtxbin = vertexBins;
 285
 286   if (fMode == 0 && fHistos->GetUEHist(2)->GetEventHist()->GetNVar() > 2)
 287   {
 288     nZvtxBins = fHistos->GetUEHist(2)->GetEventHist()->GetNBins(2);
 289     zvtxbin = (Double_t*) fHistos->GetUEHist(2)->GetEventHist()->GetAxis(2, 0)->GetXbins()->GetArray();
 290   }
 291
 292   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins, centralityBins, nZvtxBins, zvtxbin);
 293 }
 294
 295 //____________________________________________________________________
 296 void  AliAnalysisTaskPhiCorrelations::Exec(Option_t */*option*/)
 297 {
 298   // receive ESD pointer if we are not running AOD analysis
 299   if (!fAOD)
 300   {
 301     AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
 302     if (handler && handler->InheritsFrom("AliESDInputHandler"))
 303       fESD = ((AliESDInputHandler*)handler)->GetEvent();
 304   }
 305
 306   if (fMode)
 307   {
 308     // correction mode
 309
 310     if (fMcHandler)
 311       fMcEvent = fMcHandler->MCEvent();
 312
 313     if (fAOD)
 314     {
 315       // array of MC particles
 316       fArrayMC = dynamic_cast<TClonesArray*>(fAOD->FindListObject(AliAODMCParticle::StdBranchName()));
 317       if (!fArrayMC)
 318         AliFatal("No array of MC particles found !!!");
 319     }
 320
 321     AnalyseCorrectionMode();
 322   }
 323   else AnalyseDataMode();
 324 }
 325
 326 /******************** ANALYSIS METHODS *****************************/
 327
 328 //____________________________________________________________________
 329 void  AliAnalysisTaskPhiCorrelations::AddSettingsTree()
 330 {
 331   //Write settings to output list
 332   TTree *settingsTree   = new TTree("UEAnalysisSettings","Analysis Settings in UE estimation");
 333   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
 334   settingsTree->Branch("fZVertex", &fZVertex,"ZVertex/D");
 335   //settingsTree->Branch("fCentralityMethod", fCentralityMethod.Data(),"CentralityMethod/C");
 336   settingsTree->Branch("fTrackEtaCut", &fTrackEtaCut, "TrackEtaCut/D");
 337   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
 338   settingsTree->Branch("fPtMin", &fPtMin, "PtMin/D");
 339   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
 340   settingsTree->Branch("fSelectBit", &fSelectBit,"EventSelectionBit/I");
 341   settingsTree->Branch("fUseChargeHadrons", &fUseChargeHadrons,"UseChHadrons/O");
 342   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
 343   settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
 344   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
 345   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
 346   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
 347   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
 348   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
 349   settingsTree->Branch("fMixingTracks", &fMixingTracks,"MixingTracks/I");
 350   settingsTree->Branch("fSkipTrigger", &fSkipTrigger,"SkipTrigger/O");
 351   settingsTree->Branch("fInjectedSignals", &fInjectedSignals,"SkipTrigger/O");
 352   settingsTree->Branch("fRejectCentralityOutliers", &fRejectCentralityOutliers,"RejectCentralityOutliers/O");
 353
 354   settingsTree->Fill();
 355   fListOfHistos->Add(settingsTree);
 356 }
 357
 358 //____________________________________________________________________
 359 void  AliAnalysisTaskPhiCorrelations::AnalyseCorrectionMode()
 360 {
 361   // Run the analysis on MC to get the correction maps
 362   //
 363
 364   if ( fDebug > 3 ) AliInfo( " Processing event in Corrections mode ..." );
 365
 366   Double_t centrality = 0;
 367
 368   if (fCentralityMethod.Length() > 0)
 369   {
 370     AliCentrality *centralityObj = 0;
 371     if (fAOD)
 372       centralityObj = fAOD->GetHeader()->GetCentralityP();
 373     else if (fESD)
 374       centralityObj = fESD->GetCentrality();
 375
 376     if (centralityObj)
 377     {
 378       centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
 379       AliInfo(Form("Centrality is %f", centrality));
 380     }
 381     else
 382     {
 383       Printf("WARNING: Centrality object is 0");
 384       centrality = -1;
 385      }
 386   }
 387
 388   // Support for ESD and AOD based analysis
 389   AliVEvent* inputEvent = fAOD;
 390   if (!inputEvent)
 391     inputEvent = fESD;
 392
 393   Float_t bSign = 0;
 394
 395   if (inputEvent)
 396   {
 397     fHistos->SetRunNumber(inputEvent->GetRunNumber());
 398     bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
 399   }
 400
 401   TObject* mc = fArrayMC;
 402   if (!mc)
 403     mc = fMcEvent;
 404
 405   // count all events
 406   fHistos->FillEvent(centrality, -1);
 407
 408   if (centrality < 0)
 409     return;
 410
 411   // Only consider MC events within the vtx-z region used also as cut on the reconstructed vertex
 412   TObject* vertexSupplier = fMcEvent;
 413   if (fAOD) // AOD
 414     vertexSupplier = fAOD->FindListObject(AliAODMCHeader::StdBranchName());
 415
 416   if (!fAnalyseUE->VertexSelection(vertexSupplier, 0, fZVertex))
 417     return;
 418
 419   Float_t zVtx = 0;
 420   if (fAOD)
 421     zVtx = ((AliAODMCHeader*) vertexSupplier)->GetVtxZ();
 422   else
 423     zVtx = fMcEvent->GetPrimaryVertex()->GetZ();
 424   Float_t weight = 1;
 425   if (fFillpT)
 426     weight = -1;
 427
 428   // For productions with injected signals, figure out above which label to skip particles/tracks
 429   Int_t skipParticlesAbove = 0;
 430   if (fInjectedSignals)
 431   {
 432     AliGenEventHeader* eventHeader = 0;
 433     Int_t headers = 0;
 434
 435     if (fMcEvent)
 436     {
 437       // ESD
 438       AliHeader* header = (AliHeader*) fMcEvent->Header();
 439       if (!header)
 440         AliFatal("fInjectedSignals set but no MC header found");
 441
 442       AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
 443       if (!cocktailHeader)
 444       {
 445         header->Dump();
 446         AliFatal("fInjectedSignals set but no MC cocktail header found");
 447       }
 448
 449       headers = cocktailHeader->GetHeaders()->GetEntries();
 450       eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
 451     }
 452     else
 453     {
 454       // AOD
 455       AliAODMCHeader* header = (AliAODMCHeader*) fAOD->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 456       if (!header)
 457         AliFatal("fInjectedSignals set but no MC header found");
 458
 459       headers = header->GetNCocktailHeaders();
 460       eventHeader = header->GetCocktailHeader(0);
 461     }
 462
 463     if (!eventHeader)
 464     {
 465       // We avoid AliFatal here, because the AOD productions sometimes have events where the MC header is missing
 466       // (due to unreadable Kinematics) and we don't want to loose the whole job because of a few events
 467       AliError("First event header not found. Skipping this event.");
 468       fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
 469       return;
 470     }
 471
 472     skipParticlesAbove = eventHeader->NProduced();
 473     AliInfo(Form("Injected signals in this event (%d headers). Keeping events of %s. Will skip particles/tracks above %d.", headers, eventHeader->ClassName(), skipParticlesAbove));
 474   }
 475
 476   // Get MC primaries
 477   TObjArray* tmpList = fAnalyseUE->GetAcceptedParticles(mc, 0, kTRUE, -1, kTRUE);
 478   if (fInjectedSignals)
 479     fAnalyseUE->RemoveInjectedSignals(tmpList, mc, skipParticlesAbove);
 480   TObjArray* tracksMC = CloneAndReduceTrackList(tmpList);
 481   delete tmpList;
 482
 483   /*
 484   if (fAOD)
 485   {
 486     for (Int_t i=0; i<fArrayMC->GetEntriesFast(); i++)
 487       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(((AliAODMCParticle*) fArrayMC->At(i))->PdgCode());
 488   }
 489   else
 490   {
 491     for (Int_t i=0; i<fMcEvent->GetNumberOfTracks(); i++)
 492       ((TH1F*) fListOfHistos->FindObject("pids"))->Fill(fMcEvent->GetTrack(i)->PdgCode());
 493   }
 494   */
 495
 496   if (fFillOnlyStep0)
 497     zVtx = 0;
 498
 499   // (MC-true all particles)
 500   // STEP 0
 501   fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, 0, weight);
 502
 503   // mixed event
 504   if (fFillMixed)
 505   {
 506     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
 507     //pool->PrintInfo();
 508     if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
 509       for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 510         fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepAll, tracksMC, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 511     pool->UpdatePool(CloneAndReduceTrackList(tracksMC));
 512   }
 513
 514 //   Printf("trigger: %d", ((AliInputEventHandler*)fInputHandler)->IsEventSelected());
 515
 516   // Trigger selection ************************************************
 517   if (!fFillOnlyStep0 && (fSkipTrigger || fAnalyseUE->TriggerSelection(fInputHandler)))
 518   {
 519     // (MC-true all particles)
 520     // STEP 1
 521     if (!fReduceMemoryFootprint)
 522       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTriggered, tracksMC, 0, weight);
 523     else
 524       fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
 525
 526     if (!inputEvent) {
 527       AliFatal("UNEXPECTED: inputEvent is 0. Trigger selection should have failed");
 528       return;
 529     }
 530
 531     // Vertex selection *************************************************
 532     if (fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex))
 533     {
 534       // fill here for tracking efficiency
 535       // loop over particle species
 536
 537       for (Int_t particleSpecies = 0; particleSpecies < 4; particleSpecies++)
 538       {
 539         TObjArray* primMCParticles = fAnalyseUE->GetAcceptedParticles(mc, 0x0, kTRUE, particleSpecies, kTRUE);
 540         TObjArray* primRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, particleSpecies, kTRUE);
 541         TObjArray* allRecoTracksMatched = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, particleSpecies, kTRUE);
 542
 543         if (fInjectedSignals)
 544         {
 545           fAnalyseUE->RemoveInjectedSignals(primMCParticles, mc, skipParticlesAbove);
 546           fAnalyseUE->RemoveInjectedSignals(primRecoTracksMatched, mc, skipParticlesAbove);
 547           fAnalyseUE->RemoveInjectedSignals(allRecoTracksMatched, mc, skipParticlesAbove);
 548         }
 549
 550         fHistos->FillTrackingEfficiency(primMCParticles, primRecoTracksMatched, allRecoTracksMatched, 0, particleSpecies, centrality, zVtx);
 551
 552 //      Printf("%d --> %d %d %d", particleSpecies, primMCParticles->GetEntries(), primRecoTracksMatched->GetEntries(), allRecoTracksMatched->GetEntries());
 553
 554         delete primMCParticles;
 555         delete primRecoTracksMatched;
 556         delete allRecoTracksMatched;
 557       }
 558       TObjArray* fakeParticles = fAnalyseUE->GetFakeParticles(inputEvent, mc, kFALSE, -1, kTRUE);
 559       if (fInjectedSignals)
 560       {
 561         fAnalyseUE->RemoveInjectedSignals((TObjArray*) fakeParticles->At(0), mc, skipParticlesAbove);
 562         fAnalyseUE->RemoveInjectedSignals((TObjArray*) fakeParticles->At(1), mc, skipParticlesAbove);
 563       }
 564       fHistos->FillTrackingEfficiency(0, 0, 0, (TObjArray*) fakeParticles->At(2), -1, centrality, zVtx);
 565       fHistos->FillFakePt(fakeParticles, centrality);
 566 //       Printf(">>>>> %d %d %d fakes", ((TObjArray*) fakeParticles->At(0))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(1))->GetEntriesFast(), ((TObjArray*) fakeParticles->At(2))->GetEntriesFast());
 567       delete fakeParticles;
 568
 569       // (MC-true all particles)
 570       // STEP 2
 571       if (!fReduceMemoryFootprint)
 572         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepVertex, tracksMC, 0, weight);
 573       else
 574         fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
 575
 576       // Get MC primaries that match reconstructed track
 577       TObjArray* tracksRecoMatchedPrim = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kTRUE, -1, kTRUE);
 578       if (fInjectedSignals)
 579         fAnalyseUE->RemoveInjectedSignals(tracksRecoMatchedPrim, mc, skipParticlesAbove);
 580
 581       // (RECO-matched (quantities from MC particle) primary particles)
 582       // STEP 4
 583       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, 0, weight);
 584
 585       // mixed event
 586       if (fFillMixed)
 587       {
 588         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 200);
 589         //pool->PrintInfo();
 590         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
 591           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 592             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTrackedOnlyPrim, tracksRecoMatchedPrim, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 593         pool->UpdatePool(CloneAndReduceTrackList(tracksRecoMatchedPrim));
 594       }
 595
 596       // Get MC primaries + secondaries that match reconstructed track
 597       TObjArray* tracksRecoMatchedAll = fAnalyseUE->GetAcceptedParticles(inputEvent, mc, kFALSE, -1, kTRUE);
 598       if (fInjectedSignals)
 599         fAnalyseUE->RemoveInjectedSignals(tracksRecoMatchedAll, mc, skipParticlesAbove);
 600
 601       // (RECO-matched (quantities from MC particle) all particles)
 602       // STEP 5
 603       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, 0, weight);
 604
 605       // mixed event
 606       if (fFillMixed)
 607       {
 608         AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx + 300);
 609         //pool->PrintInfo();
 610         if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
 611           for (Int_t jMix=0; jMix<pool->GetCurrentNEvents(); jMix++)
 612             fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepTracked, tracksRecoMatchedAll, pool->GetEvent(jMix), 1.0 / pool->GetCurrentNEvents(), (jMix == 0));
 613         pool->UpdatePool(CloneAndReduceTrackList(tracksRecoMatchedAll));
 614       }
 615
 616       // Get RECO tracks
 617       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
 618       if (fInjectedSignals)
 619         fAnalyseUE->RemoveInjectedSignals(tracks, mc, skipParticlesAbove);
 620
 621       // (RECO all tracks)
 622       // STEP 6
 623       if (!fSkipStep6)
 624         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, 0, weight);
 625
 626       // two track cut, STEP 8
 627       if (fTwoTrackEfficiencyCut > 0)
 628         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, 0, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut);
 629
 630       // apply correction efficiency, STEP 10
 631       if (fEfficiencyCorrection)
 632       {
 633         // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
 634         Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
 635
 636         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, 0, weight, kTRUE, twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
 637       }
 638
 639       // mixed event
 640       if (fFillMixed)
 641       {
 642         AliEventPool* pool2 = fPoolMgr->GetEventPool(centrality, zVtx + 100);
 643         //pool2->PrintInfo();
 644         if (pool2->IsReady() || pool2->NTracksInPool() > fMixingTracks / 10 || pool2->GetCurrentNEvents() >= 5)
 645         {
 646           for (Int_t jMix=0; jMix<pool2->GetCurrentNEvents(); jMix++)
 647           {
 648             // STEP 6
 649             if (!fSkipStep6)
 650               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0));
 651
 652             // two track cut, STEP 8
 653             if (fTwoTrackEfficiencyCut > 0)
 654               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut);
 655
 656             // apply correction efficiency, STEP 10
 657             if (fEfficiencyCorrection)
 658             {
 659               // with or without two track efficiency depending on if fTwoTrackEfficiencyCut is set
 660               Bool_t twoTrackCut = (fTwoTrackEfficiencyCut > 0);
 661
 662               fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepCorrected, tracks, pool2->GetEvent(jMix), 1.0 / pool2->GetCurrentNEvents(), (jMix == 0), twoTrackCut, bSign, fTwoTrackEfficiencyCut, kTRUE);
 663             }
 664           }
 665         }
 666         pool2->UpdatePool(CloneAndReduceTrackList(tracks));
 667       }
 668
 669       if (0 && !fReduceMemoryFootprint)
 670       {
 671         // make list of secondaries (matched with MC)
 672         TObjArray* tracksRecoMatchedSecondaries = new TObjArray;
 673         for (Int_t i=0; i<tracksRecoMatchedAll->GetEntriesFast(); i++)
 674           if (((AliAODMCParticle*)tracksRecoMatchedAll->At(i))->IsPhysicalPrimary() == kFALSE)
 675             tracksRecoMatchedSecondaries->Add(tracksRecoMatchedAll->At(i));
 676
 677         // Study: Use only secondaries as trigger particles and plot the correlation vs. all particles; store in step 9
 678         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy2, tracksRecoMatchedSecondaries, tracksRecoMatchedAll, weight);
 679
 680         // Study: Use only primaries as trigger particles and plot the correlation vs. secondaries; store in step 8
 681         fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksRecoMatchedPrim, tracksRecoMatchedSecondaries, weight);
 682
 683         // plot delta phi vs process id of secondaries
 684         // trigger particles: primaries in 4 < pT < 10
 685         // associated particles: secondaries in 1 < pT < 10
 686
 687         for (Int_t i=0; i<tracksRecoMatchedPrim->GetEntriesFast(); i++)
 688         {
 689           AliVParticle* triggerParticle = (AliVParticle*) tracksRecoMatchedPrim->At(i);
 690
 691           if (triggerParticle->Pt() < 4 || triggerParticle->Pt() > 10)
 692             continue;
 693
 694           for (Int_t j=0; j<tracksRecoMatchedSecondaries->GetEntriesFast(); j++)
 695           {
 696             AliAODMCParticle* particle = (AliAODMCParticle*) tracksRecoMatchedSecondaries->At(j);
 697
 698             if (particle->Pt() < 1 || particle->Pt() > 10)
 699               continue;
 700
 701             if (particle->Pt() > triggerParticle->Pt())
 702               continue;
 703
 704             Double_t deltaPhi = triggerParticle->Phi() - particle->Phi();
 705             if (deltaPhi > 1.5 * TMath::Pi())
 706               deltaPhi -= TMath::TwoPi();
 707             if (deltaPhi < -0.5 * TMath::Pi())
 708               deltaPhi += TMath::TwoPi();
 709
 710             Int_t processID = fMcEvent->Stack()->Particle(particle->GetLabel())->GetUniqueID();
 711
 712             ((TH2F*) fListOfHistos->FindObject("processIDs"))->Fill(deltaPhi, processID);
 713           }
 714         }
 715
 716         delete tracksRecoMatchedSecondaries;
 717       }
 718
 719       delete tracksRecoMatchedPrim;
 720       delete tracksRecoMatchedAll;
 721       delete tracks;
 722     }
 723   }
 724
 725   delete tracksMC;
 726 }
 727
 728 //____________________________________________________________________
 729 void  AliAnalysisTaskPhiCorrelations::AnalyseDataMode()
 730 {
 731
 732   // Run the analysis on DATA or MC to get raw distributions
 733
 734   if ( fDebug > 3 ) AliInfo( " Processing event in Data mode ..." );
 735
 736   ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(0);
 737
 738   if (!fInputHandler)
 739     return;
 740
 741   // skip not selected events here (the AOD is not updated for those)
 742   if (!fSkipTrigger && !(fInputHandler->IsEventSelected() & fSelectBit))
 743     return;
 744
 745   // Support for ESD and AOD based analysis
 746   AliVEvent* inputEvent = fAOD;
 747   if (!inputEvent)
 748     inputEvent = fESD;
 749
 750   Double_t centrality = 0;
 751
 752   AliCentrality *centralityObj = 0;
 753   if (fCentralityMethod.Length() > 0)
 754   {
 755     if (fCentralityMethod == "ZNA_MANUAL")
 756     {
 757       Bool_t zna = kFALSE;
 758       for(Int_t j = 0; j < 4; ++j) {
 759         if (fESD->GetZDCData()->GetZDCTDCData(12,j) != 0) {
 760           zna = kTRUE;
 761         }
 762       }
 763
 764 //       Printf("%d %f", zna, fZNAtower[0]);
 765       if (zna)
 766       {
 767         // code from Chiara O (23.10.12)
 768         const Double_t *fZNAtower = fESD->GetZDCData()->GetZN2TowerEnergy();
 769         Float_t znacut[4] = {681., 563., 413., 191.};
 770
 771         if(fZNAtower[0]>znacut[0]) centrality = 1;
 772         else if(fZNAtower[0]>znacut[1]) centrality = 21;
 773         else if(fZNAtower[0]>znacut[2]) centrality = 41;
 774         else if(fZNAtower[0]>znacut[3]) centrality = 61;
 775         else centrality = 81;
 776       }
 777       else
 778         centrality = -1;
 779     }
 780     else if (fCentralityMethod == "TRACKS_MANUAL")
 781     {
 782       // for pp
 783       TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
 784       centrality = tracks->GetEntriesFast();
 785       if (centrality > 40)
 786         centrality = 41;
 787 //       Printf("%d %f", tracks->GetEntriesFast(), centrality);
 788       delete tracks;
 789     }
 790     else
 791     {
 792       if (fAOD)
 793         centralityObj = fAOD->GetHeader()->GetCentralityP();
 794       else if (fESD)
 795         centralityObj = fESD->GetCentrality();
 796
 797       if (centralityObj)
 798         centrality = centralityObj->GetCentralityPercentile(fCentralityMethod);
 799         //centrality = centralityObj->GetCentralityPercentileUnchecked(fCentralityMethod);
 800       else
 801         centrality = -1;
 802
 803       if (fAOD)
 804       {
 805         // remove outliers
 806         if (centrality == 0)
 807         {
 808           if (fAOD->GetVZEROData())
 809           {
 810             Float_t multV0 = 0;
 811             for (Int_t i=0; i<64; i++)
 812               multV0 += fAOD->GetVZEROData()->GetMultiplicity(i);
 813             if (multV0 < 19500)
 814             {
 815               centrality = -1;
 816               AliInfo("Rejecting event due to too small V0 multiplicity");
 817             }
 818           }
 819         }
 820       }
 821     }
 822
 823     AliInfo(Form("Centrality is %f", centrality));
 824   }
 825
 826   Float_t bSign = (inputEvent->GetMagneticField() > 0) ? 1 : -1;
 827
 828   fHistos->SetRunNumber(inputEvent->GetRunNumber());
 829
 830   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 831   fHistos->FillEvent(centrality, AliUEHist::kCFStepAll);
 832
 833   // Trigger selection ************************************************
 834   if (!fSkipTrigger && !fAnalyseUE->TriggerSelection(fInputHandler)) return;
 835
 836   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 837   fHistos->FillEvent(centrality, AliUEHist::kCFStepTriggered);
 838
 839   // Vertex selection *************************************************
 840   if(!fAnalyseUE->VertexSelection(inputEvent, fnTracksVertex, fZVertex)) return;
 841
 842   // Fill the "event-counting-container", it is needed to get the number of events remaining after each event-selection cut
 843   fHistos->FillEvent(centrality, AliUEHist::kCFStepVertex);
 844
 845   // optimization
 846   if (centrality < 0 && !fCompareCentralities)
 847     return;
 848
 849   TObjArray* tracks = fAnalyseUE->GetAcceptedParticles(inputEvent, 0, kTRUE, -1, kTRUE);
 850
 851   // check for outlier in centrality vs number of tracks (rough constants extracted from correlation histgram)
 852   Bool_t reject = kFALSE;
 853   if (fRejectCentralityOutliers)
 854   {
 855     if (centrality > 40 && centrality <= 50 && tracks->GetEntriesFast() > 1160)
 856       reject = kTRUE;
 857     if (centrality > 50 && centrality <= 60 && tracks->GetEntriesFast() > 650)
 858       reject = kTRUE;
 859     if (centrality > 60 && centrality <= 70 && tracks->GetEntriesFast() > 370)
 860       reject = kTRUE;
 861     if (centrality > 70 && centrality <= 80 && tracks->GetEntriesFast() > 220)
 862       reject = kTRUE;
 863     if (centrality > 80 && centrality <= 90 && tracks->GetEntriesFast() > 130)
 864       reject = kTRUE;
 865     if (centrality > 90 && tracks->GetEntriesFast() > 75)
 866       reject = kTRUE;
 867   }
 868
 869   if (reject)
 870   {
 871     AliInfo(Form("Rejecting event due to centrality vs tracks correlation: %f %d", centrality, tracks->GetEntriesFast()));
 872     fHistos->FillEvent(centrality, AliUEHist::kCFStepAnaTopology);
 873     delete tracks;
 874     return;
 875   }
 876
 877   // reference multiplicity
 878   Int_t referenceMultiplicity = -1;
 879   if (fESD)
 880     referenceMultiplicity = AliESDtrackCuts::GetReferenceMultiplicity(fESD);
 881
 882   ((TH2F*) fListOfHistos->FindObject("referenceMultiplicity"))->Fill(centrality, referenceMultiplicity);
 883
 884   // create a list of reduced objects. This speeds up processing and reduces memory consumption for the event pool
 885   TObjArray* tracksClone = CloneAndReduceTrackList(tracks);
 886   delete tracks;
 887
 888   //Printf("Accepted %d tracks", tracks->GetEntries());
 889
 890   const AliVVertex* vertex = inputEvent->GetPrimaryVertex();
 891   Double_t zVtx = vertex->GetZ();
 892
 893   Float_t weight = 1;
 894   if (fFillpT)
 895     weight = -1;
 896
 897   // fill two different centralities (syst study)
 898   // the zvtx axis is used to distinguish the results of both centralities: configured centrality in zvtx = 0, SPD in zvtx = 2
 899   if (fCompareCentralities)
 900     zVtx = 0;
 901
 902   // Fill containers at STEP 6 (reconstructed)
 903   if (centrality >= 0)
 904   {
 905     if (!fSkipStep6)
 906       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
 907
 908     ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(1);
 909
 910     if (fTwoTrackEfficiencyCut > 0)
 911       fHistos->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, 0, weight, kTRUE, kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
 912   }
 913
 914   // fill second time with SPD centrality
 915   if (fCompareCentralities && centralityObj)
 916   {
 917     centrality = centralityObj->GetCentralityPercentile("CL1");
 918     if (centrality >= 0 && !fSkipStep6)
 919       fHistos->FillCorrelations(centrality, 2, AliUEHist::kCFStepReconstructed, tracksClone, 0, weight, kTRUE, kFALSE, 0, 0.02, kTRUE);
 920   }
 921
 922   if (fFillMixed)
 923   {
 924     // event mixing
 925
 926     // 1. First get an event pool corresponding in mult (cent) and
 927     //    zvertex to the current event. Once initialized, the pool
 928     //    should contain nMix (reduced) events. This routine does not
 929     //    pre-scan the chain. The first several events of every chain
 930     //    will be skipped until the needed pools are filled to the
 931     //    specified depth. If the pool categories are not too rare, this
 932     //    should not be a problem. If they are rare, you could lose
 933     //    statistics.
 934
 935     // 2. Collect the whole pool's content of tracks into one TObjArray
 936     //    (bgTracks), which is effectively a single background super-event.
 937
 938     // 3. The reduced and bgTracks arrays must both be passed into
 939     //    FillCorrelations(). Also nMix should be passed in, so a weight
 940     //    of 1./nMix can be applied.
 941
 942     AliEventPool* pool = fPoolMgr->GetEventPool(centrality, zVtx);
 943
 944     if (!pool)
 945       AliFatal(Form("No pool found for centrality = %f, zVtx = %f", centrality, zVtx));
 946
 947     //pool->SetDebug(1);
 948
 949     if (pool->IsReady() || pool->NTracksInPool() > fMixingTracks / 10 || pool->GetCurrentNEvents() >= 5)
 950     {
 951
 952       Int_t nMix = pool->GetCurrentNEvents();
 953 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
 954
 955       ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(2);
 956       ((TH2F*) fListOfHistos->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
 957       if (pool->IsReady())
 958         ((TH1F*) fListOfHistos->FindObject("eventStat"))->Fill(3);
 959
 960       // Fill mixed-event histos here
 961       for (Int_t jMix=0; jMix<nMix; jMix++)
 962       {
 963         TObjArray* bgTracks = pool->GetEvent(jMix);
 964
 965         if (!fSkipStep6)
 966           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepReconstructed, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kFALSE, 0, 0.02, kTRUE);
 967
 968         if (fTwoTrackEfficiencyCut > 0)
 969           fHistosMixed->FillCorrelations(centrality, zVtx, AliUEHist::kCFStepBiasStudy, tracksClone, bgTracks, 1.0 / nMix, (jMix == 0), kTRUE, bSign, fTwoTrackEfficiencyCut, kTRUE);
 970       }
 971     }
 972
 973     // ownership is with the pool now
 974     pool->UpdatePool(tracksClone);
 975     //pool->PrintInfo();
 976   }
 977   else
 978     delete tracksClone;
 979 }
 980
 981 TObjArray* AliAnalysisTaskPhiCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
 982 {
 983   // clones a track list by using AliDPhiBasicParticle which uses much less memory (used for event mixing)
 984
 985   TObjArray* tracksClone = new TObjArray;
 986   tracksClone->SetOwner(kTRUE);
 987
 988   for (Int_t i=0; i<tracks->GetEntriesFast(); i++)
 989   {
 990     AliVParticle* particle = (AliVParticle*) tracks->At(i);
 991     tracksClone->Add(new AliDPhiBasicParticle(particle->Eta(), particle->Phi(), particle->Pt(), particle->Charge()));
 992   }
 993
 994   return tracksClone;
 995 }
 996
 997 //____________________________________________________________________
 998 void  AliAnalysisTaskPhiCorrelations::Initialize()
 999 {
1000   // input handler
1001   fInputHandler = (AliInputEventHandler*)
1002          ((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
1003   // MC handler
1004   fMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
1005 }