PWGCF/Correlations/DPhi/PidPid/AliAnalysisTaskPidPidCorrelations.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 /*  AliAnalysisTaskPidPidCorrelations.cxx
  19  *
  20  *  Analysis task for PID-PID two-particle correlations using PID form  ITS, TPC, TOF, HMPID
  21  *
  22  *  Authors:    Gyula BENCEDI (WIGNER RCP) gyula.bencedi@cern.ch
  23  *              Levente MOLNAR (CNRS-IPHC) levente.molnar@cern.ch
  24  *
  25  *  NOTE: running only on AODs, in PbPb
  26  *
  27  *
  28  *  Last modified: Wed Jan 29 15:08:56 CET 2014
  29  *
  30  *
  31 */
  32
  33 //_____ ROOT fMyAODHeaders
  34 #include "Riostream.h"
  35 #include "TChain.h"
  36 #include "TTree.h"
  37 #include "TAxis.h"
  38 #include "TF1.h"
  39 #include "TH1F.h"
  40 #include "TH2F.h"
  41 #include "TH3D.h"
  42 #include "THnSparse.h"
  43 #include "THashList.h"
  44 #include "TCanvas.h"
  45 #include "TROOT.h"
  46 #include "TList.h"
  47 #include "TFormula.h"
  48 #include "TObject.h"
  49 #include "TObjArray.h"
  50 #include "TParticle.h"
  51 #include "TPDGCode.h"
  52 #include "TParticlePDG.h"
  53 #include "TDatabasePDG.h"
  54
  55 #include <vector>
  56 #include <algorithm>
  57 #include <functional>
  58
  59 //_____ ALIROOT headers
  60 #include "AliAnalysisFilter.h"
  61
  62 //_____  Manager, Handler
  63 #include "AliAnalysisManager.h"
  64 #include "AliInputEventHandler.h"
  65
  66 //_____ Event pool includes
  67 #include "AliEventPoolManager.h"
  68 // #include "AliPool.h"
  69
  70 //_____ AOD includes
  71 #include "AliAODEvent.h"
  72 #include "AliAODTrack.h"
  73 #include "AliAODHMPIDrings.h"
  74 #include "AliAODHandler.h"
  75 #include "AliAODVertex.h"
  76 #include "AliAODInputHandler.h"
  77 #include "AliAODpidUtil.h"
  78 #include "AliAODMCParticle.h"
  79 #include "AliAODMCHeader.h"
  80
  81 //_____ PID includes
  82 #include "AliAODPid.h"
  83 #include "AliPID.h"
  84 #include "AliPIDResponse.h"
  85 #include "AliTPCPIDResponse.h"
  86 #include "AliPIDCombined.h"
  87
  88 //_____ Additional includes
  89 #include "AliLog.h"
  90 #include "AliStack.h"
  91 #include "AliVEvent.h"
  92 #include "AliMultiplicity.h"
  93 #include "AliCentrality.h"
  94 #include "AliGenEventHeader.h"
  95 #include "AliTHn.h"
  96 #include "AliAnalysisUtils.h"
  97
  98 //_____ CORRFW includes
  99 #include "AliCFContainer.h"
 100 // #include "AliCFGridSparse.h"
 101 // #include "AliCFEffGrid.h"
 102 // #include "AliCFManager.h"
 103 // #include "AliCFCutBase.h"
 104
 105 //_____ AnalysisTask headers
 106 #include "AliAnalysisTaskSE.h"
 107 #include "AliAnalysisTaskPidPidCorrelations.h"
 108
 109 using namespace myAliPID;
 110 using namespace std;
 111
 112 ClassImp( AliAnalysisTaskPidPidCorrelations )
 113 ClassImp( AliPidPidCorrelationReducedTrack )
 114
 115 #define sign(x) = (x>0) ? 1 : -1
 116
 117 #define length(a) ( sizeof(a) / sizeof(*a) )
 118
 119 #define pi TMath::Pi()
 120
 121 //________________________________________________________________________
 122 AliAnalysisTaskPidPidCorrelations::AliAnalysisTaskPidPidCorrelations() // All data members should be initialised here
 123   : AliAnalysisTaskSE()
 124
 125   , fUseMC ( kFALSE )
 126
 127   , fMyAODEvent ( 0 )
 128   , fMyAODHeader ( 0 )
 129   , fMyAODTrack ( 0 )
 130   , fPIDResponse ( 0 )
 131   , fMyPrimVtx ( 0 )
 132   , fMyMcArray ( 0 )
 133   , fMyMcHeader ( 0 )
 134   , fMyMcHandler ( 0 )
 135   , fPoolMgr ( 0x0 )
 136   , fMyCFCont ( 0x0 )
 137   , fMyprimRecoTracksPID ( 0x0 )
 138   , fMyprimMCParticlesPID ( 0x0 )
 139   , fMyprimRecoTracksMatchedPID ( 0x0 )
 140
 141 //   , fWriteCorrTree ( kTRUE )
 142 //   , fVariablesTreeCorr ( 0x0 )
 143 //   , fCorrVariables ( 0 )
 144
 145   , fTriggerType ( 1 )
 146   , fMyMcType ( -1 )
 147 //   , fFilterBit ( 128 )
 148   , fFilterType ( 2 )
 149   , fCentrality ( -1 )
 150   , fCentralityPercentileMin ( 0. )
 151   , fCentralityPercentileMax ( 10. )
 152   , fNbinsCent ( 1 )
 153   , fCentAxis ( 0x0 )
 154   , fNbinsZvtx ( 1 )
 155   , fZvtxAxis ( 0x0 )
 156   , fNbinsPt ( 1 )
 157   , fPtAxis ( 0x0 )
 158   , fNbinsEta ( 1 )
 159   , fEtaAxis ( 0x0 )
 160   , fCentralityEstimator ( "V0M" )
 161
 162   , fTrackEtaMin ( -0.9 )
 163   , fTrackEtaMax ( 0.9 )
 164   , fTrackPtMin ( 0.2 )
 165   , fTrackPtMax ( 10.0 )
 166   , fTrackStatus ( 0 )
 167   , fnTracksVertex ( 1 )
 168   , fRejectZeroTrackEvents ( kTRUE )
 169   , fEtaCut ( 0.9 )
 170   , fVxMax ( 0.3 )
 171   , fVyMax ( 0.3 )
 172   , fVzMax ( 10. )
 173   , fRemoveWeakDecays ( kFALSE )
 174   , fRemoveDuplicates ( kFALSE )
 175   , fDeltaEtaMax ( 2. )
 176
 177   , fSelectCharge ( 0 )
 178   , fTriggerSelectCharge ( 0 )
 179   , fAssociatedSelectCharge ( 0 )
 180   , fTriggerRestrictEta ( -1 )
 181   , fCutConversions ( kFALSE )
 182   , fCutResonances ( kFALSE )
 183   , fRejectResonanceDaughters ( -1 )
 184   , fOnlyOneEtaSide ( 0 )
 185   , fWeightPerEvent ( kFALSE )
 186   , fPtOrder ( kTRUE )
 187   , fQCut ( kFALSE )
 188   , fDeltaPtMin ( 0. )
 189
 190   , fPIDtrigType ( 999 )
 191   , fPIDassocType ( 999 )
 192
 193   , fCustomBinning ( "" )
 194   , fBinningString ( "" )
 195
 196   , fMinNumTrack ( 2000 )
 197   , fPoolSize ( 1000 )
 198   , fMinNEventsToMix ( 5 )
 199
 200   , fFillpT ( kFALSE )
 201   , fTwoTrackEfficiencyCut ( 1 )
 202   , twoTrackEfficiencyCutValue ( 0.02 )
 203   , fTwoTrackCutMinRadius ( 0.8 )
 204   , fEtaOrdering ( kFALSE )
 205   , fFillMixed ( kTRUE )
 206
 207   , fList ( 0x0 )
 208   , fHistNev ( 0x0 )
 209   , fHistTriggerStats ( 0x0 )
 210   , fHistTriggerRun ( 0x0 )
 211   , fHistEventStat ( 0x0 )
 212   , fHistRefTracks ( 0x0 )
 213   , fHistCentStats ( 0x0 )
 214   , fHistCentralityPercentile ( 0x0 )
 215   , fHistCentralityClass10 ( 0x0 )
 216   , fHistCentralityClass5 ( 0x0 )
 217   , fHistV0M ( 0x0 )
 218   , fHistMcGenerator ( 0x0 )
 219   , fHist_HijingBg ( 0x0 )
 220   , fHistNumOfPartPerEvt ( 0x0 )
 221   , fHistMcStats ( 0x0 )
 222   , fHistMcAllPt ( 0x0 )
 223   , fHistMcAllPt_Hijing ( 0x0 )
 224   , fHistMcAllPt_Dec ( 0x0 )
 225   , fHistMcAllPt_Inj ( 0x0 )
 226   , fHistMcAllEta_NotHijing ( 0x0 )
 227   , fHistMcAllEta_Hijing ( 0x0 )
 228   , fHistMcAllEta ( 0x0 )
 229 //   , fHistCorrSingle ( 0x0 )
 230 //   , fHistPoolMgrQA ( 0x0 )
 231 //   , fHistCorrPairSame ( 0x0 )
 232 //   , fHistCorrPairMixed ( 0x0 )
 233   , fHistControlConvResoncances ( 0x0 )
 234
 235   , fHistTriggerWeighting ( 0x0 )
 236   , fTriggerWeighting ( 0x0 )
 237
 238   , fHistHBTbefore ( 0x0 )
 239   , fHistHBTafter ( 0x0 )
 240 //   , fHistConversionbefore ( 0x0 )
 241 //   , fHistConversionafter ( 0x0 )
 242 //   , fHistPsiMinusPhi ( 0x0 )
 243 //   , fHistResonancesBefore ( 0x0 )
 244 //   , fHistResonancesRho ( 0x0 )
 245 //   , fHistResonancesK0 ( 0x0 )
 246 //   , fHistResonancesLambda ( 0x0 )
 247 //   , fHistQbefore ( 0x0 )
 248 //   , fHistQafter ( 0x0 )
 249
 250 {
 251   // Dummy constructor ALWAYS needed for I/O.
 252 }
 253
 254 //________________________________________________________________________
 255 AliAnalysisTaskPidPidCorrelations::AliAnalysisTaskPidPidCorrelations(const char *name) // All data members should be initialised here
 256   : AliAnalysisTaskSE(name)
 257
 258   , fUseMC ( kFALSE )
 259
 260   , fMyAODEvent ( 0 )
 261   , fMyAODHeader ( 0 )
 262   , fMyAODTrack ( 0 )
 263   , fPIDResponse ( 0 )
 264   , fMyPrimVtx ( 0 )
 265   , fMyMcArray ( 0 )
 266   , fMyMcHeader ( 0 )
 267   , fMyMcHandler ( 0 )
 268   , fPoolMgr ( 0x0 )
 269   , fMyCFCont ( 0x0 )
 270   , fMyprimRecoTracksPID ( 0x0 )
 271   , fMyprimMCParticlesPID ( 0x0 )
 272   , fMyprimRecoTracksMatchedPID ( 0x0 )
 273
 274 //   , fWriteCorrTree ( kTRUE )
 275 //   , fVariablesTreeCorr ( 0x0 )
 276 //   , fCorrVariables ( 0 )
 277
 278   , fTriggerType ( 1 )
 279   , fMyMcType ( -1 )
 280 //   , fFilterBit ( 128 )
 281   , fFilterType ( 2 )
 282   , fCentrality ( -1 )
 283   , fCentralityPercentileMin ( 0. )
 284   , fCentralityPercentileMax ( 10. )
 285   , fNbinsCent ( 1 )
 286   , fCentAxis ( 0x0 )
 287   , fNbinsZvtx ( 1 )
 288   , fZvtxAxis ( 0x0 )
 289   , fNbinsPt ( 1 )
 290   , fPtAxis ( 0x0 )
 291   , fNbinsEta ( 1 )
 292   , fEtaAxis ( 0x0 )
 293   , fCentralityEstimator ( "V0M" )
 294
 295   , fTrackEtaMin ( -0.9 )
 296   , fTrackEtaMax ( 0.9 )
 297   , fTrackPtMin ( 0.2 )
 298   , fTrackPtMax ( 10.0 )
 299   , fTrackStatus ( 0 )
 300   , fnTracksVertex ( 1 )
 301   , fRejectZeroTrackEvents ( kTRUE )
 302   , fEtaCut ( 0.9 )
 303   , fVxMax ( 0.3 )
 304   , fVyMax ( 0.3 )
 305   , fVzMax ( 10. )
 306   , fRemoveWeakDecays ( kFALSE )
 307   , fRemoveDuplicates ( kFALSE )
 308   , fDeltaEtaMax ( 2. )
 309
 310   , fSelectCharge ( 0 )
 311   , fTriggerSelectCharge ( 0 )
 312   , fAssociatedSelectCharge ( 0 )
 313   , fTriggerRestrictEta ( -1 )
 314   , fCutConversions ( kFALSE )
 315   , fCutResonances ( kFALSE )
 316   , fRejectResonanceDaughters ( -1 )
 317   , fOnlyOneEtaSide ( 0 )
 318   , fWeightPerEvent ( kFALSE )
 319   , fPtOrder ( kTRUE )
 320   , fQCut ( kFALSE )
 321   , fDeltaPtMin ( 0. )
 322
 323   , fPIDtrigType ( 999 )
 324   , fPIDassocType ( 999 )
 325
 326   , fCustomBinning ( "" )
 327   , fBinningString ( "" )
 328
 329   , fMinNumTrack ( 2000 )
 330   , fPoolSize ( 1000 )
 331   , fMinNEventsToMix ( 5 )
 332
 333   , fFillpT ( kFALSE )
 334   , fTwoTrackEfficiencyCut ( 1 )
 335   , twoTrackEfficiencyCutValue ( 0.02 )
 336   , fTwoTrackCutMinRadius ( 0.8 )
 337   , fEtaOrdering ( kFALSE )
 338   , fFillMixed ( kTRUE )
 339
 340   , fList ( 0x0 )
 341   , fHistNev ( 0x0 )
 342   , fHistTriggerStats ( 0x0 )
 343   , fHistTriggerRun ( 0x0 )
 344   , fHistEventStat ( 0x0 )
 345   , fHistRefTracks ( 0x0 )
 346   , fHistCentStats ( 0x0 )
 347   , fHistCentralityPercentile ( 0x0 )
 348   , fHistCentralityClass10 ( 0x0 )
 349   , fHistCentralityClass5 ( 0x0 )
 350   , fHistV0M ( 0x0 )
 351   , fHistMcGenerator ( 0x0 )
 352   , fHist_HijingBg ( 0x0 )
 353   , fHistNumOfPartPerEvt ( 0x0 )
 354   , fHistMcStats ( 0x0 )
 355   , fHistMcAllPt ( 0x0 )
 356   , fHistMcAllPt_Hijing ( 0x0 )
 357   , fHistMcAllPt_Dec ( 0x0 )
 358   , fHistMcAllPt_Inj ( 0x0 )
 359   , fHistMcAllEta_NotHijing ( 0x0 )
 360   , fHistMcAllEta_Hijing ( 0x0 )
 361   , fHistMcAllEta ( 0x0 )
 362 //   , fHistCorrSingle ( 0x0 )
 363 //   , fHistPoolMgrQA ( 0x0 )
 364 //   , fHistCorrPairSame ( 0x0 )
 365 //   , fHistCorrPairMixed ( 0x0 )
 366   , fHistControlConvResoncances ( 0x0 )
 367
 368   , fHistTriggerWeighting ( 0x0 )
 369   , fTriggerWeighting ( 0x0 )
 370
 371   , fHistHBTbefore ( 0x0 )
 372   , fHistHBTafter ( 0x0 )
 373 //   , fHistConversionbefore ( 0x0 )
 374 //   , fHistConversionafter ( 0x0 )
 375 //   , fHistPsiMinusPhi ( 0x0 )
 376 //   , fHistResonancesBefore ( 0x0 )
 377 //   , fHistResonancesRho ( 0x0 )
 378 //   , fHistResonancesK0 ( 0x0 )
 379 //   , fHistResonancesLambda ( 0x0 )
 380 //   , fHistQbefore ( 0x0 )
 381 //   , fHistQafter ( 0x0 )
 382
 383 {
 384   // Constructor
 385
 386   //____ THnSparse
 387   for(Int_t i=0; i<2; i++) { fHistCorrPair[i] = 0x0; }
 388
 389   for (Int_t i=0; i<fNMaxBinsPt; i++)
 390     for(Int_t j=0; j<2; j++)
 391       fHistTwoTrackDistancePt[i][j] = 0x0;
 392
 393   for (Int_t iSpecies=0; iSpecies<AliPID::kSPECIES; iSpecies++)
 394   {
 395     nsigmaITS[iSpecies] = -999.;
 396     nsigmaTPC[iSpecies] = -999.;
 397     nsigmaTOF[iSpecies] = -999.;
 398     nsigmaHMPID[iSpecies] = -999.;
 399
 400     for (Int_t i=0; i<4; i++)
 401       fnsigmas[iSpecies][i] = 0x0;
 402   }
 403
 404   for (Int_t iCent=0; iCent<fNMaxBinsCentrality; iCent++)
 405   {
 406     //_____ QA
 407     for (Int_t i=0; i < 14; i++) fHistQA[iCent][i] = 0x0;
 408     for (Int_t i=0; i<6; i++) fHistRefTracksCent[iCent][i] = 0x0;
 409
 410     //_____ PID
 411     for (Int_t iSpecies=0; iSpecies<AliPID::kSPECIES; iSpecies++)
 412     {
 413       fHistNSigmaTPC[iCent][iSpecies] = 0x0;
 414       fHistNSigmaTOF[iCent][iSpecies] = 0x0;
 415     }
 416
 417     fHistTPCdEdx[iCent]         = 0x0;
 418     fHistTOFbeta[iCent]         = 0x0;
 419     fHistTPCdEdx_selected[iCent]        = 0x0;
 420     fHistTOFbeta_selected[iCent]        = 0x0;
 421
 422     fHistoNSigma[iCent] = 0x0;
 423
 424 //     for (Int_t iPtBinTrig=0; iPtBinTrig<fNMaxBinsPt; iPtBinTrig++)
 425 //     {
 426 //       for (Int_t iPtBinAssoc=0; iPtBinAssoc<fNMaxBinsPt; iPtBinAssoc++)
 427 //       {
 428 //      fHistDeltaPhi[iCent][iPtBinTrig][iPtBinAssoc]   = 0x0;
 429 //      fHistDeltaPhiMix[iCent][iPtBinTrig][iPtBinAssoc]        = 0x0;
 430 //      fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc]   = 0x0;
 431 //      fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc]        = 0x0;
 432 //       }
 433 //     }
 434     fHistSingleHadronsPt[iCent]         = 0x0;
 435 //     fHistSingleHadronsPt_Mixed[iCent]        = 0x0;
 436     fHistSingleHadronsEtaPt[iCent]      = 0x0;
 437 //     fHistSingleHadronsEtaPt_Mixed[iCent]     = 0x0;
 438   }
 439
 440   DefineOutput(1,TList::Class());
 441   DefineOutput(2,AliCFContainer::Class());
 442   DefineOutput(3,TTree::Class());
 443 }
 444
 445 //________________________________________________________________________
 446 AliAnalysisTaskPidPidCorrelations::~AliAnalysisTaskPidPidCorrelations()
 447 {
 448   // Destructor. Clean-up the output list, but not the histograms that are put inside
 449   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 450
 451   if (fCentAxis)        { delete fCentAxis; fCentAxis = 0x0; }
 452   if (fZvtxAxis)        { delete fZvtxAxis; fZvtxAxis = 0x0; }
 453   if (fPtAxis)  { delete fPtAxis; fPtAxis = 0x0; }
 454   if (fEtaAxis)         { delete fEtaAxis; fEtaAxis = 0x0; }
 455
 456   if (fList && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fList; fList = 0x0; }
 457   if (fMyCFCont && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { delete fMyCFCont; fMyCFCont = 0x0;  }
 458   if (fPoolMgr) { delete fPoolMgr; fPoolMgr = 0x0; }
 459   if (fPIDResponse) { delete fPIDResponse; fPIDResponse = 0x0; }
 460 }
 461
 462 //_______________________________________________________________________________
 463 void AliAnalysisTaskPidPidCorrelations::UserExec(Option_t *)
 464 {
 465   // _____ get input event
 466   AliVEvent* event = InputEvent();
 467   if (!event) { AliError("UserExec(). Could not retrieve event"); return; }
 468
 469   static int iEvent = 0; iEvent++;
 470   fHistNev -> SetBinContent(1,iEvent);
 471
 472   // _____ get AOD event
 473   fMyAODEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 474   if (!fMyAODEvent) { AliError("Cannot get the AOD event"); return; }
 475   // _____ get AOD fMyAODHeader
 476   fMyAODHeader = dynamic_cast<AliAODHeader*>(fMyAODEvent->GetHeader());
 477   if(!fMyAODHeader) { AliError("Cannot get the AOD fMyAODHeader"); return; }
 478
 479   if (fUseMC)
 480   {
 481     fMyMcArray = dynamic_cast<TClonesArray*>(fMyAODEvent->FindListObject(AliAODMCParticle::StdBranchName()));
 482     if (!fMyMcArray) { AliError("AliAnalysisTaskPidPidCorrelations::UserExec: Could not find Monte-Carlo in AOD"); return; }
 483     fMyMcHeader = (AliAODMCHeader*)fMyAODEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName());
 484     if (!fMyMcHeader) { AliError("AliAnalysisTaskPidPidCorrelations::UserExec: MC header branch not found!\n"); return; }
 485     // MC handler
 486     fMyMcHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 487   }
 488
 489   //____________ get correlation and correction
 490   Analyse();
 491
 492   PostData(1, fList);
 493   PostData(2, fMyCFCont);
 494 //   if (fWriteCorrTree) PostData(3,fVariablesTreeCorr);
 495 }
 496
 497 //____________________________________________________________________
 498 void  AliAnalysisTaskPidPidCorrelations::AddSettingsTree()
 499 {
 500   //Write settings to output list
 501   TTree *settingsTree   = new TTree("PIDAnalysisSettings","Analysis Settings in PID correlation");
 502   settingsTree->Branch("fnTracksVertex", &fnTracksVertex,"nTracksVertex/I");
 503 //   settingsTree->Branch("fFilterBit", &fFilterBit,"FilterBit/I");
 504   settingsTree->Branch("fTrackStatus", &fTrackStatus,"TrackStatus/I");
 505   settingsTree->Branch("fRejectZeroTrackEvents", &fRejectZeroTrackEvents,"RejectZeroTrackEvents/O");
 506   settingsTree->Branch("fRemoveWeakDecays", &fRemoveWeakDecays,"RemoveWeakDecays/O");
 507   settingsTree->Branch("fRemoveDuplicates", &fRemoveDuplicates,"RemoveDuplicates/O");
 508   settingsTree->Branch("fVzMax", &fVzMax,"ZVertex/D");
 509   settingsTree->Branch("fTrackEtaMin", &fTrackEtaMin, "fTrackEtaMin/D");
 510   settingsTree->Branch("fOnlyOneEtaSide", &fOnlyOneEtaSide,"OnlyOneEtaSide/I");
 511   settingsTree->Branch("fTrackPtMin", &fTrackPtMin, "fTrackPtMin/D");
 512   settingsTree->Branch("fTriggerType", &fTriggerType,"fTriggerType/I");
 513   settingsTree->Branch("fSelectCharge", &fSelectCharge,"SelectCharge/I");
 514   settingsTree->Branch("fTriggerSelectCharge", &fTriggerSelectCharge,"TriggerSelectCharge/I");
 515   settingsTree->Branch("fAssociatedSelectCharge", &fAssociatedSelectCharge,"fAssociatedSelectCharge/I");
 516   settingsTree->Branch("fTriggerRestrictEta", &fTriggerRestrictEta,"TriggerRestrictEta/D");
 517   settingsTree->Branch("fEtaOrdering", &fEtaOrdering,"EtaOrdering/O");
 518   settingsTree->Branch("fCutConversions", &fCutConversions,"CutConversions/O");
 519   settingsTree->Branch("fCutResonances", &fCutResonances,"CutResonances/O");
 520   settingsTree->Branch("fRejectResonanceDaughters", &fRejectResonanceDaughters,"RejectResonanceDaughters/I");
 521   settingsTree->Branch("fFillpT", &fFillpT,"FillpT/O");
 522   settingsTree->Branch("fMinNumTrack", &fMinNumTrack,"fMinNumTrack/I");
 523   settingsTree->Branch("fWeightPerEvent", &fWeightPerEvent,"WeightPerEvent/O");
 524   settingsTree->Branch("fPtOrder", &fPtOrder,"PtOrder/O");
 525   settingsTree->Branch("fTwoTrackEfficiencyCut", &fTwoTrackEfficiencyCut,"TwoTrackEfficiencyCut/D");
 526   settingsTree->Branch("twoTrackEfficiencyCutValue", &twoTrackEfficiencyCutValue,"twoTrackEfficiencyCutValue/D");
 527   settingsTree -> Fill();
 528   fList -> Add(settingsTree);
 529 }
 530
 531 //_______________________________________________________________________________
 532 void AliAnalysisTaskPidPidCorrelations::Analyse()
 533 {
 534   fHistEventStat -> Fill("Total",1); // all events
 535
 536   // store offline trigger bits
 537   fHistTriggerStats -> Fill(((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected());
 538
 539 //_________Trigger selection
 540   Bool_t isSelected = kTRUE; //Bool_t lTrigger = kTRUE;
 541   // instead of task->SelectCollisionCandidate(mask) in AddTask macro
 542   isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & fTriggerType);
 543   /*
 544   Printf("Trigger classes: %s:", fAOD -> GetFiredTriggerClasses().Data());
 545   Bool_t lTrigger1 = (fAOD -> GetFiredTriggerClasses().Contains("CVLN_R1-B-NOPF-ALLNOTRD")) ? 1 : 0;
 546   Bool_t lTrigger2 = (fAOD -> GetFiredTriggerClasses().Contains("CINT7WU-B-NOPF-ALL")) ? 1 : 0;
 547   lTrigger = lTrigger1 && lTrigger2;
 548   */
 549   if(!isSelected/* && !lTrigger*/) { /*AliInfo(" === Event Rejected!");*/ fHistEventStat -> Fill("No trigger",1); PostData(1,fList); return; }
 550
 551   TString trgClasses = fMyAODEvent -> GetFiredTriggerClasses();
 552 //   Printf("GetFiredTriggerClasses: %s",trgClasses.Data());
 553
 554   if (fMyAODEvent->GetFiredTriggerClasses().Contains("ALL")) fHistTriggerRun -> Fill(0);
 555   if (fMyAODEvent->GetFiredTriggerClasses().Contains("ALLNOTRD")) fHistTriggerRun -> Fill(1);
 556   if (fMyAODEvent->GetFiredTriggerClasses().Contains("MUON")) fHistTriggerRun -> Fill(2);
 557   if (fMyAODEvent->GetFiredTriggerClasses().Contains("CENT")) fHistTriggerRun -> Fill(3);
 558   if (fMyAODEvent->GetFiredTriggerClasses().Contains("CENTNOTRD")) fHistTriggerRun -> Fill(4);
 559   if (fMyAODEvent->GetFiredTriggerClasses().Contains("TPC")) fHistTriggerRun -> Fill(5);
 560
 561
 562   //_________Centrality
 563   if (fCentralityEstimator.Length() > 0)
 564   {
 565     AliCentrality* centralityObj = 0x0;
 566
 567     Int_t lCentralityClass10    = 0;
 568     Int_t lCentralityClass5     = 0;
 569
 570     fMyAODHeader = dynamic_cast<AliAODHeader*>(fMyAODEvent->GetHeader());
 571
 572     // QA for centrality estimators
 573     fHistCentStats->Fill(0.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("V0M"));
 574     fHistCentStats->Fill(1.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("FMD"));
 575     fHistCentStats->Fill(2.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("TRK"));
 576     fHistCentStats->Fill(3.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("TKL"));
 577     fHistCentStats->Fill(4.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("CL0"));
 578     fHistCentStats->Fill(5.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("CL1"));
 579     fHistCentStats->Fill(6.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("V0MvsFMD"));
 580     fHistCentStats->Fill(7.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("TKLvsV0M"));
 581     fHistCentStats->Fill(8.,fMyAODHeader->GetCentralityP()->GetCentralityPercentile("ZEMvsZDC"));
 582
 583     // centrality QA (V0M)
 584     fHistV0M->Fill(fMyAODEvent->GetVZEROData()->GetMTotV0A(), fMyAODEvent->GetVZEROData()->GetMTotV0C());
 585
 586     // centrality QA (reference tracks)
 587     fHistRefTracks->Fill(0.,fMyAODHeader->GetRefMultiplicity());
 588     fHistRefTracks->Fill(1.,fMyAODHeader->GetRefMultiplicityPos());
 589     fHistRefTracks->Fill(2.,fMyAODHeader->GetRefMultiplicityNeg());
 590     fHistRefTracks->Fill(3.,fMyAODHeader->GetTPConlyRefMultiplicity());
 591     fHistRefTracks->Fill(4.,fMyAODHeader->GetNumberOfITSClusters(0));
 592     fHistRefTracks->Fill(5.,fMyAODHeader->GetNumberOfITSClusters(1));
 593
 594 //     if (fMyAODEvent)
 595 //       centralityObj = fMyAODEvent -> GetCentrality();
 596     centralityObj = fMyAODHeader -> GetCentralityP();
 597
 598     Bool_t lCentralityInClass = kFALSE;
 599     if (centralityObj)
 600     {
 601 //       Printf("======= GetQuality = %f\n",centralityObj->GetQuality());
 602 //       fCentrality = centralityObj -> GetCentralityPercentileUnchecked(fCentralityEstimator.Data());
 603       fCentrality = centralityObj -> GetCentralityPercentile(fCentralityEstimator.Data());      // returns the centrality percentile, a float Form 0 to 100 (or to the trigger efficiency)
 604       lCentralityClass10        = centralityObj -> GetCentralityClass10(fCentralityEstimator.Data());//"FMD"); // returns centrality class for 10% width (a integer Form 0 to 10)
 605       lCentralityClass5 = centralityObj -> GetCentralityClass5(fCentralityEstimator.Data());//"TKL"); // returns centrality class for 5% width (a integer Form 0 to 20)
 606       lCentralityInClass        = centralityObj -> IsEventInCentralityClass(fCentralityPercentileMin,fCentralityPercentileMax,fCentralityEstimator.Data()); // returns kTRUE if the centrality of the event is between a and b, otherwise kFALSE
 607       if ( !lCentralityInClass ) { fHistEventStat -> Fill("Wrong centrality", 1); /*Printf("fCentrality = %f\n",fCentrality);*/ PostData(1, fList); return; }
 608
 609       //AliInfo(Form("Centrality is %f",fCentrality));
 610       fHistCentralityPercentile -> Fill( fCentrality );
 611       fHistCentralityClass10 -> Fill( lCentralityClass10 );
 612       fHistCentralityClass5 -> Fill( lCentralityClass5 );
 613     }
 614     else
 615     {
 616       Printf("WARNING: Centrality object is 0");
 617       fCentrality = -1.;
 618     }
 619   }
 620
 621   if (fCentrality < 0) return;
 622   Int_t centBin = GetCentBin(fCentrality);
 623   if (centBin<0) return;
 624
 625   //___ QA
 626   fHistRefTracksCent[centBin][0]->Fill(fCentrality,fMyAODHeader->GetRefMultiplicity());
 627   fHistRefTracksCent[centBin][1]->Fill(fCentrality,fMyAODHeader->GetRefMultiplicityPos());
 628   fHistRefTracksCent[centBin][2]->Fill(fCentrality,fMyAODHeader->GetRefMultiplicityNeg());
 629   fHistRefTracksCent[centBin][3]->Fill(fCentrality,fMyAODHeader->GetTPConlyRefMultiplicity());
 630   fHistRefTracksCent[centBin][4]->Fill(fCentrality,fMyAODHeader->GetNumberOfITSClusters(0));
 631   fHistRefTracksCent[centBin][5]->Fill(fCentrality,fMyAODHeader->GetNumberOfITSClusters(1));
 632
 633   //_________Vertex selection
 634   if(!VertexSelection(fMyAODEvent, fnTracksVertex, centBin, fVxMax, fVyMax, fVzMax)) return;
 635
 636   //_________check the PIDResponse handler
 637   if (!fPIDResponse) return;
 638
 639   fHistEventStat -> Fill("Analyzed",1);
 640
 641
 642   if (!fUseMC)
 643   {
 644     //_________get magnetic field
 645     Double_t bSign              = (fMyAODEvent -> GetMagneticField() > 0) ? 1 : -1;  // for two track efficiency cut
 646     //Double_t bSignAux = fMyAODHeader -> GetMagneticField() ; // for dca cut
 647
 648     //____ Get AOD RECO tracks
 649     fMyprimRecoTracksPID = AcceptTracks(centBin,0x0,/*kFALSE,*/kTRUE); // !FIXME onlyprimaries part
 650     //Printf("Accepted %d fMyprimRecoTracksPID", tracks->GetEntries());
 651     if(fMyprimRecoTracksPID==0x0) { AliInfo(" ==== fMyprimRecoTracksPID: Zero track pointer"); return; }
 652     //AliDebug(1, Form("Single Event analysis : nTracks = %4d", fMyprimRecoTracksPID -> GetEntries()));
 653
 654     if (fRejectZeroTrackEvents && fMyprimRecoTracksPID->GetEntriesFast() == 0)
 655     {
 656       AliInfo(Form("========== Rejecting event because it has no tracks: %f %d", fCentrality, fMyprimRecoTracksPID->GetEntriesFast()));
 657       fHistEventStat -> Fill("NO tracks per event",1);
 658       //delete fMyprimRecoTracksPID; fMyprimRecoTracksPID=0x0;
 659       return;
 660     }
 661     //delete fMyprimRecoTracksPID;
 662
 663     const AliVVertex* vertex = fMyAODEvent->GetPrimaryVertex();
 664     Double_t zVtx = vertex->GetZ();
 665
 666     //_________ Fill Correction AliCFContainer
 667     FillCFcontainers(0x0, fMyprimRecoTracksPID, 0x0, fCentrality/*, zVtx*/);
 668
 669     //_________ Fill Correlations
 670     Double_t weight = 1.;
 671     if (fFillpT)
 672       weight = -1.;
 673
 674     //____ step 0
 675     //FillCorrelations(fMyprimRecoTracksPID, 0x0, fCentrality, zVtx, 0.0, kFALSE, 0.02, kFALSE, /*0,*/ weight);
 676
 677     //____ step 1
 678     if (fTwoTrackEfficiencyCut > 0)
 679       FillCorrelations(fMyprimRecoTracksPID, 0x0, fCentrality, zVtx, bSign, kTRUE, twoTrackEfficiencyCutValue, /*1,*/ weight);
 680
 681     if (fFillMixed)
 682     {
 683       // event mixing
 684
 685       // 1. First get an event pool corresponding in mult (cent) and
 686       //    zvertex to the current event. Once initialized, the pool
 687       //    should contain nMix (reduced) events. This routine does not
 688       //    pre-scan the chain. The first several events of every chain
 689       //    will be skipped until the needed pools are filled to the
 690       //    specified depth. If the pool categories are not too rare, this
 691       //    should not be a problem. If they are rare, you could lose
 692       //    statistics.
 693
 694       // 2. Collect the whole pool's content of tracks into one TObjArray
 695       //    (bgTracks), which is effectively a single background super-event.
 696
 697       // 3. The reduced and bgTracks arrays must both be passed into
 698       //    FillCorrelations(). Also nMix should be passed in, so a weight
 699       //    of 1./nMix can be applied.
 700
 701       AliEventPool* pool = fPoolMgr->GetEventPool(fCentrality, zVtx);
 702
 703       if (!pool)
 704         AliInfo(Form("No pool found for centrality = %f, zVtx = %f", fCentrality, zVtx));
 705       else
 706       {
 707 //       pool->SetDebug(1);
 708 //       pool->PrintInfo();
 709 //       PrintPoolManagerContents();
 710
 711         //if (pool->IsReady())
 712         if (pool->IsReady() || pool->NTracksInPool() > fMinNumTrack / 10 || pool->GetCurrentNEvents() >= fMinNEventsToMix)
 713         {
 714           Int_t nMix = pool->GetCurrentNEvents();
 715 //       cout << "nMix = " << nMix << " tracks in pool = " << pool->NTracksInPool() << endl;
 716
 717           // ! FIXME later
 718 //       ((TH1F*) fList->FindObject("eventStat"))->Fill(2);
 719 //       ((TH2F*) fList->FindObject("mixedDist"))->Fill(centrality, pool->NTracksInPool());
 720 //       if (pool->IsReady())
 721 //      ((TH1F*) fList->FindObject("eventStat"))->Fill(3);
 722
 723           // Fill mixed-event histos here
 724           for (Int_t jMix=0; jMix<nMix; jMix++)
 725           {
 726             TObjArray* bgTracks = pool->GetEvent(jMix);
 727             if (!bgTracks) continue;
 728
 729             if (fTwoTrackEfficiencyCut > 0)
 730               FillCorrelations(fMyprimRecoTracksPID, bgTracks, fCentrality, zVtx, bSign, kTRUE, twoTrackEfficiencyCutValue, /*2,*/ 1./nMix);
 731           }
 732         }
 733         pool->UpdatePool(fMyprimRecoTracksPID);
 734         //pool->PrintInfo();
 735       }//_____ pool NULL check
 736     }//_______ run mixing
 737     else
 738       delete fMyprimRecoTracksPID;
 739   }
 740   else
 741   {
 742     //______ check MC
 743     Bool_t isMCok = kFALSE;
 744     isMCok = CheckMcDistributions(fMyMcArray,fMyMcHeader);
 745     if (isMCok == kFALSE) return;
 746
 747     //____ Get MC primaries
 748     fMyprimMCParticlesPID = AcceptMcTracks(centBin,kTRUE,kTRUE);
 749     if (fMyprimMCParticlesPID==0x0) { AliInfo(" ==== fMyprimMCParticlesPID: Zero track pointer"); return; }
 750     //CleanUp(primMCParticles, fMyMcArray, skipParticlesAbove);
 751 //     delete fMyprimMCParticlesPID; fMyprimMCParticlesPID=0x0;
 752     //____ Get MC RECO-matched
 753     fMyprimRecoTracksMatchedPID = AcceptTracks(centBin,fMyMcArray,/*kTRUE,*/kTRUE);
 754     if (fMyprimRecoTracksMatchedPID==0x0) { AliInfo(" ==== Zero track pointer"); return; }
 755     //TObjArray* fMyprimRecoTracksMatchedPID    = AcceptMcRecoMachedTracks(centBin,kTRUE,kTRUE);
 756     //CleanUp(fMyprimRecoTracksMatchedPID, fMyMcArray, skipParticlesAbove);
 757 //     delete fMyprimRecoTracksMatchedPID; fMyprimRecoTracksMatchedPID=0x0;
 758
 759     //___ check pdg
 760     for (Int_t i=0; i<fMyMcArray->GetEntriesFast(); i++)
 761       ((TH1F*) fList->FindObject("pids"))->Fill(((AliAODMCParticle*)fMyMcArray->At(i))->GetPdgCode());
 762
 763     //______ Fill Corrections
 764     FillCFcontainers(fMyprimMCParticlesPID, 0x0, fMyprimRecoTracksMatchedPID, fCentrality/*, zVtx*/);
 765 //     Printf("%d %d %d" fMyprimMCParticlesPID->GetEntries(), fMyprimRecoTracksPID->GetEntries(), fMyprimRecoTracksMatchedPID->GetEntries());
 766 //     delete fMyprimMCParticlesPID;
 767 //     delete fMyprimRecoTracksMatchedPID;
 768   }
 769   return;
 770 }
 771
 772
 773 //________________________________________________________________________
 774 void AliAnalysisTaskPidPidCorrelations::FillCorrelations(TObjArray* particles, TObjArray* particlesMixed, Double_t centrality, Double_t zVtx, Double_t bSign, Bool_t twoTrackEfficiencyCut, Double_t fTwoTrackEfficiencyCutValue, /*Int_t step,*/ Double_t weight)
 775 {
 776   if (!particles){ AliInfo(" ==============  particles TObjArray is NULL pointer, return"); return; }
 777
 778 //   Double_t* trackVariablesSingle = new Double_t[kTrackVariablesSingle];
 779   Double_t* trackVariablesPair = new Double_t[kTrackVariablesPair];
 780
 781   Bool_t fillpT = kFALSE;
 782   if (weight < 0) fillpT = kTRUE;
 783
 784   // define end of particle loops
 785   Int_t iMax = particles -> GetEntriesFast();
 786   Int_t jMax = iMax;
 787   if (particlesMixed)
 788     jMax = particlesMixed -> GetEntriesFast();
 789
 790   // Eta() is extremely time consuming, therefore cache it for the inner loop here:
 791   TObjArray* input = (particlesMixed) ? particlesMixed : particles;
 792
 793   TArrayF secondEta(jMax);
 794   TArrayF secondPhi(jMax);
 795   TArrayF secondPt(jMax);
 796 //   TArrayS secondCharge(jMax);
 797   TArrayF secondPID(jMax);
 798
 799   for (Int_t i=0;i < jMax;i++)
 800   {
 801     secondEta[i]                = ((AliPidPidCorrelationReducedTrack*)input->At(i))->Eta();
 802     secondPhi[i]                = ((AliPidPidCorrelationReducedTrack*)input->At(i))->Phi();
 803     secondPt[i]         = ((AliPidPidCorrelationReducedTrack*)input->At(i))->Pt();
 804 //     secondCharge[i]  = (Short_t)((AliPidPidCorrelationReducedTrack*)input->At(i))->Charge();
 805     secondPID[i]                = ((AliPidPidCorrelationReducedTrack*)input->At(i))->GetMyPartID();
 806   }
 807
 808
 809   //_______ Reject ResonanceDaughters
 810
 811   // identify K, Lambda candidates and flag those particles
 812   // a TObject bit is used for this
 813   const UInt_t kResonanceDaughterFlag = 1 << 14;
 814   if (fRejectResonanceDaughters > 0)
 815   {
 816       Double_t resonanceMass = -1;
 817       Double_t massDaughter1 = -1;
 818       Double_t massDaughter2 = -1;
 819       const Double_t interval = 0.02;
 820
 821       switch (fRejectResonanceDaughters)
 822       {
 823         case 1: resonanceMass = 1.2;    massDaughter1 = 0.1396; massDaughter2 = 0.9383;                 break; // method test
 824         case 2: resonanceMass = 0.4976; massDaughter1 = 0.1396; massDaughter2 = massDaughter1;  break; // k0
 825         case 3: resonanceMass = 1.115;  massDaughter1 = 0.1396; massDaughter2 = 0.9383;                 break; // lambda
 826         default: AliFatal(Form("Invalid setting %d", fRejectResonanceDaughters));
 827       }
 828
 829       for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 830         particles->At(i)->ResetBit(kResonanceDaughterFlag);
 831
 832       if (particlesMixed)
 833         for (Int_t i=0; i<jMax; i++)
 834           particlesMixed->At(i)->ResetBit(kResonanceDaughterFlag);
 835
 836       for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 837       {
 838         //AliVParticle* triggerParticle = (AliVParticle*) particles->UncheckedAt(i);
 839         AliPidPidCorrelationReducedTrack* triggerParticle = (AliPidPidCorrelationReducedTrack*) particles->At(i);
 840
 841         for (Int_t j=0; j<jMax; j++)
 842         {
 843           if (!particlesMixed && i == j)        continue;
 844
 845           AliVParticle* particle = 0;
 846           if (!particlesMixed)
 847           {
 848             particle = (AliPidPidCorrelationReducedTrack*) particles->At(j);
 849             //particle = (AliVParticle*) particles->UncheckedAt(j);
 850           }
 851           else
 852           {
 853             particle = (AliPidPidCorrelationReducedTrack*) particlesMixed->At(j);
 854             //particle = (AliVParticle*) particlesMixed->UncheckedAt(j);
 855           }
 856
 857           // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
 858 //        if (particlesMixed && triggerParticle->IsEqual(particle)) continue;
 859
 860           if (triggerParticle->Charge() * particle->Charge() > 0) continue;
 861
 862           Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerParticle->Eta(), triggerParticle->Phi(), particle->Pt(), particle->Eta(), particle->Phi(), massDaughter1, massDaughter2);
 863
 864           if (TMath::Abs(mass - resonanceMass*resonanceMass) < interval*5)
 865           {
 866             mass = GetInvMassSquared(triggerParticle->Pt(), triggerParticle->Eta(), triggerParticle->Phi(), particle->Pt(), particle->Eta(), particle->Phi(), massDaughter1, massDaughter2);
 867
 868             if (mass > (resonanceMass-interval)*(resonanceMass-interval) && mass < (resonanceMass+interval)*(resonanceMass+interval))
 869             {
 870               triggerParticle->SetBit(kResonanceDaughterFlag);
 871               particle->SetBit(kResonanceDaughterFlag);
 872               //Printf("Flagged %d %d %f", i, j, TMath::Sqrt(mass));
 873             }
 874           }
 875         }
 876       }
 877     }//________ fRejectResonanceDaughters
 878
 879   //_______ fHistTriggerWeighting
 880   if (fWeightPerEvent)
 881   {
 882     for (Int_t i=0; i<particles->GetEntriesFast(); i++)
 883     {
 884       AliVParticle* triggerParticle = (AliVParticle*) particles->At(i);
 885
 886       // some optimization
 887       Double_t triggerEta = triggerParticle->Eta();
 888
 889       if (fTriggerRestrictEta > 0 && TMath::Abs(triggerEta) > fTriggerRestrictEta)
 890         continue;
 891
 892       if (fOnlyOneEtaSide != 0)
 893         if (fOnlyOneEtaSide * triggerEta < 0)
 894           continue;
 895       if (fTriggerSelectCharge != 0)
 896         if (triggerParticle->Charge() * fTriggerSelectCharge < 0)
 897           continue;
 898
 899       fHistTriggerWeighting->Fill(triggerParticle->Pt());
 900      }
 901   }
 902
 903   //_______ trigger particle loop
 904   for (Int_t i=0; i < iMax; i++)
 905   {
 906 //     AliVParticle* triggerParticle = (AliVParticle*) particles -> UncheckedAt(i);
 907     AliPidPidCorrelationReducedTrack* triggerParticle = dynamic_cast<AliPidPidCorrelationReducedTrack*>(particles -> At(i));
 908     if (!triggerParticle) { AliInfo(" ===== triggerParticle: zero pointer. Continue...\n"); continue; }
 909
 910     Double_t triggerEta = triggerParticle -> Eta();
 911     Double_t triggerPhi = triggerParticle -> Phi();
 912     Double_t triggerPt  = triggerParticle -> Pt();
 913     Int_t triggerPID    = triggerParticle -> GetMyPartID();
 914
 915     if (fTriggerRestrictEta > 0 && TMath::Abs(triggerEta) > fTriggerRestrictEta)
 916       continue;
 917
 918     if (fOnlyOneEtaSide != 0)
 919       if (fOnlyOneEtaSide * triggerEta < 0)
 920         continue;
 921
 922     if (fTriggerSelectCharge != 0)
 923       if (triggerParticle->Charge() * fTriggerSelectCharge < 0)
 924         continue;
 925
 926     if (fRejectResonanceDaughters > 0)
 927       if (triggerParticle->TestBit(kResonanceDaughterFlag))
 928       {
 929         //Printf("Skipped i=%d", i);
 930         continue;
 931       }
 932
 933       //Int_t CentBin = fHistCorrPairSame->GetGrid(0)->GetGrid()->GetAxis(7)->FindBin(centrality);//GetCentBin(centrality);
 934       //Int_t ZvtxBin = fHistCorrPairSame->GetGrid(0)->GetGrid()->GetAxis(6)->FindBin(zVtx);//GetZvtxBin(zVtx);
 935
 936     //fEtaSpectraTrig[CentBin][ZvtxBin] -> Fill(triggerEta);
 937
 938     // ___ trigger eta PID
 939     //if ( triggerPID == fPartPionPlus || triggerPID == fPartPionMinus )                fEtaSpectraTrigPion[CentBin][ZvtxBin] -> Fill(triggerEta);
 940     //if ( triggerPID == fPartKaonPlus || triggerPID == fPartKaonMinus )                fEtaSpectraTrigKaon[CentBin][ZvtxBin] -> Fill(triggerEta);
 941     //if ( triggerPID == fPartProtonPlus || triggerPID == fPartProtonMinus )    fEtaSpectraTrigProton[CentBin][ZvtxBin] -> Fill(triggerEta);
 942
 943     //___ fill single particle histograms
 944 /*
 945     trackVariablesSingle[0]    = triggerPID;
 946     trackVariablesSingle[1]    = triggerPt;
 947     trackVariablesSingle[2]    = zVtx;
 948     trackVariablesSingle[3]    = centrality;
 949     trackVariablesSingle[4]    = triggerParticle->Charge()>0 ? 7 : 8;
 950     */
 951 //     if( fWriteCorrTree == kTRUE && !particlesMixed)
 952 //     {
 953 //       fCorrVariables[0] = triggerPID;
 954 //       fCorrVariables[1] = triggerPt;
 955 //       fCorrVariables[2] = zVtx;
 956 //       fCorrVariables[3] = centrality;
 957 //       fCorrVariables[4] = triggerParticle->Charge()>0 ? 7 : 8;
 958 //
 959 //       fVariablesTreeCorr->Fill();
 960 //     }
 961
 962     if (fillpT) weight = triggerPt;
 963     Double_t useWeight = weight;
 964
 965     if (!particlesMixed)
 966     {
 967 //       fHistCorrSingle -> Fill(trackVariablesSingle,step,useWeight);
 968 //       fHistCorrSingle -> Fill(trackVariablesSingle,useWeight);
 969
 970       fHistSingleHadronsPt[GetCentBin(centrality)]->Fill(triggerPt);
 971       fHistSingleHadronsEtaPt[GetCentBin(centrality)]->Fill(triggerPt,triggerEta);
 972     }
 973
 974     //_______ assoc particle loop
 975     for (Int_t j=0; j < jMax; j++)
 976     {
 977       // no auto correlations (only for non mixing)
 978       if (!particlesMixed && i == j) continue;
 979
 980       //AliVParticle* assocParticle = 0x0;
 981       AliPidPidCorrelationReducedTrack* assocParticle = 0x0;
 982       if (!particlesMixed)
 983       {
 984         //assocParticle = (AliVParticle*) particles->At(j);
 985         assocParticle = (AliPidPidCorrelationReducedTrack*) particles->At(j);
 986       }
 987       else
 988       {
 989         //assocParticle = (AliVParticle*) mixed->At(j);
 990         assocParticle = (AliPidPidCorrelationReducedTrack*) particlesMixed->At(j);
 991       }
 992
 993       if (!assocParticle) continue;
 994
 995       // check if both particles point to the same element (does not occur for mixed events, but if subsets are mixed within the same event)
 996       //! FIXME
 997       //       if (particlesMixed && triggerParticle -> IsEqual(assocParticle))
 998       //        continue;
 999
1000       if (fPtOrder)
1001         if (assocParticle->Pt() >= triggerParticle->Pt())
1002           continue;
1003
1004       if (fAssociatedSelectCharge != 0)
1005         if (assocParticle->Charge() * fAssociatedSelectCharge < 0)
1006           continue;
1007
1008       if (fSelectCharge > 0)
1009       {
1010         // skip like sign
1011         if (fSelectCharge == 1 && assocParticle->Charge() * triggerParticle->Charge() > 0)      continue;
1012         // skip unlike sign
1013         if (fSelectCharge == 2 && assocParticle->Charge() * triggerParticle->Charge() < 0)      continue;
1014       }
1015
1016       if (fEtaOrdering)
1017       {
1018         if (triggerEta < 0 && secondEta[j] < triggerEta)        continue;
1019         if (triggerEta > 0 && secondEta[j] > triggerEta)        continue;
1020       }
1021
1022       if (fRejectResonanceDaughters > 0)
1023         if (triggerParticle->TestBit(kResonanceDaughterFlag))
1024         {
1025           //Printf("Skipped i=%d", i);
1026           continue;
1027         }
1028
1029       // conversions
1030       if (fCutConversions && assocParticle->Charge() * triggerParticle->Charge() < 0)
1031       {
1032         Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.510e-3, 0.510e-3);
1033
1034         if (mass < 0.1)
1035         {
1036           mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.510e-3, 0.510e-3);
1037           fHistControlConvResoncances->Fill(0.0, mass);
1038
1039           if (mass < 0.04*0.04) continue;
1040         }
1041       }
1042
1043       // K0s
1044       if (fCutResonances && assocParticle->Charge() * triggerParticle->Charge() < 0)
1045       {
1046         Float_t mass = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.1396, 0.1396);
1047
1048         const Float_t kK0smass = 0.4976;
1049
1050         if (TMath::Abs(mass - kK0smass*kK0smass) < 0.1)
1051         {
1052           mass = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.1396, 0.1396);
1053           fHistControlConvResoncances->Fill(1, mass - kK0smass*kK0smass);
1054
1055           if (mass > (kK0smass-0.02)*(kK0smass-0.02) && mass < (kK0smass+0.02)*(kK0smass+0.02))
1056             continue;
1057         }
1058       }
1059
1060       // Lambda
1061       if (fCutResonances && assocParticle->Charge() * triggerParticle->Charge() < 0)
1062       {
1063         Float_t mass1 = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.1396, 0.9383);
1064         Float_t mass2 = GetInvMassSquaredCheap(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.9383, 0.1396);
1065
1066         const Float_t kLambdaMass = 1.115;
1067
1068         if (TMath::Abs(mass1 - kLambdaMass*kLambdaMass) < 0.1)
1069         {
1070           mass1 = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.1396, 0.9383);
1071           fHistControlConvResoncances->Fill(2, mass1 - kLambdaMass*kLambdaMass);
1072
1073           if (mass1 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass1 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
1074             continue;
1075         }
1076         if (TMath::Abs(mass2 - kLambdaMass*kLambdaMass) < 0.1)
1077         {
1078           mass2 = GetInvMassSquared(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), assocParticle->Pt(), secondEta[j], assocParticle->Phi(), 0.9383, 0.1396);
1079           fHistControlConvResoncances->Fill(2, mass2 - kLambdaMass*kLambdaMass);
1080
1081           if (mass2 > (kLambdaMass-0.02)*(kLambdaMass-0.02) && mass2 < (kLambdaMass+0.02)*(kLambdaMass+0.02))
1082             continue;
1083         }
1084       }
1085
1086       //_______ two-Track Efficiency Cut , HBT-like cut
1087       if (twoTrackEfficiencyCut)
1088       {
1089         // the variables & cuthave been developed by the HBT group
1090         // see e.g. https://indico.cern.ch/materialDisplay.py?contribId=36&sessionId=6&materialId=slides&confId=142700
1091
1092         Double_t phi1   = triggerParticle -> Phi();
1093         //Double_t phi1 = triggerParticle -> Phi()*TMath::DegToRad();
1094         Double_t pt1    = triggerParticle -> Pt();
1095         Double_t charge1        = triggerParticle -> Charge();
1096
1097         //Double_t phi2 = assocParticle -> Phi()*TMath::DegToRad();
1098         Double_t phi2   = assocParticle -> Phi();
1099         Double_t pt2    = assocParticle -> Pt();
1100         Double_t charge2        = assocParticle -> Charge();
1101
1102         Double_t deta = triggerEta - secondEta[j];
1103         Double_t dphi = triggerPhi - secondPhi[j];
1104
1105         fHistHBTbefore -> Fill(deta,dphi);
1106
1107         // optimization
1108         if (TMath::Abs(deta) < fTwoTrackEfficiencyCutValue * 2.5 * 3) // twoTrackEfficiencyCutValue = 0.02 [default for dphicorrelations]
1109         {
1110           // check first boundaries to see if is worth to loop and find the minimum
1111           Double_t dphistar1 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 0.8, bSign);
1112           Double_t dphistar2 = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, 2.5, bSign);
1113
1114           const Double_t kLimit = fTwoTrackEfficiencyCutValue * 3;
1115
1116           Double_t dphistarminabs = 1e5;
1117           Double_t dphistarmin = 1e5;
1118           if (TMath::Abs(dphistar1) < kLimit || TMath::Abs(dphistar2) < kLimit || dphistar1 * dphistar2 < 0)
1119           {
1120             for (Double_t rad=0.8; rad<2.51; rad+=0.01)
1121             {
1122               Double_t dphistar = GetDPhiStar(phi1, pt1, charge1, phi2, pt2, charge2, rad, bSign);
1123               Double_t dphistarabs = TMath::Abs(dphistar);
1124
1125               if (dphistarabs < dphistarminabs)
1126               {
1127                 dphistarmin = dphistar;
1128                 dphistarminabs = dphistarabs;
1129               }
1130             }
1131
1132             Int_t PtBin = GetPtBin(pt2);
1133 //          Printf("================ pt = %f \t <--> \t PtBin = %d\n",pt2,PtBin);
1134             //fHistTwoTrackDistancePt[0] -> Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2)); // !FIXME
1135             fHistTwoTrackDistancePt[PtBin][0] -> Fill(deta, dphistarmin);
1136
1137             if (dphistarminabs < fTwoTrackEfficiencyCutValue && TMath::Abs(deta) < fTwoTrackEfficiencyCutValue)
1138             {
1139               // Printf("Removed track pair %d %d with %f %f %f %f %f %f %f %f %f", i, j, deta, dphistarminabs, phi1, pt1, charge1, phi2, pt2, charge2, bSign);
1140               // AliInfo(Form("HBT: Removed track pair %d %d with [[%f %f]] %f %f %f | %f %f %d %f %f %d %f", i, j, deta, dphi, dphistarminabs, dphistar1, dphistar2, phi1rad, pt1, charge1, phi2rad, pt2, charge2, bSign));
1141               continue;
1142             }
1143
1144             //fHistTwoTrackDistancePt[1] -> Fill(deta, dphistarmin, TMath::Abs(pt1 - pt2)); // !FIXME
1145             fHistTwoTrackDistancePt[PtBin][1] -> Fill(deta, dphistarmin);
1146           }
1147         }
1148         fHistHBTafter -> Fill(deta,dphi);
1149       }//____ HBT cut
1150
1151       trackVariablesPair[0] =  triggerPt;                               // pt trigger
1152       trackVariablesPair[1] =  secondPt[j];                     // pt assoc
1153       trackVariablesPair[2] =  DeltaPhi(triggerPhi-secondPhi[j]);       // delta phi
1154       trackVariablesPair[3] =  triggerEta - secondEta[j];               // delta eta
1155       trackVariablesPair[4] =  zVtx;                            // z of the primary vertex
1156
1157 //       trackVariablesPair[0] =  triggerPID;                   // trigger PID
1158 // //       trackVariablesPair[1] =  secondPID[j];                      // assoc PID
1159 //       trackVariablesPair[1] =  triggerPt;                            // pt trigger
1160 //       trackVariablesPair[2] =  secondPt[j];                  // pt assoc
1161 //       trackVariablesPair[3] =  DeltaPhi(triggerPhi-secondPhi[j]);    // delta phi
1162 //       trackVariablesPair[4] =  triggerEta - secondEta[j];            // delta eta
1163 //       trackVariablesPair[5] =  zVtx;                         // z of the primary vertex
1164 // //       trackVariablesPair[7] =  triggerParticle->Charge()>0 ? 7 : 8;
1165 //       trackVariablesPair[6] =  assocParticle->Charge()>0 ? 7 : 8;
1166
1167       //_______ momentum difference cut - suppress femtoscopic effects
1168       if (fQCut)
1169       {
1170         //Double_t ptMin        = 0.1; //const for the time being (should be changeable later on)
1171         Double_t ptDifference = TMath::Abs( triggerPt - secondPt[j]);
1172 //      fHistQbefore -> Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
1173         if (ptDifference < fDeltaPtMin) continue;
1174 //      fHistQafter -> Fill(trackVariablesPair[1],trackVariablesPair[2],ptDifference);
1175       }
1176
1177       if (fWeightPerEvent)
1178       {
1179         Int_t weightBin = fHistTriggerWeighting->GetXaxis()->FindBin(trackVariablesPair[0]);
1180         //Printf("Using weight %f", fHistTriggerWeighting->GetBinContent(weightBin));
1181         useWeight /= fHistTriggerWeighting->GetBinContent(weightBin);
1182       }
1183
1184       if (!particlesMixed)
1185       {
1186         if (triggerPID == fPIDtrigType && secondPID[j] == fPIDassocType) fHistCorrPair[0] -> Fill(trackVariablesPair,useWeight);
1187 //      fHistCorrPairSame -> Fill(trackVariablesPair,step,useWeight);
1188 /*
1189         Int_t ptBinTrackTrig  = fPtAxis -> FindBin( triggerParticle->Pt() );
1190         Int_t ptBinTrackAssoc = fPtAxis -> FindBin( assocParticle->Pt() );
1191
1192         if (ptBinTrackTrig<1 || ptBinTrackTrig>fNbinsPt || ptBinTrackAssoc<1 || ptBinTrackAssoc>fNbinsPt) continue;
1193
1194         fHistDeltaPhi[GetCentBin(centrality)][ptBinTrackTrig-1][ptBinTrackAssoc-1]      -> Fill(DeltaPhi(triggerParticle->Phi()-assocParticle->Phi()));
1195         fHistDphiDeta[GetCentBin(centrality)][ptBinTrackTrig-1][ptBinTrackAssoc-1]      -> Fill(DeltaPhi(triggerParticle->Phi()-assocParticle->Phi()),triggerParticle->Eta()-assocParticle->Eta());
1196         fHistTracksEtaTrigVsEtaAssoc[GetCentBin(centrality)]                    -> Fill(triggerParticle -> Eta(),assocParticle -> Eta());*/
1197       }
1198       else if (particlesMixed)
1199       {
1200         if (triggerPID == fPIDtrigType && secondPID[j] == fPIDassocType) fHistCorrPair[1] -> Fill(trackVariablesPair,useWeight);
1201 //      fHistCorrPairMixed -> Fill(trackVariablesPair,step,useWeight);
1202
1203 //      Int_t ptBinTrackTrigM  = fPtAxis -> FindBin( triggerParticle->Pt() );
1204 //      Int_t ptBinTrackAssocM = fPtAxis -> FindBin( assocParticle->Pt() );
1205 //
1206 //      if (ptBinTrackTrigM<1 || ptBinTrackTrigM>fNbinsPt || ptBinTrackAssocM<1 || ptBinTrackAssocM>fNbinsPt) continue;
1207 //
1208 //      fHistDeltaPhiMix[GetCentBin(centrality)][ptBinTrackTrigM-1][ptBinTrackAssocM-1] -> Fill(DeltaPhi(triggerParticle->Phi()-assocParticle->Phi()));
1209 //      fHistDphiDetaMix[GetCentBin(centrality)][ptBinTrackTrigM-1][ptBinTrackAssocM-1] -> Fill(DeltaPhi(triggerParticle->Phi()-assocParticle->Phi()),triggerParticle->Eta()-assocParticle->Eta());
1210 //      fHistTracksEtaTrigVsEtaAssocMixed[GetCentBin(centrality)]                               -> Fill(triggerParticle->Eta(),assocParticle->Eta());
1211       }
1212     }//_____ assoc loop
1213   }//_____ trigger loop
1214 //   delete[] trackVariablesSingle;
1215   delete[] trackVariablesPair;
1216 }
1217
1218 //_____________________________________________________
1219 Bool_t AliAnalysisTaskPidPidCorrelations::CheckMcDistributions(TClonesArray* arrayMC, AliAODMCHeader* mcHeader)
1220 {
1221   Int_t myPDGproton = 2212;
1222   Int_t myPDGkaon = 321;
1223   Int_t myPDGpion = 211;
1224 //   Int_t tmpPDG = -1;
1225
1226   //____ counters
1227   Int_t kCntr =0;
1228   Int_t kCntrInjected = 0;
1229   Int_t kCntrFromDecay = 0;
1230
1231   Int_t entrMC = arrayMC->GetEntriesFast();
1232   if ( entrMC < 1 ) return kFALSE;
1233     //
1234   AliAODMCParticle* mcPart = 0x0;
1235
1236   TString kGenName="";
1237 //     Bool_t isFromHijing = kFALSE;
1238
1239   for(Int_t ie=0 ;ie < entrMC;ie++)
1240   {
1241 //     isFromHijing = kFALSE;
1242
1243     mcPart = dynamic_cast<AliAODMCParticle*>(arrayMC->At(ie));
1244     if(!ie) continue;
1245
1246     if ( (TMath::Abs(mcPart->GetPdgCode()) != myPDGproton) && (TMath::Abs(mcPart->GetPdgCode()) != myPDGkaon) && (TMath::Abs(mcPart->GetPdgCode()) != myPDGpion) ) continue;
1247
1248     if( mcPart->IsPrimary() == kFALSE ) continue;
1249
1250     kCntr++;
1251
1252     ((TH1F*)(fList->FindObject("fHistMcStats")))->Fill(0);
1253     //
1254     kGenName = GetGenerator(ie,mcHeader);
1255     //AliInfo(Form("====> IsPrimary: %d IsPhysicsPrimary %d Gen %s",mcPart->IsPrimary(),mcPart->IsPhysicalPrimary(),kGenName.Data()));
1256
1257     FillMcGeneratorHistos(kGenName);
1258
1259     if ( kGenName.Contains("Hijing")) fMyMcType = 0;
1260     else if ( kGenName.Contains("Pythia")) fMyMcType = 1;
1261     else fMyMcType = -1;
1262
1263     Int_t idxmother = mcPart->GetMother();
1264     if (idxmother == -1) kCntrInjected++;
1265     else kCntrFromDecay++;
1266
1267     ((TH1F*)(fList->FindObject("fHistMcAllPt")))->Fill(mcPart->Pt());
1268     if (idxmother == -1 )       ((TH1F*)(fList->FindObject("fHistMcAllPt_Inj")))->Fill(mcPart->Pt());
1269     if (idxmother > -1 )        ((TH1F*)(fList->FindObject("fHistMcAllPt_Dec")))->Fill(mcPart->Pt());
1270
1271     ((TH1F*)(fList->FindObject("fHistMcAllEta")))->Fill(mcPart->Eta());
1272
1273     if (kGenName.Contains("ijing")) //_____ from HJING
1274     {
1275 //     isFromHijing = kTRUE;
1276       ((TH1F*)(fList->FindObject("fHistMcAllPt_Hijing")))->Fill(mcPart->Pt());
1277       ((TH1F*)(fList->FindObject("fHistMcAllEta_Hijing")))->Fill(mcPart->Eta());
1278     }
1279     else
1280       ((TH1F*)(fList->FindObject("fHistMcAllEta_NotHijing")))->Fill(mcPart->Eta());
1281   } // mc particle loop
1282
1283   ((TH1F*)(fList->FindObject("fHistNumOfPartPerEvt")))->Fill(kCntr);
1284   //AliInfo(Form(" === MC: \t injected \t %d, \t Decay: \t %d\n",kCntrInjected,kCntrFromDecay));
1285
1286   return kTRUE;
1287 }
1288
1289 //________________________________________________________________________
1290 TString AliAnalysisTaskPidPidCorrelations::GetGenerator(Int_t label, AliAODMCHeader* MCheader)
1291 {
1292   // taken from .../PWGHF/vertexingHF/AliVertexingHFUtils.cxx
1293
1294   // get the name of the generator that produced a given particle
1295   Int_t nsumpart = 0;
1296   TList*lh = MCheader -> GetCocktailHeaders();
1297   Int_t nh = lh -> GetEntries();
1298   for(Int_t i=0;i < nh;i++)
1299   {
1300     AliGenEventHeader* gh = (AliGenEventHeader*)lh -> At(i);
1301     TString genname = gh -> GetName();
1302     Int_t npart = gh -> NProduced();
1303     if ( label >= nsumpart && label < (nsumpart+npart) )        return genname;
1304     nsumpart += npart;
1305   }
1306   TString empty="";
1307   return empty;
1308 }
1309
1310 //________________________________________________
1311 Bool_t AliAnalysisTaskPidPidCorrelations::IsFromHijingBg(Int_t mcTrackLabel)
1312 {
1313   Bool_t isFromHijingBg = kFALSE;
1314
1315   if( mcTrackLabel < 0 ) return isFromHijingBg;
1316
1317 //   AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(fMyMcArray->At(mcTrackLabel));
1318
1319   if ( !mcTrackLabel ) return isFromHijingBg;
1320   if ( mcTrackLabel < 0 ) return isFromHijingBg;
1321
1322   TString nameGen = GetGenerator(mcTrackLabel,fMyMcHeader);
1323
1324   if ( nameGen.Contains("ijing") ) isFromHijingBg = kTRUE;
1325   if ( isFromHijingBg == kTRUE )
1326   {
1327     AliInfo(" === isFromHijingBg = kTRUE === \n\n");
1328     ((TH1F*)(fList->FindObject("fHist_HijingBg"))) -> Fill(0);
1329   }
1330   else
1331     ((TH1F*)(fList->FindObject("fHist_HijingBg"))) -> Fill(1);
1332
1333   return isFromHijingBg;
1334 }
1335
1336 //________________________________________________________________________
1337 void AliAnalysisTaskPidPidCorrelations::FillMcGeneratorHistos(TString genLabel)
1338 {
1339   Int_t fillVal = 0;
1340
1341   if ( genLabel.Contains("hijing_0"))   fillVal = 0;
1342   else if ( genLabel.Contains("pythia_1")) fillVal = 1;
1343   else
1344   {
1345     for(Int_t i=3;i<=20;i++)
1346     {
1347       if( ((TH1F*)(fList->FindObject("fHistMcGenerator")))->GetBinContent(i) < 1 )
1348       {
1349         fillVal = i--;
1350         ((TH1F*)(fList->FindObject("fHistMcGenerator")))->GetXaxis()->SetBinLabel(i,Form("%s",genLabel.Data()));
1351         break;
1352       }
1353     }
1354   }
1355
1356   ((TH1F*)(fList->FindObject("fHistMcGenerator")))->Fill(fillVal);
1357
1358   return;
1359 }
1360
1361 //________________________________________________________________________
1362 Bool_t  AliAnalysisTaskPidPidCorrelations::VertexSelection(TObject* obj, Int_t ntracks, Int_t centBin, Double_t gVxMax, Double_t gVyMax, Double_t gVzMax)
1363 {
1364         if (obj->InheritsFrom("AliAODEvent"))
1365         {
1366         Int_t nVertex = ((AliAODEvent*)obj) -> GetNumberOfVertices();
1367
1368                 if ( nVertex > 0 )
1369         {
1370         fMyPrimVtx = (AliAODVertex*)((AliAODEvent*)obj)->GetPrimaryVertex();
1371         if(!fMyPrimVtx) return kFALSE;
1372
1373         Bool_t lReconstructedVertexPresent = kFALSE;
1374         Bool_t lVertexAcceptable = kFALSE;
1375         Double_t fCov[6] = {0.};
1376                 fMyPrimVtx -> GetCovarianceMatrix(fCov);
1377
1378         lReconstructedVertexPresent = ( ( fMyPrimVtx->GetNContributors() > 0 ) && ( fCov[5] != 0 ) );
1379         if(!lReconstructedVertexPresent) { Printf("Bad vertex params"); fHistEventStat -> Fill("Bad vertex params",1); PostData(1,fList); return kFALSE; }
1380
1381         // before vertex cut |zvtx|<10 cm, for fMyPrimVtx only
1382         fHistQA[centBin][0] -> Fill((fMyPrimVtx->GetX()));
1383         fHistQA[centBin][1] -> Fill((fMyPrimVtx->GetY()));
1384         fHistQA[centBin][2] -> Fill((fMyPrimVtx->GetZ()));
1385
1386         // count events having a proper vertex before vertex cut
1387         fHistEventStat -> Fill("Proper vertex",1);
1388         lVertexAcceptable = ( fMyPrimVtx->GetNContributors() < ntracks || (TMath::Abs(fMyPrimVtx->GetX()) < gVxMax && TMath::Abs(fMyPrimVtx->GetY()) < gVyMax && TMath::Abs(fMyPrimVtx->GetZ()) < gVzMax) );
1389         if (!lVertexAcceptable) { Printf("Vertex out of range"); fHistEventStat -> Fill("Bad vertex position",1); PostData(1,fList); return kFALSE; }
1390
1391         // do zvtx binning
1392                 Int_t vtx = GetZvtxBin(fMyPrimVtx->GetZ());
1393                 if (vtx < 0) return kFALSE;
1394
1395         // count events after vertex cut
1396         fHistEventStat -> Fill("Proper vertex within 10cm",1);
1397
1398         // after vertex cut, for fMyPrimVtx only
1399         fHistQA[centBin][3] -> Fill((fMyPrimVtx->GetX()));
1400         fHistQA[centBin][4] -> Fill((fMyPrimVtx->GetY()));
1401         fHistQA[centBin][5] -> Fill((fMyPrimVtx->GetZ()));
1402         //___ end vertex
1403         }
1404         else
1405         {
1406         AliInfo(" Primary-vertex Selection: event REJECTED ...");
1407         return kFALSE;
1408         }
1409         }
1410
1411         //_____ AOD MC , NOT needed...
1412     /*
1413   if (obj->InheritsFrom("AliAODMCHeader"))
1414   {
1415     if (TMath::Abs(((AliAODMCHeader*) obj)->GetVtxZ()) >= fVzMax)
1416     {
1417       AliInfo(" Primary-vertex Selection: event (based on MC) REJECTED ...");
1418       return kFALSE;
1419     }
1420   }
1421 */
1422         return kTRUE;
1423 }
1424
1425 //________________________________________________________________________
1426 void AliAnalysisTaskPidPidCorrelations::CleanUp(TObjArray* tracks, TObject* mcObj)
1427 {
1428   if (!tracks) return;
1429
1430   if (fRemoveWeakDecays)        RemoveWeakDecays(tracks, mcObj);
1431   if (fRemoveDuplicates)        RemoveDuplicates(tracks);
1432 }
1433
1434 //________________________________________________________________________
1435  //________________________________________________________________________
1436 TObjArray* AliAnalysisTaskPidPidCorrelations::AcceptTracks(Int_t centBin, TObject* arrayMC, /*Bool_t onlyprimaries,*/ Bool_t useCuts)
1437 {
1438   //_____ number of RECO AOD tracks
1439   Int_t nTracks = fMyAODEvent->GetNumberOfTracks();
1440 //   AliDebug(5, Form("#tracks= %d %d", nTracks, centBin));
1441 //   AliInfo(Form("There are %d tracks in this event",nTracks));
1442
1443   TObjArray* tracksAccepted = new TObjArray;
1444   tracksAccepted -> SetOwner(kTRUE);
1445
1446   //_____ Track loop
1447   for (Int_t itrk=0; itrk < nTracks; itrk++)
1448   {
1449     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fMyAODEvent->GetTrack(itrk));
1450     if (!track) { AliInfo("AcceptTracks: Could not receive track"); continue; }
1451
1452     /* FIXME
1453     // select only primary tracks
1454     if (onlyprimaries) { if (fMyAODEvent->GetType() != AliAODTrack::kPrimary) continue; }
1455
1456     AliVParticle* vtrack = dynamic_cast<AliVParticle*>(fMyAODEvent->GetTrack(itrk));
1457     //____use only tracks from primaries
1458     if(onlyprimaries)
1459     {
1460       if (vtrack->GetLabel()<=0)        continue;
1461       if (!(static_cast<AliAODMCParticle*>(fMyMcArray->At(vtrack->GetLabel()))->IsPhysicalPrimary()))           continue;
1462     }
1463     delete vtrack; // vtrack needed just for correction for underestimation of strangeness in Monte Carlos, see below, FIXME later...
1464     */
1465     //______track selection cuts, hybrid track cuts
1466
1467     Bool_t bGood = kFALSE;
1468     if          (fFilterType == 0)      bGood = kTRUE;
1469     else if     (fFilterType == 1)      bGood = track -> IsHybridTPCConstrainedGlobal();
1470     else if     (fFilterType == 2)      bGood = track -> IsHybridGlobalConstrainedGlobal();
1471 //     if ( !(track -> TestFilterBit(fFilterBit)) ) { AliInfo(" === TestFilterBit: not passed..."); return 0x0; }
1472 //     if ( (fFilterBit>0) && ( (!track->TestFilterBit(fFilterBit) || (!bGood)) )) return 0;
1473     if ( !bGood )
1474     {
1475 //       AliInfo(" ==== IsHybridGlobalConstrainedGlobal = kFALSE\n");
1476 //       Printf("%s:%d Not matching filter %d/%d %d/%d",(Char_t*)__FILE__,__LINE__,itrk,fMyAODEvent->GetNumberOfTracks(),fFilterBit,track->GetFilterMap());
1477       continue;
1478     }
1479     //if (fTrackStatus != 0 && !CheckTrack(track)) return 0; //clm
1480
1481     Double_t dxy          = 0.;
1482     Double_t dz   = 0.;
1483     Double_t chi2ndf = 0.;
1484     Double_t nCrossedRowsTPC = 0.;
1485     Double_t ratioCrossedRowsOverFindableClustersTPC = 0.;
1486
1487     dxy = track -> DCA();
1488     dz = track -> ZAtDCA();
1489     fHistQA[centBin][6] -> Fill(dxy);
1490     fHistQA[centBin][7] -> Fill(dz);
1491     chi2ndf = track -> Chi2perNDF();
1492     fHistQA[centBin][11] -> Fill(chi2ndf);
1493     nCrossedRowsTPC = track -> GetTPCClusterInfo(2,1);
1494     fHistQA[centBin][12] -> Fill(nCrossedRowsTPC);
1495     //Double_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1496
1497     if (track -> GetTPCNclsF() > 0)
1498     {
1499       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC /1.0/track->GetTPCNclsF();
1500       fHistQA[centBin][13] -> Fill(ratioCrossedRowsOverFindableClustersTPC);
1501     }
1502
1503     //_____ only charged
1504     if (track->Charge()==0) continue;
1505
1506     if (useCuts)
1507     {
1508       if ( track->Eta() < fTrackEtaMin || track->Eta() > fTrackEtaMax || track->Pt() < fTrackPtMin || track->Pt() > fTrackPtMax )
1509         continue;
1510     }
1511
1512     //______ QA
1513     fHistQA[centBin][8] -> Fill(track->Pt());
1514     fHistQA[centBin][9] -> Fill(track->Eta());
1515     fHistQA[centBin][10] -> Fill(track->Phi());
1516
1517     //!FIXME -- correction for underestimation of strangeness in Monte Carlos...
1518     //.....
1519
1520     //_____ MC-RECO matched
1521     if (arrayMC)
1522     {
1523       Int_t label = track->GetLabel();
1524       if (label<0) continue;
1525
1526       //___ check HIJING Bg
1527       Bool_t kIsFromHijingBg = kFALSE;
1528       kIsFromHijingBg = IsFromHijingBg(label);
1529       if (!kIsFromHijingBg) continue;
1530         //_____ redefine track as a matched MC-particle
1531         //track = dynamic_cast<AliAODTrack*>(arrayMC->At(TMath::Abs(label)));
1532         //if (!track) continue;
1533     }
1534
1535     //______ PID
1536     Int_t myPartID = -1;
1537     myPartID = GetParticleID((AliVTrack*)track,centBin,kTRUE);
1538 //     Printf("=========================== myPartID = %d",myPartID);
1539
1540     tracksAccepted -> Add(new AliPidPidCorrelationReducedTrack(myPartID,track->Eta(),track->Phi(),track->Pt(),track->Charge()));
1541     //tracksAccepted -> SetUniqueID(eventno * 100000 + TMath::Abs(part->GetLabel())); //______ !FIXME
1542     //tracksAccepted -> SetUniqueID(part->GetUniqueID()); //______ !FIXME
1543     //tracksAccepted -> Add(part);
1544   }//_____ track loop
1545
1546   return tracksAccepted;
1547 }
1548
1549 //________________________________________________________________________
1550 //________________________________________________________________________
1551 TObjArray* AliAnalysisTaskPidPidCorrelations::AcceptMcTracks(Int_t centBin, Bool_t onlyprimaries, Bool_t useCuts)
1552 {
1553   //_____ number of MC AOD particles
1554   Int_t nMCPart = fMyMcArray->GetEntriesFast();
1555   //if(fVerbose>0) Printf("AOD tracks: %d", nTracks); //! FIXME
1556
1557   TObjArray* tracksAccepted = new TObjArray;
1558   tracksAccepted -> SetOwner(kTRUE); //! FIXME
1559
1560   for (Int_t iMC=0; iMC < nMCPart; iMC++)
1561   {
1562     AliAODMCParticle* partMC = dynamic_cast<AliAODMCParticle*>(fMyMcArray->At(iMC));
1563     if (!partMC) { AliWarning("AcceptMcTracks: Could not receive particle"); continue; }
1564
1565     //_____ eventually only primaries
1566     if (onlyprimaries)
1567     {
1568       if (!partMC->IsPhysicalPrimary()) continue;
1569     }
1570     //_____ eventually only hadrons
1571     Int_t pdgCode = partMC->GetPdgCode();
1572     Bool_t isHadron = TMath::Abs(pdgCode)==211 || TMath::Abs(pdgCode)==2212 || TMath::Abs(pdgCode)==321;
1573     if (!isHadron) continue;
1574     //_____ eventually only charged
1575     if (!partMC->Charge()) continue;
1576
1577     if (useCuts)
1578     {
1579       if ( partMC->Eta() < fTrackEtaMin || partMC->Eta() > fTrackEtaMax || partMC->Pt() < fTrackPtMin || partMC->Pt() > fTrackPtMax )
1580         continue;
1581     }
1582
1583     //______ QA
1584     fHistQA[centBin][8] -> Fill(partMC->Pt());
1585     fHistQA[centBin][9] -> Fill(partMC->Eta());
1586     fHistQA[centBin][10] -> Fill(partMC->Phi());
1587
1588     //______ PID
1589     Int_t myPartID = fPartUndefined;
1590     myPartID = GetParticleIDMC((AliVTrack*)partMC,centBin,kTRUE);
1591 //     Printf("myPartID = %d",myPartID);
1592
1593     tracksAccepted -> Add(new AliPidPidCorrelationReducedTrack(myPartID,partMC->Eta(),partMC->Phi(),partMC->Pt(),partMC->Charge()));
1594     //tracksAccepted -> SetUniqueID(eventno * 100000 + TMath::Abs(part->GetLabel())); //______ !FIXME
1595     //tracksAccepted -> SetUniqueID(part->GetUniqueID()); //______ !FIXME
1596     //tracksAccepted -> Add(part);
1597
1598   }//_____ track loop
1599
1600   return tracksAccepted;
1601 }
1602
1603 //________________________________________________________________________
1604 //________________________________________________________________________
1605 TObjArray* AliAnalysisTaskPidPidCorrelations::AcceptMcRecoMachedTracks(Int_t centBin, Bool_t onlyprimaries, Bool_t useCuts)
1606 {
1607   //_____ number of RECO AOD tracks
1608   Int_t nTracks = fMyAODEvent->GetNumberOfTracks();
1609   //if(fVerbose>0) Printf("AOD tracks: %d", nTracks); //! FIXME
1610
1611   TObjArray* tracksAccepted = new TObjArray;
1612   tracksAccepted -> SetOwner(kFALSE); //! FIXME
1613
1614   //fMyPrimVtx = dynamic_cast<AliAODVertex*>(fMyAODEvent->GetPrimaryVertexSPD());//GetPrimaryVertex()
1615   //Double_t vzAOD = fMyPrimVtx->GetZ();
1616
1617   //______ track loop
1618   for (Int_t itrk=0; itrk < nTracks; itrk++)
1619   {
1620     AliAODTrack* track = dynamic_cast<AliAODTrack*>(fMyAODEvent->GetTrack(itrk));
1621     if (!track) { AliWarning("AcceptMcRecoMachedTracks: Could not receive track"); continue; }
1622
1623     if (track->GetLabel()<0) continue;
1624
1625     AliVParticle* vtrack = dynamic_cast<AliVParticle*>(fMyMcArray->At(track->GetLabel()));
1626     if (!vtrack) { AliWarning("AcceptMcRecoMachedTracks: Could not receive vtrack"); continue; }
1627
1628     //use only tracks from primaries
1629     if (onlyprimaries && !((static_cast<AliAODMCParticle*>(vtrack))->IsPhysicalPrimary())) continue; // !FIXME gives seg fault...
1630
1631     //______track selection cuts, hybrid track cuts
1632     Bool_t bGood = kFALSE;
1633     if  (fFilterType == 0)      bGood = kTRUE;
1634     else if     (fFilterType == 1)      bGood = track -> IsHybridTPCConstrainedGlobal();
1635     else if     (fFilterType == 2)      bGood = track -> IsHybridGlobalConstrainedGlobal();
1636     //if ( !(track -> TestFilterBit(fFilterBit)) ) return 0;
1637     //if ( (fFilterBit>0) && ((!track) -> TestFilterBit(fFilterBit) || (!bGood))) )
1638     if ( !bGood )
1639     {
1640       //AliInfo(" ==== IsHybridGlobalConstrainedGlobal = kFALSE\n");
1641       //Printf("%s:%d Not matching filter %d/%d %d/%d",(Char_t*)__FILE__,__LINE__,it,fMyAODEvent->GetNumberOfTracks(),fFilterBit,track->GetFilterMap());
1642       continue;
1643     }
1644
1645 //     if (fTrackStatus != 0 && !CheckTrack(track)) return 0; //clm
1646 /*
1647     if(vtrack->GetLabel()<0)
1648         continue;
1649     else
1650     {
1651       AliAODMCParticle* partOfTrack = dynamic_cast<AliAODMCParticle*>(fMyMcArray->At(vtrack->GetLabel()));
1652       if(!partOfTrack) Printf("label=%d", vtrack->GetLabel());
1653       delete partOfTrack;delete vtrack;
1654
1655       //correction for underestimation of strangeness in Monte Carlos !FIXME
1656       //...
1657     }
1658 */
1659     //_____ only charged
1660     if (track->Charge()==0) continue;
1661
1662     if (useCuts)
1663     {
1664       if ( track->Eta() < fTrackEtaMin || track->Eta() > fTrackEtaMax || track->Pt() < fTrackPtMin || track->Pt() > fTrackPtMax )
1665         continue;
1666     }
1667
1668     //______ QA
1669     fHistQA[centBin][8] -> Fill(track->Pt());
1670     fHistQA[centBin][9] -> Fill(track->Eta());
1671     fHistQA[centBin][10] -> Fill(track->Phi());
1672
1673     //______ PID
1674     Int_t myPartID = fPartUndefined;
1675     myPartID = GetParticleID(dynamic_cast<AliVTrack*>(track),centBin,kTRUE);
1676     //Printf("myPartID = %d",myPartID);
1677
1678     tracksAccepted -> Add(new AliPidPidCorrelationReducedTrack(myPartID,track->Eta(),track->Phi(),track->Pt(),track->Charge()));
1679     //tracksAccepted -> SetUniqueID(eventno * 100000 + TMath::Abs(part->GetLabel())); //______ !FIXME
1680     //tracksAccepted -> SetUniqueID(part->GetUniqueID()); //______ !FIXME
1681     //tracksAccepted -> Add(part);
1682   }
1683
1684   //Double_t vzMC = fMyMcHeader->GetVtxZ();
1685   return tracksAccepted;
1686 }
1687
1688 //________________________________________________________________________
1689 void AliAnalysisTaskPidPidCorrelations::UserCreateOutputObjects()
1690 {
1691 //   fList = new THashList;
1692   fList = new TList();
1693   fList -> SetOwner(kTRUE);
1694   fList -> SetName(GetOutputListName());
1695
1696   fPIDResponse = (dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager() -> GetInputEventHandler())) -> GetPIDResponse();
1697   if( !fPIDResponse ) { AliError("Cannot get the PID response"); return; }
1698
1699   AddSettingsTree();
1700
1701   const Char_t* kMyPIDTypeName[] = {"ITS","TPC","TOF","HMPID"} ;
1702   const Char_t* kMyParticleSpeciesName[] = { "Electrons","Muons","Pions","Kaons","Protons","Undefined" } ;
1703
1704
1705   fHistNev = new TH1I("fHistNev","Nev.",1,0,1); fHistNev -> Sumw2();
1706   fList -> Add(fHistNev);
1707   fHistTriggerStats = new TH1F("fHistTriggerStats","Trigger statistics;TriggerBit;N_{events}",130,0,130);
1708   fList -> Add(fHistTriggerStats);
1709   fHistTriggerRun = new TH1F("fHistTriggerRun","Trigger;trigger",10,-0.5,9.5);
1710   fList -> Add(fHistTriggerRun);
1711
1712   const Int_t nEventStatBins = 6;
1713   TString gCutName[nEventStatBins] = {"Total","No trigger","Wrong centrality","Bad vertex params","Bad vertex position","Analyzed"};
1714   fHistEventStat = new TH1F("fHistEventStat","Event statistics;;N_{events}",nEventStatBins,-0.5,nEventStatBins-0.5);
1715   for(Int_t i = 1; i <= nEventStatBins; i++) fHistEventStat -> GetXaxis() -> SetBinLabel(i,gCutName[i-1].Data());
1716   fList -> Add(fHistEventStat);
1717
1718   TString gCentName[9] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
1719   fHistCentStats = new TH2F("fHistCentStats","Centrality statistics;;Cent percentile",9,-0.5,8.5,110,-5,105);
1720   for(Int_t i = 1; i <= 9; i++) fHistCentStats->GetXaxis()->SetBinLabel(i,gCentName[i-1].Data());
1721   fList -> Add(fHistCentStats);
1722
1723   fHistCentralityPercentile = new TH1F("fHistCentralityPercentile","Centrality Percentile;Centrality;Entries",101,-1,100);
1724   fList -> Add(fHistCentralityPercentile);
1725   fHistCentralityClass10 = new TH1F("fHistCentralityClass10","Centrality Class 10;Centrality;Entries",10,-0.5,9.5);
1726   fList -> Add(fHistCentralityClass10);
1727   fHistCentralityClass5 = new TH1F("fHistCentralityClass5","Centrality Class 5;Centrality;Entries",20,-0.5,19.5);
1728   fList -> Add(fHistCentralityClass5);
1729
1730   fHistV0M = new TH2F("fHistV0M","V0 Multiplicity C vs. A",5, 0, 20000, 5, 0, 20000);
1731   fList -> Add(fHistV0M);
1732
1733   TString gRefTrackName[6] = {"tracks","tracksPos","tracksNeg","tracksTPConly","clusITS0","clusITS1"};
1734   fHistRefTracks  = new TH2F("fHistRefTracks","Nr of Ref tracks/event vs. ref track estimator;;Nr of tracks",6, 0, 6, 20, 0, 20000);
1735   for(Int_t i = 1; i <= 6; i++) fHistRefTracks->GetXaxis()->SetBinLabel(i,gRefTrackName[i-1].Data());
1736   fList -> Add(fHistRefTracks);
1737
1738   //______ MC QA
1739
1740   if (fUseMC)
1741   {
1742
1743   fHistMcGenerator = new TH1F("fHistMcGenerator","All MC generators;;Entries",20,0,20);
1744   fHistMcGenerator->GetXaxis()->SetBinLabel(1,"Hijing");
1745   fHistMcGenerator->GetXaxis()->SetBinLabel(2,"Pythia");
1746   fHistMcGenerator->GetXaxis()->SetBinLabel(3,"Unknown");
1747   fHistMcGenerator->Sumw2();
1748   fList->Add(fHistMcGenerator);
1749
1750   fHist_HijingBg = new TH1F("fHist_HijingBg","MC generator;;Entries",2,0,2);
1751   fHist_HijingBg->GetXaxis()->SetBinLabel(1,"Hijing BG");
1752   fHist_HijingBg->GetXaxis()->SetBinLabel(2,"No Hijing BG");
1753   fHist_HijingBg->Sumw2();
1754   fList->Add(fHist_HijingBg);
1755
1756   fHistNumOfPartPerEvt = new TH1F("fHistNumOfPartPerEvt","Number of part per event;;Entries",1,0,1);
1757   fHistNumOfPartPerEvt->GetXaxis()->SetBinLabel(1,"# particles per event");
1758   fHistNumOfPartPerEvt->Sumw2();
1759   fList->Add(fHistNumOfPartPerEvt);
1760
1761   fHistMcStats = new TH1F("fHistMcStats","MC stats;;Entries",6,0,6);
1762   fHistMcStats->Sumw2();
1763   fHistMcStats->GetXaxis() -> SetBinLabel(1,"MC all");
1764   fList->Add(fHistMcStats);
1765
1766   fHistMcAllPt = new TH1F("fHistMcAllPt","MC All;pT (GeV/c);Entries",1000,0,100);
1767   fHistMcAllPt->Sumw2();
1768   fList->Add(fHistMcAllPt);
1769   fHistMcAllPt_Hijing= new TH1F("fHistMcAllPt_Hijing","Hijing;pT (GeV/c);Entries",1000,0,100);
1770   fHistMcAllPt_Hijing->Sumw2();
1771   fList->Add(fHistMcAllPt_Hijing);
1772
1773   fHistMcAllPt_Dec= new TH1F("fHistMcAllPt_Dec","Decay;pT (GeV/c);Entries",1000,0,100);
1774   fHistMcAllPt_Dec->Sumw2();
1775   fList->Add(fHistMcAllPt_Dec);
1776   fHistMcAllPt_Inj= new TH1F("fHistMcAllPt_Inj","Injected;pT (GeV/c);Entries",1000,0,100);
1777   fHistMcAllPt_Inj->Sumw2();
1778   fList->Add(fHistMcAllPt_Inj);
1779
1780   fHistMcAllEta_NotHijing= new TH1F("fHistMcAllEta_NotHijing","not Hijing;#eta;Entries",2000,-10,10);
1781   fHistMcAllEta_NotHijing->Sumw2();
1782   fList->Add(fHistMcAllEta_NotHijing);
1783   fHistMcAllEta_Hijing= new TH1F("fHistMcAllEta_Hijing","Hijing;#eta;Entries",2000,-10,10);
1784   fHistMcAllEta_Hijing->Sumw2();
1785   fList->Add(fHistMcAllEta_Hijing);
1786   fHistMcAllEta= new TH1F("fHistMcAllEta","All;#eta;Entries",2000,-10,10);
1787   fHistMcAllEta->Sumw2();
1788   fList->Add(fHistMcAllEta);
1789
1790   fList->Add(new TH1F("pids", ";pdg;tracks", 6001, -3000.5, 3000.5));
1791
1792   }
1793
1794   for (Int_t iCent=0; iCent < fNbinsCent; iCent++)
1795   {
1796     fHistQA[iCent][0] = new TH1F(Form("fHistQAvx_Cent%02d",iCent), Form("Histo Vx All, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1797     fHistQA[iCent][1] = new TH1F(Form("fHistQAvy_Cent%02d",iCent), Form("Histo Vy All, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1798     fHistQA[iCent][2] = new TH1F(Form("fHistQAvz_Cent%02d",iCent), Form("Histo Vz All, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -25., 25.);
1799     fHistQA[iCent][3] = new TH1F(Form("fHistQAvxA_Cent%02d",iCent), Form("Histo Vx  After Cut, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1800     fHistQA[iCent][4] = new TH1F(Form("fHistQAvyA_Cent%02d",iCent), Form("Histo Vy After Cut, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1801     fHistQA[iCent][5] = new TH1F(Form("fHistQAvzA_Cent%02d",iCent), Form("Histo Vz After Cut, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -25., 25.);
1802     fHistQA[iCent][6] = new TH1F(Form("fHistQADcaXyC_Cent%02d",iCent), Form("Histo DCAxy after cut, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1803     fHistQA[iCent][7] = new TH1F(Form("fHistQADcaZC_Cent%02d",iCent), Form("Histo DCAz after cut, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))), 50, -5., 5.);
1804     fHistQA[iCent][8] = new TH1F(Form("fHistQAPt_Cent%02d",iCent),Form("p_{T} distribution, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),900,0.,9.);
1805     fHistQA[iCent][9] = new TH1F(Form("fHistQAEta_Cent%02d",iCent),Form("#eta distribution, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),360,-1.8,1.8);
1806     fHistQA[iCent][10] = new TH1F(Form("fHistQAPhi_Cent%02d",iCent),Form("#phi distribution, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),340,0,6.8);
1807     fHistQA[iCent][11] = new TH1F(Form("fHistQAChi2_Cent%02d",iCent),Form("Chi2 per NDF, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),100,0,10);
1808     fHistQA[iCent][12] = new TH1F(Form("fHistQAnCrossedRowsTPC_Cent%02d",iCent),Form("Number of TPC ccrossed rows, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),200,0,200);
1809     fHistQA[iCent][13] = new TH1F(Form("fHistQAratioCrossedRowsOverFindableClustersTPC_Cent%02d",iCent),Form("Number of TPC ccrossed rows find clusters, CentBin %d-%d",Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),Int_t(fCentAxis -> GetBinUpEdge(iCent+1))),200,0,2);
1810     for(Int_t i=0; i<14; i++) fList -> Add(fHistQA[iCent][i]);
1811
1812     for (Int_t i=0;i<6;i++)
1813     {
1814       fHistRefTracksCent[iCent][i] = new TH2F(Form("fHistRefTracksCent_%02d_%d",iCent,i),
1815                                               Form("Nr of tracksNr of Ref tracks/event in ref track estimator %s for CentBin %d-%d;centrality [%%];Nr of Ref tracks/event",
1816                                                    gRefTrackName[i-1].Data(),
1817                                                    Int_t(fCentAxis->GetBinLowEdge(iCent+1)),
1818                                                    Int_t(fCentAxis->GetBinUpEdge(iCent+1))),10.*(fCentralityPercentileMax-fCentralityPercentileMin),fCentralityPercentileMin,fCentralityPercentileMax,20,0,20000);
1819       fList->Add(fHistRefTracksCent[iCent][i]);
1820     }
1821
1822 /*
1823     for (Int_t iPtBinTrig=0; iPtBinTrig < fNbinsPt; iPtBinTrig++)
1824     {
1825       for (Int_t iPtBinAssoc=0; iPtBinAssoc < fNbinsPt; iPtBinAssoc++)
1826       {
1827         fHistDeltaPhi[iCent][iPtBinTrig][iPtBinAssoc] = new TH1F(Form("fHistDeltaPhi_Cent%02d_PtBin%02d_%02d",iCent,iPtBinTrig,iPtBinAssoc),
1828                                                                Form("%d-%d %%, %3.1f<p_{T,trig}<%3.1f, %3.1f<p_{T,assoc}<%3.1f",
1829                                                                     Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),
1830                                                                     Int_t(fCentAxis -> GetBinUpEdge(iCent+1)),
1831                                                                     fPtAxis -> GetBinLowEdge(iPtBinTrig+1),
1832                                                                     fPtAxis -> GetBinUpEdge(iPtBinTrig+1),
1833                                                                     fPtAxis -> GetBinLowEdge(iPtBinAssoc+1),
1834                                                                     fPtAxis -> GetBinUpEdge(iPtBinAssoc+1)),
1835                                                                     72, -0.5*TMath::Pi(), 1.5*TMath::Pi());
1836
1837         fHistDeltaPhiMix[iCent][iPtBinTrig][iPtBinAssoc] = new TH1F(Form("fHistDeltaPhiMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinTrig,iPtBinAssoc),
1838                                                                   Form("%d-%d %%, %3.1f<p_{T,trig}<%3.1f, %3.1f<p_{T,assoc}<%3.1f MIXED",
1839                                                                        Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),
1840                                                                        Int_t(fCentAxis -> GetBinUpEdge(iCent+1)),
1841                                                                        fPtAxis -> GetBinLowEdge(iPtBinTrig+1),
1842                                                                        fPtAxis -> GetBinUpEdge(iPtBinTrig+1),
1843                                                                        fPtAxis -> GetBinLowEdge(iPtBinAssoc+1),
1844                                                                        fPtAxis -> GetBinUpEdge(iPtBinAssoc+1)),
1845                                                                         72, -0.5*TMath::Pi(), 1.5*TMath::Pi());
1846
1847         fHistDeltaPhi[iCent][iPtBinTrig][iPtBinAssoc]    -> SetXTitle("#Delta#varphi [rad]");
1848         fHistDeltaPhiMix[iCent][iPtBinTrig][iPtBinAssoc] -> SetXTitle("#Delta#varphi [rad]");
1849
1850         fHistDeltaPhi[iCent][iPtBinTrig][iPtBinAssoc]    -> Sumw2();
1851         fHistDeltaPhiMix[iCent][iPtBinTrig][iPtBinAssoc] -> Sumw2();
1852
1853         fList -> Add(fHistDeltaPhi[iCent][iPtBinTrig][iPtBinAssoc]);
1854         fList -> Add(fHistDeltaPhiMix[iCent][iPtBinTrig][iPtBinAssoc]);
1855
1856         fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc] = new TH2F(Form("fHistDphiDeta_Cent%02d_PtBin%02d_%02d",iCent,iPtBinTrig,iPtBinAssoc),
1857                                                                   Form("%d-%d %%, %3.1f<p_{T,trig}<%3.1f, %3.1f<p_{T,assoc}<%3.1f",
1858                                                                        Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),
1859                                                                        Int_t(fCentAxis -> GetBinUpEdge(iCent+1)),
1860                                                                        fPtAxis -> GetBinLowEdge(iPtBinTrig+1),
1861                                                                        fPtAxis -> GetBinUpEdge(iPtBinTrig+1),
1862                                                                        fPtAxis -> GetBinLowEdge(iPtBinAssoc+1),
1863                                                                        fPtAxis -> GetBinUpEdge(iPtBinAssoc+1)),
1864                                                                         72, -0.5*TMath::Pi(), 1.5*TMath::Pi(),
1865                                                                         36, -1.8, 1.8);
1866
1867         fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc] = new TH2F(Form("fHistDphiDetaMix_Cent%02d_PtBin%02d_%02d",iCent,iPtBinTrig,iPtBinAssoc),
1868                                                                      Form("%d-%d %%, %3.1f<p_{T,trig}<%3.1f, %3.1f<p_{T,assoc}<%3.1f MIXED",
1869                                                                           Int_t(fCentAxis -> GetBinLowEdge(iCent+1)),
1870                                                                           Int_t(fCentAxis -> GetBinUpEdge(iCent+1)),
1871                                                                           fPtAxis -> GetBinLowEdge(iPtBinTrig+1),
1872                                                                           fPtAxis -> GetBinUpEdge(iPtBinTrig+1),
1873                                                                           fPtAxis -> GetBinLowEdge(iPtBinAssoc+1),
1874                                                                           fPtAxis -> GetBinUpEdge(iPtBinAssoc+1)),
1875                                                                           72, -0.5*TMath::Pi(), 1.5*TMath::Pi(),
1876                                                                           36, -1.8, 1.8);
1877
1878         fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc]    -> SetXTitle("#Delta#varphi [rad]");
1879         fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc] -> SetXTitle("#Delta#varphi [rad]");
1880         fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc]    -> SetYTitle("#Delta#eta");
1881         fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc] -> SetYTitle("#Delta#eta");
1882
1883         fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc]    -> Sumw2();
1884         fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc] -> Sumw2();
1885
1886         fList -> Add(fHistDphiDeta[iCent][iPtBinTrig][iPtBinAssoc]);
1887         fList -> Add(fHistDphiDetaMix[iCent][iPtBinTrig][iPtBinAssoc]);
1888       }
1889     }*/
1890
1891     fHistTracksEtaTrigVsEtaAssoc[iCent]         = new TH2F(Form("fHistTracksEtaTrigVsEtaAssoc_%02d",iCent),"#eta_{trig} vs #eta_{assoc}",20,-1.,1.,20,-1.,1.);
1892     fList -> Add(fHistTracksEtaTrigVsEtaAssoc[iCent]);
1893     fHistTracksEtaTrigVsEtaAssocMixed[iCent]    = new TH2F(Form("fHistTracksEtaTrigVsEtaAssocMixed_%02d",iCent),"#eta_{trig} vs #eta_{assoc}",20,-1.,1.,20,-1.,1.);
1894     fList -> Add(fHistTracksEtaTrigVsEtaAssocMixed[iCent]);
1895
1896     fHistSingleHadronsPt[iCent]         = new TH1F(Form("fHistSingleHadronsPt_Cent%02d",iCent),"p_{T}", fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray());
1897     fList -> Add(fHistSingleHadronsPt[iCent]);
1898 //     fHistSingleHadronsPt_Mixed[iCent]        = new TH1F(Form("fHistSingleHadronsPt_Mixed_Cent%02d",iCent), "p_{T}", fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray());
1899 //     fList -> Add(fHistSingleHadronsPt_Mixed[iCent]);
1900
1901     fHistSingleHadronsEtaPt[iCent]      = new TH2F(Form("fHistSingleHadronsEtaPt_Cent%02d",iCent), "#eta vs p_{T}",fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray(),fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
1902     fList -> Add(fHistSingleHadronsEtaPt[iCent]);
1903 //     fHistSingleHadronsEtaPt_Mixed[iCent]     = new TH2F(Form("fHistSingleHadronsEtaPt_Mixed_Cent%02d",iCent), "#eta vs p_{T}",fNbinsPt, (Double_t*)fPtAxis->GetXbins()->GetArray(),fEtaAxis->GetNbins(),(Double_t*)fEtaAxis->GetXbins()->GetArray());
1904 //     fList -> Add(fHistSingleHadronsEtaPt_Mixed[iCent]);
1905
1906
1907     fHistTPCdEdx[iCent] = new TH2F(Form("fHistTPCdEdx_%d",iCent), ";p_{TPCin} (GeV/c);dE/dx (au.)", 5, 0., 5., 50, 0., 500.);
1908     fList -> Add(fHistTPCdEdx[iCent]);
1909     fHistTOFbeta[iCent] = new TH2F(Form("fHistTOFbeta_%d",iCent), ";p (GeV/c);v/c", 5, 0., 5., 50, 0.1, 1.1);
1910     fList -> Add(fHistTOFbeta[iCent]);
1911
1912     fHistTPCdEdx_selected[iCent] = new TH2F(Form("fHistTPCdEdx_selected_%d",iCent), ";p_{TPCin} (GeV/c);dE/dx (au.)",  5, 0., 5., 50, 0., 500.);
1913     fList -> Add(fHistTPCdEdx_selected[iCent]);
1914     fHistTOFbeta_selected[iCent] = new TH2F(Form("fHistTOFbeta_selected_%d",iCent), ";p (GeV/c);v/c",  5, 0., 5., 50, 0.1, 1.1);
1915     fList -> Add(fHistTOFbeta_selected[iCent]);
1916
1917     for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
1918     {
1919       fHistNSigmaTPC[iCent][iSpecies] = new TH2F(Form("fHistNSigmaTPC_%d_%s",iCent,AliPID::ParticleName(iSpecies)), Form(";p_{T} (GeV/c); N_{#sigma-TPC}^{%s}", AliPID::ParticleLatexName(iSpecies)), 5, 0., 5., 20, -10., 10.);
1920       fList -> Add(fHistNSigmaTPC[iCent][iSpecies]);
1921       fHistNSigmaTOF[iCent][iSpecies] = new TH2F(Form("fHistNSigmaTOF_%d_%s",iCent,AliPID::ParticleName(iSpecies)), Form(";p_{T} (GeV/c); N_{#sigma-TOF}^{%s}", AliPID::ParticleLatexName(iSpecies)), 5, 0., 5., 20, -10., 10.);
1922       fList -> Add(fHistNSigmaTOF[iCent][iSpecies]);
1923     }
1924
1925
1926     Double_t ptmin[] = { 0.,0.,0.,0. };
1927     Double_t ptmax[] = { 2.,3.,5.,6. };
1928
1929     //_____ nsigma plot
1930     for(Int_t ipart=0;ipart < AliPID::kSPECIES;ipart++)
1931     {
1932       for(Int_t ipid=0;ipid < kMyNSigmaPIDType;ipid++)
1933       {
1934         Double_t miny = -10.;
1935         Double_t maxy =  10.;
1936
1937         //if (ipid == kMyNSigmaTPCTOF) { miny=0; maxy=50; }
1938         fHistoNSigma[iCent] = new TH2F(Form("NSigma_Cent%02d_%d_%d",iCent,ipart,ipid),
1939                                        Form("n#sigma %d-%d %% %s %s",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1)),kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]),
1940                                        10.*(ptmax[ipid]-ptmin[ipid]),ptmin[ipid],ptmax[ipid],2000,miny,maxy);
1941         fHistoNSigma[iCent] -> GetXaxis() -> SetTitle("p_{T} (GeV / c)");
1942         fHistoNSigma[iCent] -> GetYaxis() -> SetTitle(Form("n#sigma %s %s",kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]));
1943         fList -> Add(fHistoNSigma[iCent]);
1944       }
1945     }
1946
1947     if (fUseMC)
1948     {
1949       //_____ nsigmaRec plot
1950       for(Int_t ipart=0;ipart < AliPID::kSPECIES;ipart++)
1951       {
1952         for(Int_t ipid=0;ipid < kMyNSigmaPIDType;ipid++)
1953         {
1954           Double_t miny =-10.;
1955           Double_t maxy = 10.;
1956           //if(ipid == kMyNSigmaTPCTOF) { miny=0; maxy=20; }
1957         fHistoNSigma[iCent] = new TH2F(Form("NSigmaRec_Cent%02d_%d_%d",iCent,ipart,ipid),
1958                                        Form("n#sigma reconstructed %d-%d %% %s %s",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1)),kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]),
1959                                        10.*(ptmax[ipid]-ptmin[ipid]),ptmin[ipid],ptmax[ipid],2000,miny,maxy);
1960         fHistoNSigma[iCent] -> GetXaxis() -> SetTitle("p_{T} (GeV / c)");
1961         fHistoNSigma[iCent] -> GetYaxis() -> SetTitle(Form("n#sigma %s %s",kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]));
1962         fList -> Add(fHistoNSigma[iCent]);
1963         }
1964       }
1965
1966       //_____ nsigmaMC plot
1967       for(Int_t ipart=0;ipart<fPartNSpeciesQA;ipart++)
1968       {
1969         for(Int_t ipid=0;ipid < kMyNSigmaPIDType;ipid++)
1970         {
1971           Double_t miny = -10.;
1972           Double_t maxy =  10.;
1973           //if(ipid == kMyNSigmaTPCTOF) { miny=0; maxy=50; }
1974         fHistoNSigma[iCent] = new TH2F(Form("NSigmaMC_Cent%02d_%d_%d",iCent,ipart,ipid),
1975                                        Form("n#sigma %d-%d %% %s %s",Int_t(fCentAxis->GetBinLowEdge(iCent+1)),Int_t(fCentAxis->GetBinUpEdge(iCent+1)),kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]),
1976                                        10.*(ptmax[ipid]-ptmin[ipid]),ptmin[ipid],ptmax[ipid],2000,miny,maxy);
1977         fHistoNSigma[iCent] -> GetXaxis() -> SetTitle("p_{T} (GeV / c)");
1978         fHistoNSigma[iCent] -> GetYaxis() -> SetTitle(Form("n#sigma %s %s",kMyParticleSpeciesName[ipart],kMyPIDTypeName[ipid]));
1979         fList -> Add(fHistoNSigma[iCent]);
1980         }
1981       }
1982     }// MC
1983   }// iCent
1984
1985
1986   fHistControlConvResoncances = new TH2F("fHistControlConvResoncances", ";id;delta mass", 3, -0.5, 2.5, 20, -0.1, 0.1);
1987   fList -> Add(fHistControlConvResoncances);
1988
1989 /*
1990   Int_t binsH[4]   = { fNbinsZvtx               ,  fNbinsCent                   , fPoolSize     ,  fMinNumTrack };
1991   Double_t minH[4] = { -10              ,  fCentralityPercentileMin     , 0             ,  0            };
1992   Double_t maxH[4] = { 10               ,  fCentralityPercentileMax     , fPoolSize     ,  fMinNumTrack };
1993
1994   fHistPoolMgrQA = new THnSparseD("fHistPoolMgrQA","vtxz:cent:Nevents:Ntracks",4,binsH,minH,maxH);
1995   fHistPoolMgrQA->Sumw2();
1996   fList -> Add(fHistPoolMgrQA);*/
1997
1998         //_______ FillCorrelation _______
1999
2000 //   Int_t anaSteps = 3;       // analysis steps
2001
2002         // single particle histograms
2003 //   Int_t iBinSingle[kTrackVariablesSingle];        // binning for track variables
2004 //   Double_t* dBinsSingle[kTrackVariablesSingle];   // bins for track variables
2005 //   TString axisTitleSingle[kTrackVariablesSingle]; // axis titles for track variables
2006
2007   // two particle histograms
2008   Int_t iBinPair[kTrackVariablesPair];         // binning for track variables
2009   Double_t* dBinsPair[kTrackVariablesPair];    // bins for track variables
2010 //   TString axisTitlePair[kTrackVariablesPair];  // axis titles for track variables
2011
2012   //_____ dim 1,2 -- trig, assoc PID bins
2013 //   const Int_t kNpidBins = 6;
2014 //   Double_t pidBins[kNpidBins+1] = {1.,2.,3.,4.,5.,6.,7.};
2015 //   iBinSingle[0]              = kNpidBins;
2016 //   dBinsSingle[0]             = pidBins;
2017 //   axisTitleSingle[0]  = "PID_{trig.}";
2018 //   iBinPair[0]                = kNpidBins;
2019 //   dBinsPair[0]               = pidBins;
2020 //   axisTitlePair[0]    = "PID_{trig.}";
2021 //   iBinPair[1]                = kNpidBins;
2022 //   dBinsPair[1]               = pidBins;
2023 //   axisTitlePair[1]    = "PID_{assoc.}";
2024   //_____ dim 3,4 -- pT,trig; pT,assoc bins
2025   //const Int_t kNPtBins = 16;
2026   //Double_t ptBins[kNPtBins+1] = {0.2,0.6,1.0,1.5,2.0,2.5,3.0,3.5,4.0,5.0,6.0,7.0,8.0,10.,12.,15.,20.};
2027 //   iBinSingle[1]     = fNbinsPt;//kNPtBins;
2028 //   dBinsSingle[1]    = (Double_t*)fPtAxis->GetXbins()->GetArray();//ptBins;
2029 //   axisTitleSingle[1]  = "p_{T,trig.} (GeV/c)";
2030   iBinPair[0]       = fNbinsPt;//kNPtBins;
2031   dBinsPair[0]      = (Double_t*)fPtAxis->GetXbins()->GetArray();//ptBins;
2032 //   axisTitlePair[2]    = "p_{T,trig.} (GeV/c)";
2033   iBinPair[1]       = fNbinsPt;//kNPtBins;
2034   dBinsPair[1]      = (Double_t*)fPtAxis->GetXbins()->GetArray();//ptBins;
2035 //   axisTitlePair[3]    = "p_{T,assoc.} (GeV/c)";
2036   //_____ dim 5 -- dphi bins
2037   const Int_t kNDeltaPhiBins = 72;
2038   Double_t deltaPhiBins[kNDeltaPhiBins+1];
2039   for(Int_t i = 0; i < kNDeltaPhiBins+1; i++) { deltaPhiBins[i] = -TMath::Pi()/2. + i * 5.*TMath::Pi()/180.; }
2040   iBinPair[2]       = kNDeltaPhiBins;
2041   dBinsPair[2]      = deltaPhiBins;
2042 //   axisTitlePair[4]  = "#Delta#varphi (rad)";
2043   //_____ dim 6 -- deta bins
2044   const Int_t kNDeltaEtaBins = 80;
2045   Double_t deltaEtaBins[kNDeltaEtaBins+1];
2046   for(Int_t i=0; i < kNDeltaEtaBins+1; i++) { deltaEtaBins[i] = - fDeltaEtaMax + i * 2 * fDeltaEtaMax / (Double_t)kNDeltaEtaBins; }
2047   iBinPair[3]       = kNDeltaEtaBins;
2048   dBinsPair[3]      = deltaEtaBins;
2049 //   axisTitlePair[5]   = "#Delta#eta";
2050   //_____ dim 7 -- vertex z bins !! same as in AliEventPoolManager()
2051   //const Int_t kNVertexZBins  = 10;
2052   //Double_t vertexZBins[kNZvtxBins+1] = { -10,   -8,  -6,  -4,  -2,   0,   2,   4,   6,   8,  10 };
2053   //const Int_t kNVertexZBins = 9;
2054   //Double_t vertexZBins[kNVertexZBins+1] = {-10., -7., -5., -3., -1., 1., 3., 5., 7., 10.};
2055 //   iBinSingle[2]    = fNbinsZvtx;//kNVertexZBins;
2056 //   dBinsSingle[2]   = (Double_t*)fZvtxAxis->GetXbins()->GetArray();//vertexZBins;
2057   iBinPair[4]      = fNbinsZvtx;//kNVertexZBins;
2058   dBinsPair[4]     = (Double_t*)fZvtxAxis->GetXbins()->GetArray();//vertexZBins;
2059 //   axisTitleSingle[2]  = "v_{z} (cm)";
2060 //   axisTitlePair[6]    = "v_{z} (cm)";
2061   //_____ dim 8 -- centrality bins !! same as in AliEventPoolManager()
2062   //const Int_t kNCentralityBins = 9;
2063   //Double_t centralityBins[kNCentralityBins+1] = {0.,5.,10.,20.,30.,40.,50.,60.,70.,80.};
2064 //   iBinSingle[3]  = fNbinsCent;//kNCentralityBins;
2065 //   dBinsSingle[3] = (Double_t*)fCentAxis->GetXbins()->GetArray();//centralityBins;
2066 //   axisTitleSingle[3]  = "centrality [#%]";
2067 //   iBinPair[7]    = fNbinsCent;//kNCentralityBins;
2068 //   dBinsPair[7]   = (Double_t*)fCentAxis->GetXbins()->GetArray();//centralityBins;
2069 //   axisTitlePair[7]   = "centrality [#%]";
2070
2071 /*
2072         TString histName;
2073   histName = "fHistCorrSingle";
2074   if (fCentralityEstimator) histName += fCentralityEstimator.Data();
2075   const TString title = histName+";#part;#p_{T,trig};vertex Z (bin);centr (bin);";
2076   fHistCorrSingle = new AliTHn(histName.Data(),title.Data(),anaSteps,kTrackVariablesSingle,iBinSingle);
2077   for (Int_t j=0; j<kTrackVariablesSingle; j++)
2078         {
2079     fHistCorrSingle -> SetBinLimits(j, dBinsSingle[j]);
2080     fHistCorrSingle -> SetVarTitle(j, axisTitleSingle[j]);
2081   }
2082   fList->Add(fHistCorrSingle);
2083   //______ THn for triggers and assoc, pair histo SAME evt
2084   histName = "fHistCorrPairSame";
2085   if (fCentralityEstimator) histName += fCentralityEstimator.Data();
2086   fHistCorrPairSame = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
2087   for (Int_t j=0; j<kTrackVariablesPair; j++)
2088         {
2089     fHistCorrPairSame -> SetBinLimits(j, dBinsPair[j]);
2090     fHistCorrPairSame -> SetVarTitle(j, axisTitlePair[j]);
2091   }
2092   fList->Add(fHistCorrPairSame);
2093   //______ THn for triggers and assoc, pair histo MIXED evt
2094   histName = "fHistCorrPairMixed";
2095   if (fCentralityEstimator) histName += fCentralityEstimator.Data();
2096   fHistCorrPairMixed = new AliTHn(histName.Data(),histName.Data(),anaSteps, kTrackVariablesPair, iBinPair);
2097   for (Int_t j=0; j<kTrackVariablesPair; j++)
2098         {
2099     fHistCorrPairMixed -> SetBinLimits(j, dBinsPair[j]);
2100     fHistCorrPairMixed -> SetVarTitle(j, axisTitlePair[j]);
2101   }
2102   fList->Add(fHistCorrPairMixed);
2103
2104   AliInfo("Finished setting up the AliTHn");*/
2105
2106   //_______ THnSparse histos
2107
2108   Int_t ptMin   = fPtAxis->GetXmin();
2109   Int_t ptMax   = fPtAxis->GetXmax();
2110   Int_t zvtxMin = fZvtxAxis->GetXmin();
2111   Int_t zvtxMax = fZvtxAxis->GetXmax();
2112 //   Int_t centMin      = fCentAxis->GetXmin();
2113 //   Int_t centMax      = fCentAxis->GetXmax();
2114
2115   //_____ single histo
2116 /*
2117   Int_t binsHisto4[kTrackVariablesSingle]   = { iBinSingle[0]   ,  iBinSingle[1]        , iBinSingle[2]  ,  iBinSingle[3]       , 2};
2118   Double_t minHisto4[kTrackVariablesSingle] = { 1                       ,  ptMin                , zvtxMin        ,  centMin             , 7};
2119   Double_t maxHisto4[kTrackVariablesSingle] = { 8               ,  ptMax        , zvtxMax        ,  centMax             , 9};
2120 */
2121   //_____ pair histo
2122 //   Int_t binsHistoPair[kTrackVariablesPair]   = { iBinPair[0] , iBinPair[1]   , iBinPair[2]   , iBinPair[3]   ,       iBinPair[4]     ,       iBinPair[5]     ,       iBinPair[6]     , iBinPair[7]   , 2,2};
2123 //   Double_t minHistoPair[kTrackVariablesPair] = { 1           , 1             , ptMin         , ptMin         ,       -0.5*TMath::Pi()        ,       -fDeltaEtaMax   ,       zvtxMin         , centMin       , 7,7};
2124 //   Double_t maxHistoPair[kTrackVariablesPair] = { 8           , 8             , ptMax         , ptMax         ,       1.5*TMath::Pi() ,       fDeltaEtaMax    ,       zvtxMax         , centMax       , 9,9};
2125
2126 //   Int_t binsHistoPair[kTrackVariablesPair]   = { iBinPair[0] , iBinPair[1]   , iBinPair[2]   , iBinPair[3]   ,       iBinPair[4]     ,       iBinPair[5]     ,       iBinPair[6]     , 2,2};
2127 //   Double_t minHistoPair[kTrackVariablesPair] = { 1           , 1             , ptMin         , ptMin         ,       -0.5*TMath::Pi()        ,       -fDeltaEtaMax   ,       zvtxMin         , 7,7};
2128 //   Double_t maxHistoPair[kTrackVariablesPair] = { 7           , 7             , ptMax         , ptMax         ,       1.5*TMath::Pi() ,       fDeltaEtaMax    ,       zvtxMax         , 9,9};
2129
2130 //   Int_t binsHistoPair[kTrackVariablesPair]   = { iBinPair[0] , iBinPair[1]   , iBinPair[2]   ,       iBinPair[3]     ,       iBinPair[4]     ,       iBinPair[5]     , 2};
2131 //   Double_t minHistoPair[kTrackVariablesPair] = { 1           , ptMin         , ptMin         ,       -0.5*TMath::Pi()        ,       -fDeltaEtaMax   ,       zvtxMin         , 7};
2132 //   Double_t maxHistoPair[kTrackVariablesPair] = { 7           , ptMax         , ptMax         ,       1.5*TMath::Pi() ,       fDeltaEtaMax    ,       zvtxMax         , 9};
2133
2134
2135
2136   Int_t binsHistoPair[kTrackVariablesPair]   = { iBinPair[0]    , iBinPair[1]   ,       iBinPair[2]     ,       iBinPair[3]     ,       iBinPair[4]     };
2137   Double_t minHistoPair[kTrackVariablesPair] = { static_cast<Double_t>(ptMin)           , static_cast<Double_t>(ptMin)          ,       -0.5*TMath::Pi()        ,       static_cast<Double_t>(-fDeltaEtaMax)    ,       static_cast<Double_t>(zvtxMin)          };
2138   Double_t maxHistoPair[kTrackVariablesPair] = { static_cast<Double_t>(ptMax)           , static_cast<Double_t>(ptMax)          ,       1.5*TMath::Pi() ,       static_cast<Double_t>(fDeltaEtaMax)     ,       static_cast<Double_t>(zvtxMax)          };
2139
2140
2141 /*
2142   Double_t ptMin = 0.0, ptMax = 6.0;
2143 //   Double_t binsPt[]  = {0.0, 0.1, 0.2, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0};
2144   Double_t binsPt[]  = {0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0};
2145   const Int_t nPtBins = length(binsPt)-1;
2146
2147   Double_t zvtxMin = -10., zvtxMax = 10.;
2148   Double_t binsZvtx[] = {-10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.};
2149   const Int_t nZvtxBins = length(binsZvtx)-1;
2150
2151   Double_t centMin = 0., centMax = 10.;
2152 //   Double_t binsCent[] = {0.,1.,2.,4.,5.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.1}; // 90
2153   Double_t binsCent[] = {0.,1.,2.,5.,10.}; // 90
2154   const Int_t nCentBins = length(binsCent)-1;
2155
2156
2157   //_____ single histo, 4 dim matrix
2158   Int_t binsHisto4[4+1]   = { 6 ,       nPtBins , nZvtxBins     , nCentBins     ,2};
2159   Double_t minHisto4[4+1] = { 1 ,       ptMin   , zvtxMin       , centMin       ,7};
2160   Double_t maxHisto4[4+1] = { 7 ,       ptMax   , zvtxMax       , centMax       ,9};
2161
2162   //_____ pair histo, 8 dim matrix
2163   Int_t binsHistoPair[8+2]   = { 6 , 6   , nPtBins      , nPtBins ,                  72         ,  Int_t(fDeltaEtaMax*10)       ,  nZvtxBins    , nCentBins     ,2,2};
2164   Double_t minHistoPair[8+2] = { 1 , 1   , ptMin                , ptMin   ,     -0.5*TMath::Pi()        , -fDeltaEtaMax                 ,  zvtxMin      , centMin       ,7,7};
2165   Double_t maxHistoPair[8+2] = { 7 , 7   , ptMax                , ptMax   ,      1.5*TMath::Pi()        ,  fDeltaEtaMax                 ,  zvtxMax      , centMax       ,9,9};
2166 */
2167 /*
2168   fHistCorrSingle = new THnSparseD("fHistCorrSingle","pdg_trig:pt_trig:Zvtx:centr:chhdr_trig",4+1,binsHisto4,minHisto4,maxHisto4);
2169   fHistCorrSingle->GetAxis(0)->SetTitle("partType_{trig}");
2170   fHistCorrSingle->SetBinEdges(1,dBinsSingle[1]);
2171 //   fHistCorrSingle->SetBinEdges(1,binsPt);
2172   fHistCorrSingle->GetAxis(1)->SetTitle("p_{T, trig} (GeV/c)");
2173   fHistCorrSingle->SetBinEdges(2,dBinsSingle[2]);
2174 //   fHistCorrSingle->SetBinEdges(2,binsZvtx);
2175   fHistCorrSingle->GetAxis(2)->SetTitle("v_{z} (cm)");
2176   fHistCorrSingle->SetBinEdges(3,dBinsSingle[3]);
2177 //   fHistCorrSingle->SetBinEdges(3,binsCent);
2178   fHistCorrSingle->GetAxis(3)->SetTitle("Centrality [#%]");
2179   fHistCorrSingle->GetAxis(4)->SetTitle("hadron_trig");
2180
2181   fHistCorrSingle->Sumw2();
2182   fList -> Add(fHistCorrSingle);
2183 */
2184
2185   for(Int_t i=0;i<2;i++)
2186   {
2187     fHistCorrPair[i] = new THnSparseD(Form("fHistCorrPair_%d",i),"pt_trig:pt_assoc:DeltaPhi:DeltaEta:Zvtx",5,binsHistoPair,minHistoPair,maxHistoPair);
2188     fHistCorrPair[i]->SetBinEdges(0,dBinsPair[0]);
2189     fHistCorrPair[i]->SetBinEdges(1,dBinsPair[1]);
2190     fHistCorrPair[i]->GetAxis(0)->SetTitle("p_{T, trig} (GeV/c)");
2191     fHistCorrPair[i]->GetAxis(1)->SetTitle("p_{T, assoc} (GeV/c)");
2192     fHistCorrPair[i]->SetBinEdges(2,dBinsPair[2]);
2193     fHistCorrPair[i]->SetBinEdges(3,dBinsPair[3]);
2194     fHistCorrPair[i]->GetAxis(2)->SetTitle("#Delta #varphi");
2195     fHistCorrPair[i]->GetAxis(3)->SetTitle("#Delta #eta");
2196     fHistCorrPair[i]->SetBinEdges(4,dBinsPair[4]);
2197     fHistCorrPair[i]->GetAxis(4)->SetTitle("v_{z} (cm)");
2198     fHistCorrPair[i]->Sumw2();
2199     fList -> Add(fHistCorrPair[i]);
2200   }
2201
2202   //     TAxis* axis = fHistCorrPairSame->GetGrid(0)->GetGrid()->GetAxis(2);
2203   fTriggerWeighting = fHistCorrPair[0]->GetAxis(0);
2204   fHistTriggerWeighting = new TH1D("fHistTriggerWeighting","", fTriggerWeighting->GetNbins(),fTriggerWeighting->GetXbins()->GetArray());
2205
2206
2207
2208 /*
2209   for(Int_t i=0;i<2;i++)
2210   {
2211     fHistCorrPair[i] = new THnSparseD(Form("fHistCorrPair_%d",i),"pdg_trig:pt_trig:pt_assoc:DeltaPhi:DeltaEta:Zvtx:chhdr_assoc",7,binsHistoPair,minHistoPair,maxHistoPair);
2212     fHistCorrPair[i]->GetAxis(0)->SetTitle("partType_{trig}");
2213 //     fHistCorrPair[i]->GetAxis(1)->SetTitle("partType_{assoc}");
2214     fHistCorrPair[i]->SetBinEdges(0,dBinsPair[0]);
2215 //     fHistCorrPair[i]->SetBinEdges(1,dBinsPair[1]);
2216 //     fHistCorrPair[i]->SetBinEdges(2,binsPt);
2217 //     fHistCorrPair[i]->SetBinEdges(3,binsPt);
2218     fHistCorrPair[i]->GetAxis(1)->SetTitle("p_{T, trig} (GeV/c)");
2219     fHistCorrPair[i]->GetAxis(2)->SetTitle("p_{T, assoc} (GeV/c)");
2220     fHistCorrPair[i]->SetBinEdges(1,dBinsPair[1]);
2221     fHistCorrPair[i]->SetBinEdges(2,dBinsPair[2]);
2222     fHistCorrPair[i]->GetAxis(3)->SetTitle("#Delta #varphi");
2223     fHistCorrPair[i]->GetAxis(4)->SetTitle("#Delta #eta");
2224     fHistCorrPair[i]->SetBinEdges(3,dBinsPair[3]);
2225     fHistCorrPair[i]->SetBinEdges(4,dBinsPair[4]);
2226 //     fHistCorrPair[i]->SetBinEdges(6,binsZvtx);
2227 //     fHistCorrPair[i]->SetBinEdges(7,binsCent);
2228     fHistCorrPair[i]->GetAxis(5)->SetTitle("v_{z} (cm)");
2229     fHistCorrPair[i]->SetBinEdges(5,dBinsPair[5]);
2230 //     fHistCorrPair[i]->GetAxis(7)->SetTitle("hadron_trig");
2231     fHistCorrPair[i]->GetAxis(6)->SetTitle("hadron_assoc");
2232
2233     fHistCorrPair[i]->Sumw2();
2234     fList -> Add(fHistCorrPair[i]);
2235   }
2236 */
2237 /*
2238   if( fWriteCorrTree == kTRUE )
2239   {
2240     fVariablesTreeCorr = new TTree("CorrTree","2PartCorrPID tree");
2241     Int_t nVar = 5;
2242     fCorrVariables = new Float_t [nVar];
2243     TString * fCorrVariableNames = new TString[nVar];
2244     fCorrVariableNames[0] = "fMyPIDtrig";
2245     fCorrVariableNames[1] = "fMyPTtrig";
2246     fCorrVariableNames[2] = "fMyZvtx";
2247     fCorrVariableNames[3] = "fMyCent";
2248     fCorrVariableNames[4] = "fMyChardHdr";
2249
2250     for(Int_t ivar=0; ivar<nVar; ivar++)
2251       fVariablesTreeCorr->Branch(fCorrVariableNames[ivar].Data(),&fCorrVariables[ivar],Form("%s/f",fCorrVariableNames[ivar].Data()));
2252     PostData(3,fVariablesTreeCorr);
2253   }*/
2254
2255   for (Int_t i=0; i<fNMaxBinsPt; i++)
2256   {
2257     fHistTwoTrackDistancePt[i][0] = new TH2F(Form("fHistTwoTrackDistancePt[%d][0]",i),Form("fHistTwoTrackDistancePt[%d][0] -- PtBin %.2f-%.2f;#Delta#eta;#Delta#varphi^{*}_{min}",i,fPtAxis->GetBinLowEdge(i+1),fPtAxis->GetBinUpEdge(i+1)), 100, -0.05, 0.05, 400, -0.2, 0.2);
2258     fList->Add(fHistTwoTrackDistancePt[i][0]);
2259     fHistTwoTrackDistancePt[i][1] = (TH2F*)fHistTwoTrackDistancePt[i][0]->Clone(Form("fHistTwoTrackDistancePt[%d][1]",i));
2260     fList->Add(fHistTwoTrackDistancePt[i][1]);
2261   }
2262
2263   // QA
2264   fHistHBTbefore        = new TH2F("fHistHBTbefore","before HBT cut",20,0,2,20,0,2.*TMath::Pi()); fList->Add(fHistHBTbefore);
2265   fHistHBTafter         = new TH2F("fHistHBTafter","after HBT cut",20,0,2,20,0,2.*TMath::Pi()); fList->Add(fHistHBTafter);
2266 /*
2267   fHistConversionbefore = new TH2F("fHistConversionbefore","before Conversion cut",200,0,2,200,0,2.*TMath::Pi());
2268   fHistConversionafter  = new TH2F("fHistConversionafter","after Conversion cut",200,0,2,200,0,2.*TMath::Pi());
2269   fHistPsiMinusPhi      = new TH2F("fHistPsiMinusPhi","",4,-0.5,3.5,100,0,2.*TMath::Pi());
2270   fHistResonancesBefore = new TH3D("fHistResonancesBefore","before resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2271   fHistResonancesRho    = new TH3D("fHistResonancesRho","after #rho resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2272   fHistResonancesK0     = new TH3D("fHistResonancesK0","after #rho, K0 resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2273   fHistResonancesLambda = new TH3D("fHistResonancesLambda","after #rho, K0, Lambda resonance cut;#Delta#eta;#Delta#phi;M_{inv}",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2274   fHistQbefore          = new TH3D("fHistQbefore","before momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2275   fHistQafter           = new TH3D("fHistQafter","after momentum difference cut;#Delta#eta;#Delta#phi;|#Delta p_{T}| (GeV/c)",50,-2.0,2.0,50,-TMath::Pi()/2.,3.*TMath::Pi()/2.,300,0,1.5);
2276 */
2277
2278   //_____ PoolManager setup for mixing
2279   SetupForMixing();
2280
2281   //_____ AliCFContainer settings
2282   Int_t nTrackBin[kNvars];
2283   Double_t* trackBins[kNvars];
2284   TString trackAxisTitle[kNvars];
2285
2286   TString defaultBinningStr;
2287   defaultBinningStr =   "eta: -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0\n"
2288 //                      "phi: 0., 0.62831, 1.25662, 1.88493, 2.51324, 3.14156, 3.76987, 4.39818, 5.02649, 5.65480, 6.28312\n"
2289                         "p_T: 0.0, 0.5, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0\n"
2290                         "PID: 0., 1., 2., 3., 4., 5., 6., 7., 8., 9.\n"
2291                         "centrality: 0.0, 5.0\n"
2292                         "vertex: -10.,-8.,-6.,-4.,-2.,0.,2.,4.,6.,8.,10.\n";
2293  // =========================================================
2294   // Customization (adopted from AliUEHistograms)
2295   // =========================================================
2296   TObjArray* lines = defaultBinningStr.Tokenize("\n");
2297   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2298   {
2299     TString line(lines->At(i)->GetName());
2300     TString tag = line(0, line.Index(":")+1);
2301     if (!fCustomBinning.BeginsWith(tag) && !fCustomBinning.Contains(TString("\n") + tag))
2302       fBinningString += line + "\n";
2303     else
2304       AliInfo(Form("Using custom binning for %s", tag.Data()));
2305   }
2306   delete lines;
2307   fBinningString += fCustomBinning;
2308
2309   AliInfo(Form("Used THn Binning:\n%s",fBinningString.Data()));
2310
2311
2312   // eta
2313   trackBins[0] = GetBinning(fBinningString, "eta", nTrackBin[0]);
2314   trackAxisTitle[0] = "#eta";
2315 //   Int_t nEtaBins = 20;
2316 //   Double_t etaMin = -2., etaMax = -2.;
2317   // phi
2318 //   trackBins[1] = GetBinning(fBinningString, "phi", nTrackBin[1]);
2319 //   trackAxisTitle[1] = "#phi (rad)";
2320 //   Int_t nPhiBins = 36;
2321 //   Double_t phiMin = 0.; Double_t phiMax = 2.*TMath::Pi();
2322   // pT
2323   trackBins[1] = GetBinning(fBinningString, "p_T", nTrackBin[1]);
2324   trackAxisTitle[1] = "p_{T} (GeV/c)";
2325 //   Int_t nPtBins = 20;
2326 //   Double_t ptMin = 0., ptMax = 6.
2327   // particle species
2328   trackBins[2] = GetBinning(fBinningString, "PID", nTrackBin[2]);
2329   trackAxisTitle[2] = "particle species";
2330 //   Int_t nPidBins = 9;
2331 //   Double_t pidMin = -0.5, pidMax = 8.5;
2332   // centrality
2333   trackBins[3] = GetBinning(fBinningString, "centrality", nTrackBin[3]);
2334   trackAxisTitle[3] = "cent (%)";
2335 //   Double_t centMin = 0., centMax = 100.1;
2336 //   Double_t binsCent[] = {0.,1.,2.,4.,5.,10.,20.,30.,40.,50.,60.,70.,80.,90.,100.1};
2337 //   Int_t nCentBins = length(binsCent)-1;
2338   // vtx-z axis
2339   trackBins[4] = GetBinning(fBinningString, "vertex", nTrackBin[4]);
2340   trackAxisTitle[4] = "z-vtx (cm)";
2341 //   Int_t nVzBins = 40;
2342 //   Double_t vzMin = -10., vzMax = 10.;
2343
2344 //   Int_t nbins[kNvars] = {nEtaBins, nPhiBins, nPtBins, nPidBins, nCentBins, nVzBins};
2345 //   Double_t xmin[kNvars] = {etaMin, phiMin, ptMin, pidMin, centMin, vzMin};
2346 //   Double_t xmax[kNvars] = {etaMax, phiMax, ptMax, pidMax, centMax, vzMax};
2347
2348   fMyCFCont = new AliCFContainer("PidPidCorrCFCont","CF for PidPid Corr",kNsteps,kNvars,nTrackBin);
2349 //   fMyCFCont = new AliCFContainer("PidPidCorrCFCont","CF for PidPid Corr",kNsteps,kNvars,nbins);
2350
2351   for ( Int_t idim=0; idim < kNvars; idim++)
2352   {
2353     fMyCFCont->SetVarTitle(idim, Form("%s",trackAxisTitle[idim].Data()));
2354     fMyCFCont->SetBinLimits(idim, trackBins[idim]);
2355 //     fMyCFCont->SetBinLimits(idim, xmin[idim], xmax[idim]);
2356   }
2357
2358   TString stepTitle[kNsteps] = {"generated","reco","recomatch"};
2359   for (Int_t istep=0; istep<kNsteps; istep++)
2360     fMyCFCont->SetStepTitle(istep, stepTitle[istep].Data());
2361
2362   PostData(1, fList);
2363   PostData(2, fMyCFCont);
2364 }
2365
2366 //________________________________________________________________________
2367 Double_t* AliAnalysisTaskPidPidCorrelations::GetBinning(const Char_t* configuration, const Char_t* tag, Int_t& nBins)
2368 {
2369   // takes the binning from <configuration> identified by <tag>
2370   // configuration syntax example:
2371   // eta: 2.4, -2.3, -2.2, -2.1, -2.0, -1.9, -1.8, -1.7, -1.6, -1.5, -1.4, -1.3, -1.2, -1.1, -1.0, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4
2372   // phi: .....
2373   //
2374   // returns bin edges which have to be deleted by the caller
2375
2376   TString config(configuration);
2377   TObjArray* lines = config.Tokenize("\n");
2378   for (Int_t i=0; i<lines->GetEntriesFast(); i++)
2379   {
2380     TString line(lines->At(i)->GetName());
2381     if (line.BeginsWith(TString(tag) + ":"))
2382     {
2383       line.Remove(0, strlen(tag) + 1);
2384       line.ReplaceAll(" ", "");
2385       TObjArray* binning = line.Tokenize(",");
2386       Double_t* bins = new Double_t[binning->GetEntriesFast()];
2387       for (Int_t j=0; j<binning->GetEntriesFast(); j++)
2388         bins[j] = TString(binning->At(j)->GetName()).Atof();
2389
2390       nBins = binning->GetEntriesFast() - 1;
2391
2392       delete binning;
2393       delete lines;
2394       return bins;
2395     }
2396   }
2397
2398   delete lines;
2399   AliFatal(Form("Tag %s not found in %s", tag, configuration));
2400   return 0;
2401 }
2402
2403 //________________________________________________________________________
2404 Bool_t AliAnalysisTaskPidPidCorrelations::CheckTrack(AliAODTrack* track)
2405 {
2406   // check if the track status flags are set
2407   UInt_t status = track -> GetStatus();
2408   if ((status & fTrackStatus) == fTrackStatus)
2409     return kTRUE;
2410   return kFALSE;
2411 }
2412
2413 //________________________________________________________________________
2414 void AliAnalysisTaskPidPidCorrelations::CalculateNSigmas(AliAODTrack* track, Int_t centBin, Bool_t* pidFlag, Bool_t fillQA)
2415 {
2416   if(!fPIDResponse)
2417     fPIDResponse = (dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager() -> GetInputEventHandler())) -> GetPIDResponse();
2418   if(!fPIDResponse)
2419     AliFatal("AliAnalysisTaskPidPidCorrelations::CalculateNSigmas(): Cannot get pid response");
2420
2421   // ____________ set cuts on pT, nsigma for ITS, TPC, TOF, HMPID
2422   Double_t fgTPCPIDmomcut[AliPID::kSPECIES]     = { 0., 0., 0.5, 0.5, 0.7 };    // _____ TPC pT
2423   Double_t fgTPCPIDsigmacut[AliPID::kSPECIES]   = { 0., 0., 3. , 3. , 3.  };    // _____ TPC nsigma
2424   Double_t fgTOFPIDmomcut[AliPID::kSPECIES]     = { 0., 0., 1.5, 1.5, 2.0 };    // _____ TOF pT
2425   Double_t fgTOFPIDsigmacut[AliPID::kSPECIES]   = { 0., 0., 2. , 2. , 2.  };    // _____ TOF nsgima
2426 //
2427   Double_t p    = track -> P();
2428   Double_t pt   = track -> Pt();
2429   Double_t ptpc = track -> GetTPCmomentum();
2430
2431   //___ nsigmas
2432 //   Double_t nsigmaITS[AliPID::kSPECIES]       = {-999,-999,-999,-999,-999};   //      ITS
2433 //   Double_t nsigmaTPC[AliPID::kSPECIES]       = {-999,-999,-999,-999,-999};   //      TPC
2434 //   Double_t nsigmaTOF[AliPID::kSPECIES]       = {-999,-999,-999,-999,-999};   //      TOF
2435 //   Double_t nsigmaHMPID[AliPID::kSPECIES]     = {-999,-999,-999,-999,-999};   //      HMPID
2436
2437   //___ auxiliary nsigma for QA
2438   Double_t auxNsigmaTPC[AliPID::kSPECIES]       = {-999,-999,-999,-999,-999};   //      TPC
2439   Double_t auxNsigmaTOF[AliPID::kSPECIES]       = {-999,-999,-999,-999,-999};   //      TOF
2440
2441   Double_t dEdx = MakeTPCPID(track, auxNsigmaTPC);
2442   Double_t beta = MakeTOFPID(track, auxNsigmaTOF);
2443
2444   Bool_t isOnTPCPID = dEdx > 0.;
2445   Bool_t isOnTOFPID = beta > 0.;
2446
2447   Bool_t ITSnSigmaIsOk[AliPID::kSPECIES] = {kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
2448   Bool_t TPCnSigmaIsOk[AliPID::kSPECIES] = {kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
2449   Bool_t TOFnSigmaIsOk[AliPID::kSPECIES] = {kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
2450   Bool_t HMPnSigmaIsOk[AliPID::kSPECIES] = {kFALSE,kFALSE,kFALSE,kFALSE,kFALSE};
2451
2452   for(Int_t ispecies=0; ispecies < AliPID::kSPECIES; ispecies++)
2453   {
2454     //_____ ITS
2455     nsigmaITS[ispecies] = fPIDResponse -> NumberOfSigmasITS((AliVTrack*)track,(AliPID::EParticleType)ispecies);
2456     if (nsigmaITS[ispecies] != -999.) ITSnSigmaIsOk[ispecies] = kTRUE;
2457     //_____ TPC
2458     nsigmaTPC[ispecies] = fPIDResponse -> NumberOfSigmasTPC((AliVTrack*)track,(AliPID::EParticleType)ispecies);
2459     if (nsigmaTPC[ispecies] != -999.) TPCnSigmaIsOk[ispecies] = kTRUE;
2460     //_____ TOF
2461     nsigmaTOF[ispecies] = fPIDResponse -> NumberOfSigmasTOF((AliVTrack*)track,(AliPID::EParticleType)ispecies);
2462     if (nsigmaTOF[ispecies] != -999.) TOFnSigmaIsOk[ispecies] = kTRUE;
2463     //_____ HMPID
2464     nsigmaHMPID[ispecies] = fPIDResponse -> NumberOfSigmasHMPID((AliVTrack*)track,(AliPID::EParticleType)ispecies);
2465     if (nsigmaHMPID[ispecies] != -999.) HMPnSigmaIsOk[ispecies] = kTRUE;
2466   }
2467
2468   for(Int_t ispecies=0; ispecies < AliPID::kSPECIES; ispecies++)
2469   {
2470     if (ITSnSigmaIsOk[ispecies]) fnsigmas[ispecies][0] = nsigmaITS[ispecies];
2471     if (TPCnSigmaIsOk[ispecies]) fnsigmas[ispecies][1] = nsigmaTPC[ispecies];
2472     if (TOFnSigmaIsOk[ispecies]) fnsigmas[ispecies][2] = nsigmaTOF[ispecies];
2473     if (HMPnSigmaIsOk[ispecies]) fnsigmas[ispecies][3] = nsigmaHMPID[ispecies];
2474   }
2475
2476   //______ fill QA
2477   if (fillQA)
2478   {
2479     TH2F* h[fNMaxBinsCentrality];
2480
2481     for (Int_t iSpecies=0; iSpecies < AliPID::kSPECIES; iSpecies++)
2482     {
2483       pidFlag[iSpecies] = kFALSE;
2484
2485       //_______________ PID signal plots for TPC, TOF
2486       if (isOnTPCPID)
2487       {
2488         fHistNSigmaTPC[centBin][iSpecies] -> Fill(pt, nsigmaTPC[iSpecies]);
2489         fHistTPCdEdx[centBin] -> Fill(ptpc, dEdx);
2490       }
2491       if (isOnTOFPID)
2492       {
2493         fHistNSigmaTOF[centBin][iSpecies] -> Fill(pt, nsigmaTOF[iSpecies]);
2494         fHistTOFbeta[centBin] -> Fill(p, beta);
2495       }
2496
2497       //_______________ PID nsigma plots for TPCTOF, TPC, TOF
2498       if (isOnTPCPID && isOnTOFPID)
2499       {
2500         // get TPC & TOF nsigma -- combined
2501         if (pt > fgTOFPIDmomcut[iSpecies]) // cut on TOF pT
2502 //      nsigmaTPCTOF[iSpecies] = TMath::Sqrt(nsigmaTPC[iSpecies]*nsigmaTPC[iSpecies]+nsigmaTOF[iSpecies]*nsigmaTOF[iSpecies]);
2503
2504         if (pt < fgTOFPIDmomcut[iSpecies] && TMath::Abs(nsigmaTOF[iSpecies]) < fgTOFPIDsigmacut[iSpecies] && TMath::Abs(nsigmaTPC[iSpecies]) < fgTPCPIDmomcut[iSpecies])
2505         {
2506           fHistTOFbeta_selected[centBin] -> Fill(p, beta);
2507           pidFlag[iSpecies] = kTRUE;
2508         }
2509       }
2510       //_____ TPC-only PID cuts
2511       else if (isOnTPCPID && !isOnTOFPID)
2512       {
2513         if (pt < fgTPCPIDmomcut[iSpecies] && TMath::Abs(nsigmaTPC[iSpecies]) < fgTPCPIDsigmacut[iSpecies])
2514         {
2515           fHistTPCdEdx_selected[centBin] -> Fill(ptpc, dEdx);
2516           pidFlag[iSpecies] = kTRUE;
2517         }
2518       }
2519       //_____ TOF-only PID cuts
2520       else if (!isOnTPCPID && isOnTOFPID)
2521       {
2522         if (pt < fgTOFPIDmomcut[iSpecies] && TMath::Abs(nsigmaTOF[iSpecies]) < fgTOFPIDsigmacut[iSpecies])
2523         {
2524           fHistTOFbeta_selected[centBin] -> Fill(p, beta);
2525           pidFlag[iSpecies] = kTRUE;
2526         }
2527       }
2528
2529       //__________ fill nsigma vs pT for all the detectors
2530       for (int iPid=0; iPid < kMyNSigmaPIDType; iPid++)
2531       {
2532         if (fUseMC)
2533           h[centBin] = GetHisto2D(Form("NSigmaRec_Cent%02d_%d_%d",centBin,iSpecies,iPid));
2534         else
2535           h[centBin] = GetHisto2D(Form("NSigma_Cent%02d_%d_%d",centBin,iSpecies,iPid));
2536
2537         h[centBin] -> Fill(pt,fnsigmas[iSpecies][iPid]);
2538       }
2539     }//____ iSpecies
2540   }//____ fill QA
2541 }
2542
2543 //___________________________________________________________
2544 Int_t AliAnalysisTaskPidPidCorrelations::FindNSigma(AliAODTrack* track)
2545 {
2546   //_____ ITS
2547   Double_t fgITSPIDmomcutMin[AliPID::kSPECIES]  = { 0., 0., 0.2 , 0.6 , 0.4 };  // pT_min
2548   Double_t fgITSPIDsigmacutMin[AliPID::kSPECIES]        = { 0., 0.,-1.5 ,-2.0 ,-6.0 };  // nsigma_min
2549   Double_t fgITSPIDmomcutMax[AliPID::kSPECIES]  = { 0., 0., 1.2 , 1.4 , 2.0 };  // pT_max
2550   Double_t fgITSPIDsigmacutMax[AliPID::kSPECIES]        = { 0., 0., 1.5 , 1.0 , 1.0 };  // nsigma_max
2551   //_____ TPC
2552   Double_t fgTPCPIDmomcutMin[AliPID::kSPECIES]  = { 0., 0., 0.2 , 0.6 , 0.6 };  // pT_min
2553   Double_t fgTPCPIDsigmacutMin[AliPID::kSPECIES]        = { 0., 0.,-2.0 ,-2.0 ,-4.0 };  // nsigma_min
2554   Double_t fgTPCPIDmomcutMax[AliPID::kSPECIES]  = { 0., 0., 1.0 , 2.0 , 4.0 };  // pT_max
2555   Double_t fgTPCPIDsigmacutMax[AliPID::kSPECIES]        = { 0., 0., 2.0 , 3.0 , 3.0 };  // nsigma_max
2556   //_____ TOF
2557   Double_t fgTOFPIDmomcutMin[AliPID::kSPECIES]  = { 0., 0., 0.4 , 0.6 , 0.6 };  // pT_min
2558   Double_t fgTOFPIDsigmacutMin[AliPID::kSPECIES]        = { 0., 0.,-1.0 ,-2.0 ,-2.0 };  // nsgima_min
2559   Double_t fgTOFPIDmomcutMax[AliPID::kSPECIES]  = { 0., 0., 1.2 , 2.6 , 4.0 };  // pT_max
2560   Double_t fgTOFPIDsigmacutMax[AliPID::kSPECIES]        = { 0., 0., 1.0 , 2.0 , 3.0 };  // nsgima_max
2561   //_____ HMPID
2562   Double_t fgHMPIDPIDmomcutMin[AliPID::kSPECIES]                = { 0., 0., 0.5 , 0.8 , 1.4 };  // pT_min
2563   Double_t fgHMPIDPIDsigmacutMin[AliPID::kSPECIES]      = { 0., 0.,-2.0 ,-2.0 ,-2.0 };  // nsigma_min
2564   Double_t fgHMPIDPIDmomcutMax[AliPID::kSPECIES]                = { 0., 0., 1.5 , 1.5 , 4.0 };  // pT_max
2565   Double_t fgHMPIDPIDsigmacutMax[AliPID::kSPECIES]      = { 0., 0., 3.0 , 6.0 , 4.0 };  // nsigma_max
2566
2567   Bool_t isITSok[AliPID::kSPECIES]      =  {kFALSE};
2568   Bool_t isTPCok[AliPID::kSPECIES]      =  {kFALSE};
2569   Bool_t isTOFok[AliPID::kSPECIES]      =  {kFALSE};
2570   Bool_t isHMPIDok[AliPID::kSPECIES]    =  {kFALSE};
2571
2572   //Int_t trackCharge = (track->Charge()>0) ? 1 : -1;
2573 //Int_t trackCharge = 0;
2574 //trackCharge = sign(track->Charge());
2575
2576 //if (trackCharge>0)            return fPartHadronPlus;
2577 //else if (trackCharge<0)       return fPartHadronMinus;
2578
2579   //Printf(" --->>>  %f <--> %f <--> %f <--> %f <--> %f <--> %f",track->Pt(),fgITSPIDmomcutMin[2], fgITSPIDmomcutMax[2], fgITSPIDsigmacutMin[2], fgITSPIDsigmacutMax[2],fnsigmas[2][0]);
2580
2581   for(Int_t iSpecies=0; iSpecies < AliPID::kSPECIES; iSpecies++)
2582   {
2583     if ( (track->Pt() > fgITSPIDmomcutMin[iSpecies] && track->Pt() < fgITSPIDmomcutMax[iSpecies]) && (fnsigmas[iSpecies][0] > fgITSPIDsigmacutMin[iSpecies] && fnsigmas[iSpecies][0] < fgITSPIDsigmacutMax[iSpecies]) )
2584       isITSok[iSpecies] = kTRUE;
2585     if ( (track->Pt() > fgTPCPIDmomcutMin[iSpecies] && track->Pt() < fgTPCPIDmomcutMax[iSpecies]) && (fnsigmas[iSpecies][1] > fgTPCPIDsigmacutMin[iSpecies] && fnsigmas[iSpecies][1] < fgTPCPIDsigmacutMax[iSpecies]) )
2586       isTPCok[iSpecies] = kTRUE;
2587     if ( (track->Pt() > fgTOFPIDmomcutMin[iSpecies] && track->Pt() < fgTOFPIDmomcutMax[iSpecies]) && (fnsigmas[iSpecies][2] > fgTOFPIDsigmacutMin[iSpecies] && fnsigmas[iSpecies][2] < fgTOFPIDsigmacutMax[iSpecies]) )
2588       isTOFok[iSpecies] = kTRUE;
2589     if ( (track->Pt() > fgHMPIDPIDmomcutMin[iSpecies] && track->Pt() < fgHMPIDPIDmomcutMax[iSpecies]) && (fnsigmas[iSpecies][3] > fgHMPIDPIDsigmacutMin[iSpecies] && fnsigmas[iSpecies][3] < fgHMPIDPIDsigmacutMax[iSpecies]) )
2590       isHMPIDok[iSpecies] = kTRUE;
2591
2592     switch(iSpecies)
2593     {
2594       case 2:
2595         if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()>0 )
2596           return fPartPionPlus;
2597         else if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()<0 )
2598           return fPartPionMinus;
2599       break;
2600       case 3:
2601         if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()>0 )
2602           return fPartKaonPlus;
2603         else if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()<0 )
2604           return fPartKaonMinus;
2605       break;
2606       case 4:
2607         if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()>0 )
2608           return fPartProtonPlus;
2609         else if ( (isITSok[iSpecies] == kTRUE || isTPCok[iSpecies] == kTRUE || isTOFok[iSpecies] == kTRUE || isHMPIDok[iSpecies] == kTRUE) && track->Charge()<0 )
2610           return fPartProtonMinus;
2611       break;
2612     }
2613   }//iSpecies
2614
2615   return fPartUndefined;
2616 }
2617
2618 //___________________________________________________________
2619 Int_t AliAnalysisTaskPidPidCorrelations::GetParticleID(AliVTrack* trk, Int_t centbin, Bool_t fillQA)
2620 {
2621   if (!trk) { AliInfo(" ==== zero track pointer."); return fPartUndefined; }
2622
2623   Bool_t* pidFlag; pidFlag = new Bool_t[AliPID::kSPECIES];
2624   CalculateNSigmas((AliAODTrack*)trk, centbin, pidFlag, fillQA);
2625   delete[] pidFlag;
2626   Int_t mypid = -1;
2627   mypid = FindNSigma((AliAODTrack*)trk);
2628 //     Printf(" >>>>>>>>>>>>>>>>> mypid = %d",mypid);
2629   return mypid;
2630 }
2631
2632 //___________________________________________________________
2633 Int_t AliAnalysisTaskPidPidCorrelations::GetParticleIDMC(AliVTrack* trk, Int_t centbin, Bool_t fillQA)
2634 {
2635   if (!trk) { AliInfo(" ==== zero track pointer."); return fPartUndefined; }
2636
2637   if (fUseMC)//______ MC
2638   {
2639     if(!fMyMcArray) AliFatal("Error: AliAnalysisTaskPidPidCorrelations::GetParticleID called on data\n");
2640
2641     Int_t pdgCode = 999;
2642 //     AliAODMCParticle* partMC = dynamic_cast<AliAODMCParticle*>(fMyMcArray->At(TMath::Abs(trk->GetLabel())));
2643     AliAODMCParticle* partMC = (AliAODMCParticle*)fMyMcArray->At(TMath::Abs(trk->GetLabel()));
2644     if (!partMC){ AliError("Cannot get MC particle"); return fPartUndefined; }
2645     pdgCode = partMC->GetPdgCode();
2646 //     delete partMC; partMC=0x0;
2647
2648     //Int_t qPart = 0;
2649     //qPart  = (Int_t)partMC->GetPDG()->Charge();
2650     //if (qPart<0) return fPartHadronMinus;
2651     //else if (qPart<0) return fPartHadronPlus;
2652
2653     switch(pdgCode)
2654     {
2655       case 2212:
2656         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartProtonQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartProtonQA][ipid]);}}
2657           return fPartProtonPlus;
2658       break;
2659       case -2212:
2660         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartProtonQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartProtonQA][ipid]);}}
2661           return fPartProtonPlus;
2662       break;
2663       case 321:
2664         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartKaonQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartKaonQA][ipid]);}}
2665           return fPartKaonPlus;
2666       break;
2667       case -321:
2668         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartKaonQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartKaonQA][ipid]);}}
2669           return fPartKaonMinus;
2670       break;
2671       case 211:
2672         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartPionQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartPionQA][ipid]);}}
2673           return fPartPionPlus;
2674       break;
2675       case -211:
2676         if (fillQA) {for(Int_t ipid=0;ipid<kMyNSigmaPIDType;ipid++){TH2F* h = GetHisto2D(Form("NSigmaMC_Cent%02d_%d_%d",centbin,fPartPionQA,ipid));h -> Fill(trk->Pt(),fnsigmas[fPartPionQA][ipid]);}}
2677           return fPartPionMinus;
2678       break;
2679       default:
2680         return fPartUndefined;
2681     }
2682   } // MC
2683   return fPartUndefined;
2684 }
2685
2686 //___________________________________________________________
2687 Double_t AliAnalysisTaskPidPidCorrelations::MakeTPCPID(AliAODTrack* track, Double_t* nSigma)
2688 {
2689         //___ fills nSigma array with TPC nsigmas for e, mu, pi, K, p
2690
2691         /* check TPC PID */
2692         if (!HasTPCPID(track)) return -1.;
2693
2694         /* get TPC info */
2695         Double_t ptpc   = track -> GetTPCmomentum();
2696         Double_t dEdx   = track -> GetTPCsignal();
2697         Int_t dEdxN     = track -> GetTPCsignalN();
2698
2699         /* loop over particles */
2700         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
2701         {
2702                 Double_t bethe = fPIDResponse -> GetTPCResponse().GetExpectedSignal(ptpc, (AliPID::EParticleType)iSpecies);
2703                 Double_t diff = dEdx - bethe;
2704                 Double_t sigma = fPIDResponse -> GetTPCResponse().GetExpectedSigma(ptpc, dEdxN, (AliPID::EParticleType)iSpecies);
2705                 nSigma[iSpecies] = diff/sigma;
2706         }
2707         return dEdx;
2708 }
2709
2710 //___________________________________________________________
2711 Double_t AliAnalysisTaskPidPidCorrelations::MakeTOFPID(AliAODTrack* track, Double_t* nSigma)
2712 {
2713         //___ fills nSigma array with TOF nsigmas for e, mu, pi, K, p
2714
2715         /* check TOF PID */
2716         if (!HasTOFPID(track)) return -1.;
2717
2718         /* get TOF info */
2719         Double_t p      = track -> P();
2720         Double_t time   = track -> GetTOFsignal() - fPIDResponse -> GetTOFResponse().GetStartTime(p);
2721         Double_t timei[5];
2722         track -> GetIntegratedTimes(timei,AliPID::kSPECIES);
2723
2724         /* loop over particles */
2725         for (Int_t iSpecies = 0; iSpecies < AliPID::kSPECIES; iSpecies++)
2726         {
2727                 Double_t timez = time - timei[iSpecies];
2728                 Double_t sigma = fPIDResponse -> GetTOFResponse().GetExpectedSigma(p, timei[iSpecies], AliPID::ParticleMass(iSpecies));
2729                 nSigma[iSpecies] = timez/sigma;
2730         }
2731         //return timei[0]/time;
2732         return GetBeta(track);
2733 }
2734
2735 //___________________________________________________________
2736 Bool_t AliAnalysisTaskPidPidCorrelations::HasTPCPID(AliAODTrack* track) const
2737 {
2738   // check PID signal
2739   AliPIDResponse::EDetPidStatus statusTPC = fPIDResponse -> CheckPIDStatus(AliPIDResponse::kTPC,track);
2740   if (statusTPC != AliPIDResponse::kDetPidOk)
2741     return kFALSE;
2742   return kTRUE;
2743 }
2744
2745 //___________________________________________________________
2746 Bool_t AliAnalysisTaskPidPidCorrelations::HasTOFPID(AliAODTrack* track) const
2747 {
2748   AliPIDResponse::EDetPidStatus statusTOF = fPIDResponse -> CheckPIDStatus(AliPIDResponse::kTOF,track);
2749   if(statusTOF != AliPIDResponse::kDetPidOk)
2750     return kFALSE;
2751   if(0)
2752   {
2753     Int_t startTimeMask = fPIDResponse -> GetTOFResponse().GetStartTimeMask(track->P());
2754     if (startTimeMask < 0) return kFALSE;
2755   }
2756   return kTRUE;
2757 }
2758
2759 //___________________________________________________________
2760 Double_t AliAnalysisTaskPidPidCorrelations ::GetBeta(AliAODTrack* track)
2761 {
2762         //it is called only when TOF PID is available
2763         Double_t p      = track -> P();
2764         Double_t time   = track -> GetTOFsignal() - fPIDResponse -> GetTOFResponse().GetStartTime(p);
2765         Double_t timei[5];
2766         track -> GetIntegratedTimes(timei,AliPID::kSPECIES);
2767         return timei[0]/time;
2768 }
2769
2770 //___________________________________________________________
2771 void AliAnalysisTaskPidPidCorrelations::RemoveDuplicates(TObjArray* tracks)
2772 {
2773         // remove particles with the same label
2774         Int_t before = tracks -> GetEntriesFast();
2775
2776         for (Int_t i=0; i<before; ++i)
2777         {
2778                 AliVParticle* part = (AliVParticle*) tracks -> At(i);
2779
2780                 for (Int_t j=i+1; j<before; ++j)
2781                 {
2782         AliVParticle* part2 = (AliVParticle*) tracks -> At(j);
2783
2784         if (part -> GetLabel() == part2 -> GetLabel())
2785         {
2786                         Printf("Removing %d with label %d (duplicated in %d)", i, part -> GetLabel(), j); part -> Dump(); part2 -> Dump();
2787                         TObject* object = tracks -> RemoveAt(i);
2788                         if (tracks -> IsOwner())
2789                         delete object;
2790                 break;
2791         }
2792                 }
2793         }
2794         tracks -> Compress();
2795
2796         if (before > tracks -> GetEntriesFast())
2797                 AliInfo(Form("Reduced Form %d to %d", before, tracks -> GetEntriesFast()));
2798 }
2799
2800 //___________________________________________________________
2801 void AliAnalysisTaskPidPidCorrelations::RemoveWeakDecays(TObjArray* tracks, TObject* mcObj)
2802 {
2803   // remove particles Form weak decays
2804   // <tracks> can be the following cases:
2805   // a. tracks: in this case the label is taken and then case b.
2806   // b. particles: it is checked if IsSecondaryFromWeakDecay is true
2807   // <mcObj> can be AOD (TClonesArray) or ESD (AliMCEvent)
2808
2809   TClonesArray* arrayMC = 0;
2810   AliMCEvent* mcEvent = 0;
2811   if (mcObj -> InheritsFrom("AliMCEvent"))
2812     mcEvent = static_cast<AliMCEvent*>(mcObj);
2813   else if (mcObj -> InheritsFrom("TClonesArray"))
2814     arrayMC = static_cast<TClonesArray*>(mcObj);
2815   else
2816   {
2817     mcObj -> Dump();
2818     AliFatal("Invalid object passed");
2819   }
2820
2821   Int_t before = tracks->GetEntriesFast();
2822
2823   for (Int_t i=0; i<before; ++i)
2824   {
2825     AliVParticle* part = (AliVParticle*) tracks->At(i);
2826
2827     if (part->InheritsFrom("AliESDtrack") || part->InheritsFrom("AliAODTrack"))
2828       part = ((mcEvent) ? mcEvent->GetTrack(TMath::Abs(part->GetLabel())) : (AliVParticle*)arrayMC->At(TMath::Abs(part->GetLabel())));
2829
2830     if (part->InheritsFrom("AliAODMCParticle"))
2831     {
2832       if (!((AliAODMCParticle*) part)->IsSecondaryFromWeakDecay())
2833                 continue;
2834     }
2835     else if (part->InheritsFrom("AliMCParticle") && mcEvent)
2836     {
2837       if (!(mcEvent->Stack()->IsSecondaryFromWeakDecay(part->GetLabel())))
2838                 continue;
2839     }
2840     else
2841     {
2842       part -> Dump();
2843       AliFatal("Unknown particle");
2844     }
2845
2846 //     Printf("Removing %d with label %d", i, part->GetLabel()); part->Dump();
2847     TObject* object = tracks -> RemoveAt(i);
2848     if (tracks -> IsOwner())
2849       delete object;
2850   }
2851
2852   tracks -> Compress();
2853
2854   if (before > tracks->GetEntriesFast())
2855     AliInfo(Form("Reduced Form %d to %d", before, tracks->GetEntriesFast()));
2856 }
2857
2858 //________________________________________________________________________
2859 Double_t AliAnalysisTaskPidPidCorrelations::DeltaPhi(Double_t Dphi) const
2860 {
2861   if (Dphi < -0.5*TMath::Pi())  Dphi += TMath::TwoPi();
2862   if (Dphi > 1.5*TMath::Pi())   Dphi -= TMath::TwoPi();
2863   return Dphi;
2864 }
2865
2866 //___________________________________________________________
2867 TH2F* AliAnalysisTaskPidPidCorrelations::GetHisto2D(const Char_t* name)
2868 {
2869   return (TH2F*) fList -> FindObject(name);
2870 }
2871
2872 //___________________________________________________________
2873 Int_t AliAnalysisTaskPidPidCorrelations::GetCentBin(Double_t percentile)
2874 {
2875 //   Double_t percentile = fMyAODEvent -> GetCentrality() -> GetCentralityPercentile(fCentralityEstimator.Data());
2876   Int_t bin = fCentAxis -> FindBin(percentile) - 1;
2877   if (bin >= fNbinsCent) bin = -1;
2878   return bin;
2879 }
2880
2881 //___________________________________________________________
2882 void AliAnalysisTaskPidPidCorrelations::SetCentBinning(Int_t nBins, Double_t* limits)
2883 {
2884   if (nBins > fNMaxBinsCentrality)
2885   {
2886     AliInfo(Form("WARNING : only %d centrality bins (out of the %d proposed) will be considered",fNMaxBinsCentrality,nBins));
2887     nBins = fNMaxBinsCentrality;
2888   }
2889   if (nBins <= 0)
2890   {
2891     AliInfo("WARNING : at least one centrality bin must be considered");
2892     nBins = 1;
2893   }
2894   fNbinsCent = nBins;
2895   fCentAxis  = new TAxis(fNbinsCent, limits);
2896 }
2897
2898 //___________________________________________________________
2899 Int_t AliAnalysisTaskPidPidCorrelations::GetZvtxBin(Double_t zvtx)
2900 {
2901   Int_t bin = fZvtxAxis -> FindBin(zvtx) - 1;
2902   if (bin >= fNbinsZvtx) bin = -1;
2903   return bin;
2904 }
2905
2906 //___________________________________________________________
2907 void AliAnalysisTaskPidPidCorrelations::SetZvtxBinning(Int_t nBins, Double_t* limits)
2908 {
2909   if (nBins > fNMaxBinsZvtx)
2910   {
2911     AliInfo(Form("WARNING : only %d z-vertex bins (out of the %d proposed) will be considered",fNMaxBinsZvtx,nBins));
2912     nBins = fNMaxBinsZvtx;
2913   }
2914   if (nBins <= 0)
2915   {
2916     AliInfo("WARNING : at least one centrality bin must be considered");
2917     nBins = 1;
2918   }
2919   fNbinsZvtx = nBins;
2920   fZvtxAxis  = new TAxis(fNbinsZvtx, limits);
2921 }
2922
2923 //___________________________________________________________
2924 Int_t AliAnalysisTaskPidPidCorrelations::GetPtBin(Double_t valPt)
2925 {
2926   Int_t bin = fPtAxis -> FindBin(valPt) - 1;
2927   if (bin >= fNbinsPt) bin = -1;
2928   return bin;
2929 }
2930
2931 //__________________________________________________________________________________________________
2932 void AliAnalysisTaskPidPidCorrelations::SetPtBinning(Int_t nBins, Double_t* limits)
2933 {
2934   if (nBins > fNMaxBinsPt)
2935   {
2936     AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsPt,nBins));
2937     nBins = fNMaxBinsPt;
2938   }
2939   if (nBins <= 0)
2940   {
2941     AliInfo("WARNING : at least one pt bin must be considered");
2942     nBins = 1;
2943   }
2944   fNbinsPt = nBins;
2945   fPtAxis  = new TAxis(fNbinsPt, limits);
2946 }
2947
2948 //___________________________________________________________
2949 Int_t AliAnalysisTaskPidPidCorrelations::GetEtaBin(Double_t valEta)
2950 {
2951   Int_t bin = fEtaAxis -> FindBin(valEta) - 1;
2952   if (bin >= fNbinsEta) bin = -1;
2953   return bin;
2954 }
2955
2956 //__________________________________________________________________________________________________
2957 void AliAnalysisTaskPidPidCorrelations::SetEtaBinning(Int_t nBins, Double_t* limits)
2958 {
2959   if (nBins > fNMaxBinsEta)
2960   {
2961     AliInfo(Form("WARNING : only %d pt bins (out of the %d proposed) will be considered",fNMaxBinsEta,nBins));
2962     nBins = fNMaxBinsEta;
2963   }
2964   if (nBins <= 0)
2965   {
2966     AliInfo("WARNING : at least one pt bin must be considered");
2967     nBins = 1;
2968   }
2969   fNbinsEta = nBins;
2970   fEtaAxis  = new TAxis(nBins, limits);
2971 }
2972
2973 //________________________________________________________________________
2974 Double_t AliAnalysisTaskPidPidCorrelations::GetDPhiStar(Double_t phi1, Double_t pt1, Double_t charge1, Double_t phi2, Double_t pt2, Double_t charge2, Double_t radius, Double_t bSign)
2975 {
2976   //
2977   // calculates dphistar
2978   //
2979
2980   Double_t dphistar = phi1 - phi2 - charge1 * bSign * TMath::ASin(0.075 * radius / pt1) + charge2 * bSign * TMath::ASin(0.075 * radius / pt2);
2981
2982   static const Double_t kPi = TMath::Pi();
2983
2984   // circularity
2985 //   if (dphistar > 2 * kPi)
2986 //     dphistar -= 2 * kPi;
2987 //   if (dphistar < -2 * kPi)
2988 //     dphistar += 2 * kPi;
2989
2990   if (dphistar > kPi)
2991     dphistar = kPi * 2 - dphistar;
2992   if (dphistar < -kPi)
2993     dphistar = -kPi * 2 - dphistar;
2994   if (dphistar > kPi) // might look funny but is needed
2995     dphistar = kPi * 2 - dphistar;
2996
2997   return dphistar;
2998 }
2999
3000 //________________________________________________________________________ copy&paste from $ALICE_ROOT/PWGCF/Correlations/Base/AliUEHistograms.h
3001 Float_t AliAnalysisTaskPidPidCorrelations::GetInvMassSquared(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
3002 {
3003   // calculate inv mass squared
3004   // same can be achieved, but with more computing time with
3005   /*TLorentzVector photon, p1, p2;
3006   p1.SetPtEtaPhiM(triggerParticle->Pt(), triggerEta, triggerParticle->Phi(), 0.510e-3);
3007   p2.SetPtEtaPhiM(particle->Pt(), eta[j], particle->Phi(), 0.510e-3);
3008   photon = p1+p2;
3009   photon.M()*/
3010
3011   Float_t tantheta1 = 1e10;
3012
3013   if (eta1 < -1e-10 || eta1 > 1e-10)
3014   {
3015     Float_t expTmp = TMath::Exp(-eta1);
3016     tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3017   }
3018
3019   Float_t tantheta2 = 1e10;
3020   if (eta2 < -1e-10 || eta2 > 1e-10)
3021   {
3022     Float_t expTmp = TMath::Exp(-eta2);
3023     tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3024   }
3025
3026   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3027   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3028
3029   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( TMath::Cos(phi1 - phi2) + 1.0 / tantheta1 / tantheta2 ) ) );
3030
3031 //   Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3032
3033   return mass2;
3034 }
3035
3036 //________________________________________________________________________ copy&paste from $ALICE_ROOT/PWGCF/Correlations/Base/AliUEHistograms.h
3037 Float_t AliAnalysisTaskPidPidCorrelations::GetInvMassSquaredCheap(Float_t pt1, Float_t eta1, Float_t phi1, Float_t pt2, Float_t eta2, Float_t phi2, Float_t m0_1, Float_t m0_2)
3038 {
3039   // calculate inv mass squared approximately
3040
3041   Float_t tantheta1 = 1e10;
3042
3043   if (eta1 < -1e-10 || eta1 > 1e-10)
3044   {
3045     Float_t expTmp = 1.0-eta1+eta1*eta1/2-eta1*eta1*eta1/6+eta1*eta1*eta1*eta1/24;
3046     tantheta1 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3047   }
3048
3049   Float_t tantheta2 = 1e10;
3050   if (eta2 < -1e-10 || eta2 > 1e-10)
3051   {
3052     Float_t expTmp = 1.0-eta2+eta2*eta2/2-eta2*eta2*eta2/6+eta2*eta2*eta2*eta2/24;
3053     tantheta2 = 2.0 * expTmp / ( 1.0 - expTmp*expTmp);
3054   }
3055
3056   Float_t e1squ = m0_1 * m0_1 + pt1 * pt1 * (1.0 + 1.0 / tantheta1 / tantheta1);
3057   Float_t e2squ = m0_2 * m0_2 + pt2 * pt2 * (1.0 + 1.0 / tantheta2 / tantheta2);
3058
3059   // fold onto 0...pi
3060   Float_t deltaPhi = TMath::Abs(phi1 - phi2);
3061   while (deltaPhi > TMath::TwoPi())
3062     deltaPhi -= TMath::TwoPi();
3063   if (deltaPhi > TMath::Pi())
3064     deltaPhi = TMath::TwoPi() - deltaPhi;
3065
3066   Float_t cosDeltaPhi = 0;
3067   if (deltaPhi < TMath::Pi()/3)
3068     cosDeltaPhi = 1.0 - deltaPhi*deltaPhi/2 + deltaPhi*deltaPhi*deltaPhi*deltaPhi/24;
3069   else if (deltaPhi < 2*TMath::Pi()/3)
3070     cosDeltaPhi = -(deltaPhi - TMath::Pi()/2) + 1.0/6 * TMath::Power((deltaPhi - TMath::Pi()/2), 3);
3071   else
3072     cosDeltaPhi = -1.0 + 1.0/2.0*(deltaPhi - TMath::Pi())*(deltaPhi - TMath::Pi()) - 1.0/24.0 * TMath::Power(deltaPhi - TMath::Pi(), 4);
3073
3074   Float_t mass2 = m0_1 * m0_1 + m0_2 * m0_2 + 2 * ( TMath::Sqrt(e1squ * e2squ) - ( pt1 * pt2 * ( cosDeltaPhi + 1.0 / tantheta1 / tantheta2 ) ) );
3075
3076 //   Printf(Form("%f %f %f %f %f %f %f %f %f", pt1, eta1, phi1, pt2, eta2, phi2, m0_1, m0_2, mass2));
3077
3078   return mass2;
3079 }
3080
3081 //___________________________________________________________
3082 void AliAnalysisTaskPidPidCorrelations::PrintPoolManagerContents()
3083 {
3084   // Determine the number of pools in the pool manager.
3085   AliEventPool* poolin = fPoolMgr -> GetEventPool(0,0);
3086   Int_t NPoolsCentrality = 0;
3087   while (poolin)
3088   {
3089     NPoolsCentrality++;
3090     poolin = fPoolMgr -> GetEventPool(NPoolsCentrality,0);
3091   }
3092
3093   poolin = fPoolMgr -> GetEventPool(0,0);
3094   Int_t NPoolsVtxZ = 0;
3095   while (poolin)
3096   {
3097     NPoolsVtxZ++;
3098     poolin = fPoolMgr -> GetEventPool(0,NPoolsVtxZ);
3099   }
3100
3101   // Loop over all Pools in the matrix of the pool manager.
3102   std::cout<<" Pool manager contents: (Nevt,NTrack)"<<std::endl;
3103   for (Int_t iCentrality = 0; iCentrality < NPoolsCentrality; iCentrality++)
3104   {
3105     std::cout<<Form("Centrality Bin: %2i --> ", iCentrality);
3106
3107     for (Int_t iVtxZ = 0; iVtxZ < NPoolsVtxZ; iVtxZ++)
3108     {
3109       poolin = fPoolMgr->GetEventPool(iCentrality, iVtxZ);
3110       std::cout<<Form("(%2i,%4i) ",poolin->GetCurrentNEvents(), poolin->NTracksInPool());
3111
3112 //       Double_t vars[] = { iCentrality, iVtxZ, poolin->GetCurrentNEvents(), poolin->NTracksInPool() };
3113 //       fHistPoolMgrQA -> Fill(vars);
3114     }
3115     std::cout<<std::endl;
3116   }
3117 }
3118
3119 //___________________________________________________________
3120 void AliAnalysisTaskPidPidCorrelations::SetupForMixing()
3121 {
3122   const Int_t trackDepth        = fMinNumTrack;
3123   const Int_t poolsize          = fPoolSize;
3124   Double_t centralityBins[] = { fCentralityPercentileMin, fCentralityPercentileMax };
3125   const Int_t nCentralityBins = length(centralityBins)-1;
3126   fPoolMgr = new AliEventPoolManager(poolsize, trackDepth, nCentralityBins/*fCentAxis->GetNbins()+1*/, centralityBins/*const_cast<Double_t*>(fCentAxis->GetXbins()->GetArray())*/, fZvtxAxis->GetNbins()+1, const_cast<Double_t*>(fZvtxAxis->GetXbins()->GetArray()));
3127 //   fPoolMgr -> SetTargetValues(fMinNumTrack, 0.1, 5);
3128 }
3129
3130 //___________________________________________________________
3131 TObjArray* AliAnalysisTaskPidPidCorrelations::CloneAndReduceTrackList(TObjArray* tracks)
3132 {
3133   // clones a track list by using AliPidPidCorrelationReducedTrack which uses much less memory (used for event mixing)
3134
3135   TObjArray* tracksClone = new TObjArray;
3136   tracksClone -> SetOwner(kTRUE);
3137
3138   for (Int_t i=0; i < tracks->GetEntriesFast(); i++)
3139   {
3140     AliPidPidCorrelationReducedTrack* particle = (AliPidPidCorrelationReducedTrack*) tracks->UncheckedAt(i);
3141     AliPidPidCorrelationReducedTrack* copy = new AliPidPidCorrelationReducedTrack(particle->GetMyPartID(),particle->Eta(),particle->Phi(),particle->Pt(),particle->Charge());
3142     copy -> SetUniqueID(particle->GetUniqueID());
3143     tracksClone -> Add(copy);
3144   }
3145   return tracksClone;
3146 }
3147 //___________________________________________________________
3148 void AliAnalysisTaskPidPidCorrelations::FillCFcontainers(TObjArray* mc, TObjArray* reco, TObjArray* recomatch, Double_t cent/*, Double_t zvtx*//*const AliVEvent* event*/)
3149 {
3150   if (!fMyCFCont) return;
3151 /*
3152   Double_t zvtx=0., zvtxMC=0.;
3153   if ( event->IsA() == AliAODEvent::Class() )
3154   {
3155     const AliVVertex* vertex = (AliVVertex*)static_cast<const AliAODEvent*>(event)->GetPrimaryVertexSPD();
3156     zvtx = vertex->GetZ();
3157   }
3158   if ( event->IsA() == AliAODEvent::Class() )
3159   {
3160     AliAODMCHeader* aodMCHeader = static_cast<AliAODMCHeader*> (static_cast<const AliAODEvent*>(event)->FindListObject(AliAODMCHeader::StdBranchName()));
3161     zvtxMC = aodMCHeader->GetVtxZ();
3162   }
3163   */
3164   for (Int_t istep=0; istep < kNsteps; istep++)
3165   {
3166     TObjArray* list = mc;
3167     if (istep == 1)
3168       list = reco;
3169     else if (istep == 2)
3170       list = recomatch;
3171
3172     if (!list)
3173       continue;
3174
3175     Double_t cfValue[kNvars];
3176     for(Int_t ie = 0 ; ie < list->GetEntriesFast(); ie++)
3177     {
3178       AliPidPidCorrelationReducedTrack* mytrack = dynamic_cast<AliPidPidCorrelationReducedTrack*>(list->At(ie));
3179
3180       cfValue[kVarEta]  = mytrack->Eta();
3181 //       cfValue[kVarPhi]       = mytrack->Phi();
3182       cfValue[kVarPt]   = mytrack->Pt();
3183       cfValue[kVarPID]  = mytrack->GetMyPartID();
3184       cfValue[kVarCent] = cent;
3185       cfValue[kVarZvtx] = ( istep == kStepRec ) ? (fMyAODEvent->GetPrimaryVertex()->GetZ()) : (fMyMcHeader->GetVtxZ()); //! FIXME
3186 //       cfValue[kVarZvtx]      = zvtx;
3187
3188       fMyCFCont->Fill(cfValue,istep);
3189     }
3190   }
3191 }
3192
3193 //___________________________________________________________
3194 TString AliAnalysisTaskPidPidCorrelations::GetOutputListName() const
3195 {
3196   TString listName("listPIDCorr");
3197   listName += TString::Format("_%smix",         fFillMixed ? "" : "no");
3198   listName += TString::Format("_cent%.0f%.0f", fCentralityPercentileMin, fCentralityPercentileMax);
3199   listName += TString::Format("_ptMin%.0fMeV",  1e3*fTrackPtMin);
3200   return listName;
3201 }
3202
3203 //___________________________________________________________
3204 void AliAnalysisTaskPidPidCorrelations::Terminate(Option_t *)
3205 {
3206   //AliInfo("Terminate","");
3207   AliAnalysisTaskSE::Terminate();
3208
3209   fList   = dynamic_cast<TList*>(GetOutputData(1));
3210   if(!fList) { Printf("ERROR: could not retrieve TList fList"); return; }
3211
3212   fMyCFCont = dynamic_cast<AliCFContainer*> (GetOutputData(2));
3213   if ( !fMyCFCont ) { AliError("Cannot find container in file"); return; }
3214
3215   if(fPoolMgr) delete fPoolMgr;
3216
3217   return;
3218 }
3219