PWGGA/PHOSTasks/PHOS_Correlation/AliPHOSCorrelations.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2014, 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 // Analysis task for identifion PHOS cluster from Pi0 and extracting pi0-hadron correlation.
  17 // Author:     Daniil Ponomarenko <Daniil.Ponomarenko@cern.ch>
  18 // 20-Sept-2014
  19
  20 #include <Riostream.h>
  21 #include "THashList.h"
  22 #include "TObjArray.h"
  23 #include "TClonesArray.h"
  24 #include "TH1F.h"
  25 #include "TH2F.h"
  26 #include "TH2I.h"
  27 #include "TH3F.h"
  28 #include "TH3D.h"
  29 #include "TMath.h"
  30 #include "TVector3.h"
  31 #include "TProfile.h"
  32
  33 #include "AliAnalysisManager.h"
  34 #include "AliAnalysisTaskSE.h"
  35 #include "AliPHOSCorrelations.h"
  36 #include "AliPHOSGeometry.h"
  37 #include "AliESDEvent.h"
  38 #include "AliESDCaloCluster.h"
  39 #include "AliESDVertex.h"
  40 #include "AliESDtrackCuts.h"
  41 #include "AliESDtrack.h"
  42 #include "AliAODTrack.h"
  43 #include "AliVTrack.h"
  44 #include "AliPID.h"
  45 #include "AliTriggerAnalysis.h"
  46 #include "AliPIDResponse.h"
  47 #include "AliPHOSEsdCluster.h"
  48 #include "AliCDBManager.h"
  49 #include "AliPHOSCalibData.h"
  50 #include "AliCentrality.h"
  51 #include "AliEventplane.h"
  52 #include "AliOADBContainer.h"
  53 #include "AliAODEvent.h"
  54 #include "AliAODCaloCells.h"
  55 #include "AliAODCaloCluster.h"
  56 #include "AliCaloPhoton.h"
  57 #include "AliAODVertex.h"
  58 #include "AliInputEventHandler.h"
  59
  60 using std::cout;
  61 using std::endl;
  62
  63 ClassImp(AliPHOSCorrelations)
  64
  65 //_______________________________________________________________________________
  66 AliPHOSCorrelations::AliPHOSCorrelations()
  67 :AliAnalysisTaskSE(),
  68     fPHOSGeo(0x0),
  69     fOutputContainer(0x0),
  70     fEvent(0x0),
  71     fEventESD(0x0),
  72     fEventAOD(0x0),
  73     fEventHandler(0),
  74     fCaloPhotonsPHOS(0x0),
  75     fTracksTPC(0x0),
  76     fCaloPhotonsPHOSLists(0x0),
  77     fTracksTPCLists(0x0),
  78     fRunNumber(-999),
  79     fInternalRunNumber(0),
  80     fPeriod("0x0"),
  81     fPHOSEvent(false),
  82     fMBEvent(false),
  83     fNVtxZBins(10),
  84     fVtxEdges(10),
  85     fCentEdges(10),
  86     fCentNMixed(),
  87     fNEMRPBins(6),
  88     fAssocBins(),
  89     fVertexVector(),
  90     fVtxBin(0),
  91     fCentralityEstimator("V0M"),
  92     fCentrality(0.),
  93     fCentBin(0),
  94     fHaveTPCRP(0),
  95     fRP(0.),
  96     fEMRPBin(0),
  97     fMaxAbsVertexZ(10.),
  98     fCentralityLowLimit(0.),
  99     fCentralityHightLimit(90),
 100     fMinClusterEnergy(0.3),
 101     fMinBCDistance(0),
 102     fMinNCells(3),
 103     fMinM02(0.2),
 104     fTOFCutEnabled(1),
 105     fTOFCut(100.e-9),
 106     fUseMassWindowParametrisation(true),
 107     fMassInvMeanMin(0.13),
 108     fMassInvMeanMax(0.145),
 109     fNSigmaWidth(0.),
 110     fUseEfficiency(true),
 111     fESDtrackCuts(0x0),
 112     fSelectHybridTracks(0),
 113     fTrackStatus(0),
 114     fTrackFilterMask(786),
 115     fSelectSPDHitTracks(0),
 116     fSelectFractionTPCSharedClusters(kTRUE),
 117     fCutTPCSharedClustersFraction(0.4)
 118 {
 119      // Constructor
 120
 121     fMassMean[0]  = 1.00796e-05 ;
 122     fMassMean[1]  = 0.136096    ;
 123     fMassSigma[0] = 0.00383029 ;
 124     fMassSigma[1] = 0.0041709 ;
 125     fMassSigma[2] = 0.00468736 ;
 126
 127     const Int_t nPtAssoc = 10 ;
 128     Double_t ptAssocBins[nPtAssoc] = {0.,0.5,1.0,1.5,2.0,3.,5.,7.,10.,16} ;
 129     fAssocBins.Set(nPtAssoc,ptAssocBins) ;
 130
 131     const int nbins = 9;
 132     Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
 133     TArrayD centEdges( nbins+1, edges );
 134     Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
 135     TArrayI centNMixed(nbins, nMixed);
 136     SetCentralityBinning(centEdges, centNMixed);
 137     SetVertexBinning();
 138
 139     fVertex[0] = 0;
 140     fVertex[1] = 0;
 141     fVertex[2] = 0;
 142
 143     //SetGeometry();
 144
 145     ZeroingVariables();
 146 }
 147
 148 //_______________________________________________________________________________
 149 AliPHOSCorrelations::AliPHOSCorrelations(const char *name)
 150 :AliAnalysisTaskSE(name),
 151     fPHOSGeo(0x0),
 152     fOutputContainer(0x0),
 153     fEvent(0x0),
 154     fEventESD(0x0),
 155     fEventAOD(0x0),
 156     fEventHandler(0),
 157     fCaloPhotonsPHOS(0x0),
 158     fTracksTPC(0x0),
 159     fCaloPhotonsPHOSLists(0x0),
 160     fTracksTPCLists(0x0),
 161     fRunNumber(-999),
 162     fInternalRunNumber(0),
 163     fPeriod("0x0"),
 164     fPHOSEvent(false),
 165     fMBEvent(false),
 166     fNVtxZBins(10),
 167     fVtxEdges(10),
 168     fCentEdges(10),
 169     fCentNMixed(),
 170     fNEMRPBins(6),
 171     fAssocBins(),
 172     fVertexVector(),
 173     fVtxBin(0),
 174     fCentralityEstimator("V0M"),
 175     fCentrality(0.),
 176     fCentBin(0),
 177     fHaveTPCRP(0),
 178     fRP(0.),
 179     fEMRPBin(0),
 180     fMaxAbsVertexZ(10.),
 181     fCentralityLowLimit(0.),
 182     fCentralityHightLimit(90),
 183     fMinClusterEnergy(0.3),
 184     fMinBCDistance(0),
 185     fMinNCells(3),
 186     fMinM02(0.2),
 187     fTOFCutEnabled(1),
 188     fTOFCut(100.e-9),
 189     fUseMassWindowParametrisation(true),
 190     fMassInvMeanMin(0.13),
 191     fMassInvMeanMax(0.145),
 192     fNSigmaWidth(0.),
 193     fUseEfficiency(true),
 194     fESDtrackCuts(0x0),
 195     fSelectHybridTracks(0),
 196     fTrackStatus(0),
 197     fTrackFilterMask(786),
 198     fSelectSPDHitTracks(0),
 199     fSelectFractionTPCSharedClusters(kTRUE),
 200     fCutTPCSharedClustersFraction(0.4)
 201 {
 202     // Constructor
 203     // Output slots #0 write into a TH1 container
 204     DefineOutput(1,THashList::Class());
 205
 206     fMassMean[0]  = 1.00796e-05 ;
 207     fMassMean[1]  = 0.136096    ;
 208     fMassSigma[0] = 0.00383029 ;
 209     fMassSigma[1] = 0.0041709 ;
 210     fMassSigma[2] = 0.00468736 ;
 211
 212     const Int_t nPtAssoc = 10 ;
 213     Double_t ptAssocBins[nPtAssoc] = {0.,0.5,1.0,1.5,2.0,3.,5.,7.,10.,16} ;
 214     fAssocBins.Set(nPtAssoc,ptAssocBins) ;
 215
 216     const int nbins = 9;
 217     Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
 218     TArrayD centEdges( nbins+1, edges );
 219     Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
 220     TArrayI centNMixed(nbins, nMixed);
 221     SetCentralityBinning(centEdges, centNMixed);
 222     SetVertexBinning();
 223
 224     fVertex[0] = 0;
 225     fVertex[1] = 0;
 226     fVertex[2] = 0;
 227
 228     //SetGeometry();
 229
 230     ZeroingVariables();
 231 }
 232
 233 //_______________________________________________________________________________
 234 AliPHOSCorrelations::~AliPHOSCorrelations()
 235 {
 236     if(fCaloPhotonsPHOS){
 237       delete fCaloPhotonsPHOS;
 238       fCaloPhotonsPHOS=0x0;
 239     }
 240
 241     if(fTracksTPC){
 242       delete fTracksTPC;
 243       fTracksTPC=0x0;
 244     }
 245
 246     if(fCaloPhotonsPHOSLists){
 247       fCaloPhotonsPHOSLists->SetOwner() ;
 248       delete fCaloPhotonsPHOSLists;
 249       fCaloPhotonsPHOSLists=0x0;
 250     }
 251
 252     if(fTracksTPCLists){
 253       fTracksTPCLists->SetOwner() ;
 254       delete fTracksTPCLists;
 255       fTracksTPCLists=0x0 ;
 256     }
 257
 258     if( fESDtrackCuts){
 259       delete fESDtrackCuts;
 260       fESDtrackCuts=0x0 ;
 261     }
 262
 263     if(fOutputContainer){
 264       delete fOutputContainer;
 265       fOutputContainer=0x0;
 266     }
 267 }
 268
 269 //_______________________________________________________________________________
 270 void AliPHOSCorrelations::UserCreateOutputObjects()
 271 {
 272     // Create histograms
 273     // Called once
 274     const Int_t     nRuns   = 200 ;
 275     const Int_t     ptMult  = 300 ;
 276     const Double_t  ptMin   = 0.  ;
 277     const Double_t  ptMax   = 30. ;
 278
 279     // Create histograms
 280     if(fOutputContainer != NULL) { delete fOutputContainer; }
 281     fOutputContainer = new THashList();
 282     fOutputContainer->SetOwner(kTRUE);
 283
 284     // Event selection
 285     fOutputContainer->Add(new TH1F( "hTriggerPassedEvents", "Event selection passed Cuts",                  60, 0., 60.  )) ;
 286
 287     fOutputContainer->Add(new TH1F( "hMBTriggerperRunNumber", "MB trigger vs run number",                   200, 0., 200.)) ;
 288     fOutputContainer->Add(new TH1F( "hCentralTriggerperRunNumber", "Central trigger vs run number",         200, 0., 200.)) ;
 289     fOutputContainer->Add(new TH1F( "hSemiCentralTriggerperRunNumber", "SemiCentral trigger vs run number", 200, 0., 200.)) ;
 290
 291     // Analysis event's progress
 292     fOutputContainer->Add(new TH1F( "hTotSelEvents", "Event selection", 15, 0., 15 )) ;
 293     fOutputContainer->Add(new TH2F( "hSelEvents", "Event selection", kTotalSelected+1, 0., double(kTotalSelected+1), nRuns,0., float(nRuns) )) ;
 294
 295     // Centrality, Reaction plane and Vertex selection
 296     fOutputContainer->Add(new TH2F( "hCentrality", "Event centrality of all events",                100, 0., 100., nRuns,0., float(nRuns)   )) ;
 297     fOutputContainer->Add(new TH2F( "hCentralityTriggerEvent", "Event centrality trigger events",   100, 0., 100., nRuns,0., float(nRuns)   )) ;
 298     fOutputContainer->Add(new TH2F( "hCentralityMBEvent", "Event centrality MB events",             100, 0., 100., nRuns,0., float(nRuns)   )) ;
 299     fOutputContainer->Add(new TH2F( "phiRPflat", "RP distribution with TPC flat",                   100, 0., 2.*TMath::Pi(), 20, 0., 100.   )) ;
 300
 301     fOutputContainer->Add(new TH1F( "hCentralityBining", "Event centrality bining of all events",            GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
 302     fOutputContainer->Add(new TH1F( "hCentralityBiningTrigger", "Event centrality bining of trigger events", GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
 303     fOutputContainer->Add(new TH1F( "hCentralityBiningMB", "Event centrality bining of MB events",           GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins() )) ;
 304
 305     fOutputContainer->Add(new TH1F( "phiRPflatBining", "Event RP bining of all events",                     GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
 306     fOutputContainer->Add(new TH1F( "phiRPflatBiningTrigger", "Event RP bining of trigger events",          GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
 307     fOutputContainer->Add(new TH1F( "phiRPflatBiningMB", "Event RP bining of MB events",                    GetNumberOfRPBins(), 0., GetNumberOfRPBins() )) ;
 308
 309     fOutputContainer->Add(new TH1F( "hVertexZBining", "Event vertex bining of all events",                  GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
 310     fOutputContainer->Add(new TH1F( "hVertexZBiningTrigger", "Event vertex bining of trigger events",       GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
 311     fOutputContainer->Add(new TH1F( "hVertexZBiningMB", "Event vertex bining of MB events",                 GetNumberOfVertexBins(), 0., GetNumberOfVertexBins() )) ;
 312
 313     // Mass selection
 314     fOutputContainer->Add(new TH1F( "massWindowPass", "Mass selection", 10,     0.,     10. )) ;
 315     fOutputContainer->Add(new TH2F( "massWindow", "mean & sigma",       100.,   0.085,  0.185,  50,    0., 0.05 )) ;
 316     // Cluster multiplisity
 317     fOutputContainer->Add(new TH2F( "hCluEvsClu", "ClusterMult vs E",   200,    0.,     10.,    100,    0., 100. )) ;
 318
 319     // Mixing progress
 320     // TODO: Fix filling maximum mixing depth value
 321     fOutputContainer->Add(new TH2F( "hCentralityMixingProgress", "Centrality Mixing Progress",  GetNumberOfCentralityBins(), 0., GetNumberOfCentralityBins(), 110, 0, 110 )) ;
 322     fOutputContainer->Add(new TH2F( "hVertexMixingProgress", "Vertex Mixing Progress",          GetNumberOfVertexBins(),     0., GetNumberOfVertexBins(),     110, 0, 110 )) ;
 323     fOutputContainer->Add(new TH2F( "hRPMixingProgress", "RP Mixing Progress",                  GetNumberOfRPBins(),         0., GetNumberOfRPBins(),         110, 0, 110 )) ;
 324
 325     // Time Of Flight
 326     fOutputContainer->Add(new TH1F( "hTOFcut", "TOF PHOS's clastrers distribution", 10000, -5000, 5000 )) ;
 327
 328     // Set hists, with track's and cluster's angle distributions.
 329     SetHistPtNumTrigger (ptMult, ptMin, ptMax);
 330     SetHistEtaPhi       (ptMult, ptMin, ptMax);
 331     SetHistPHOSClusterMap();
 332     SetHistMass         (ptMult, ptMin, ptMax);
 333     SetHistPtAssoc      (ptMult, ptMin, ptMax);
 334
 335     // Setup photon lists
 336     Int_t kapacity = fNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
 337     fCaloPhotonsPHOSLists = new TObjArray(kapacity);
 338     fCaloPhotonsPHOSLists->SetOwner();
 339
 340     fTracksTPCLists = new TObjArray(kapacity);
 341     fTracksTPCLists->SetOwner();
 342
 343     PostData(1, fOutputContainer);
 344 }
 345
 346 //_______________________________________________________________________________
 347 void AliPHOSCorrelations::SetHistPtNumTrigger(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
 348 {
 349     TString spid[4]={"all","cpv","disp","both"} ;
 350     for(Int_t ipid=0; ipid<4; ipid++)
 351     {
 352         fOutputContainer->Add(new TH1F(    Form("nTrigger_%s", spid[ipid].Data()),
 353                                         Form("Num of trigger particle %s", spid[ipid].Data()),
 354                                         ptMult+300, ptMin, ptMax ) );
 355         TH1F *h = static_cast<TH1F*>(fOutputContainer->Last()) ;
 356         h->Sumw2();
 357         h->GetXaxis()->SetTitle("Pt [GEV]");
 358     }
 359     for(Int_t ipid=0; ipid<4; ipid++)
 360     {
 361         fOutputContainer->Add(new TH1F( Form("nTrigger_%s_MB", spid[ipid].Data()),
 362                                         Form("Num of trigger particle %s", spid[ipid].Data()),
 363                                         ptMult+300, ptMin, ptMax ) );
 364         TH1F *h = static_cast<TH1F*>(fOutputContainer->Last()) ;
 365         h->Sumw2();
 366         h->GetXaxis()->SetTitle("Pt [GEV]");
 367     }
 368 }
 369
 370 //_______________________________________________________________________________
 371 void AliPHOSCorrelations::SetHistEtaPhi(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
 372 {
 373     // Set hists, with track's and cluster's angle distributions.
 374
 375     Float_t pi = TMath::Pi();
 376
 377     //===
 378     fOutputContainer->Add(new TH2F( "clu_phieta","Cluster's #phi & #eta distribution",
 379                             300, double(-1.8), double(-0.6),
 380                             300, double(-0.2), double(0.2) ) );
 381     TH2F * h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 382     h->GetXaxis()->SetTitle("#phi [rad]");
 383     h->GetYaxis()->SetTitle("#eta");
 384
 385     //===
 386     fOutputContainer->Add(new TH2F( "clusingle_phieta","Cluster's  #phi & #eta distribution",
 387                             300, double(-1.8), double(-0.6),
 388                             300, double(-0.2), double(0.2) ) );
 389     h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 390     h->GetXaxis()->SetTitle("#phi [rad]");
 391     h->GetYaxis()->SetTitle("#eta");
 392
 393     //===
 394     fOutputContainer->Add(new TH3F( "track_phieta","TPC track's  #phi & #eta distribution",
 395                             200, double(-pi-0.3), double(pi+0.3),
 396                             200, double(-0.9), double(0.9),
 397     ptMult, ptMin, ptMax ) );
 398     TH3F* h3 = static_cast<TH3F*>(fOutputContainer->FindObject("track_phieta")) ;
 399     h3->GetXaxis()->SetTitle("#phi [rad]");
 400     h3->GetYaxis()->SetTitle("#eta");
 401     h3->GetZaxis()->SetTitle("Pt [GeV/c]");
 402 }
 403
 404 //_______________________________________________________________________________
 405 void AliPHOSCorrelations::SetHistMass(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
 406 {
 407     // Set mass histograms.
 408
 409     Double_t binMult = 400;
 410     Double_t massMin = 0.0;
 411     Double_t massMax = 0.4;
 412
 413     TString spid[4]={"all","cpv","disp","both"} ;
 414
 415     TH2F * h;
 416
 417     for(Int_t ipid=0; ipid<4; ipid++)
 418     {
 419         // Real ++++++++++++++++++++++++++++++
 420
 421         fOutputContainer->Add(new TH2F(Form("%s_mpt", spid[ipid].Data() ), "Real",
 422                                 binMult, massMin, massMax,
 423                                 ptMult, ptMin, ptMax ) );
 424         h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 425         h->Sumw2();
 426         h->GetXaxis()->SetTitle("Mass [GeV]");
 427         h->GetYaxis()->SetTitle("Pt [GEV]");
 428
 429         // MIX +++++++++++++++++++++++++
 430
 431         fOutputContainer->Add(new TH2F(Form("mix_%s_mpt", spid[ipid].Data() ), "Mix",
 432                                 binMult, massMin, massMax,
 433                                 ptMult, ptMin, ptMax ) );
 434         h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 435         h->Sumw2();
 436         h->GetXaxis()->SetTitle("Mass [GeV]");
 437         h->GetYaxis()->SetTitle("Pt [GEV]");
 438     }
 439
 440     // Calibration PHOS Module Pi0peak {REAL}
 441     for(Int_t mod=1; mod<4; mod++)
 442     {
 443         fOutputContainer->Add(new TH2F(Form(  "both%d_mpt",mod), Form("Both cuts (CPV + Disp) mod[%d]",mod),
 444                                 binMult, massMin, massMax,
 445                                 ptMult, ptMin, ptMax ) );
 446         h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 447         h->Sumw2();
 448         h->GetXaxis()->SetTitle("Mass [GeV]");
 449         h->GetYaxis()->SetTitle("Pt [GEV]");
 450
 451         // Calibration PHOS Module Pi0peak {MIX}
 452         fOutputContainer->Add(new TH2F(Form(    "mix_both%d_mpt",mod), Form(" Both cuts (CPV + Disp) mod[%d]",mod),
 453                                 binMult, massMin, massMax,
 454                                 ptMult, ptMin, ptMax ) );
 455         h = static_cast<TH2F*>(fOutputContainer->Last()) ;
 456         h->Sumw2();
 457         h->GetXaxis()->SetTitle("Mass [GeV]");
 458         h->GetYaxis()->SetTitle("Pt [GEV]");
 459
 460     }
 461 }
 462
 463 //_______________________________________________________________________________
 464 void AliPHOSCorrelations::SetHistPtAssoc(const Int_t& ptMult, const Double_t& ptMin, const Double_t& ptMax) const
 465 {
 466     Double_t pi = TMath::Pi();
 467
 468     Int_t PhiMult  =  100        ;
 469     Float_t PhiMin =  -0.5*pi    ;
 470     Float_t PhiMax =  1.5*pi    ;
 471     Int_t EtaMult  =  20        ;
 472     Float_t EtaMin = -1.        ;
 473     Float_t EtaMax =  1.        ;
 474
 475     TString spid[4]={"all","cpv","disp","both"} ;
 476
 477     for (int i = 0; i<fAssocBins.GetSize()-1; i++)
 478     {
 479         for(Int_t ipid=0; ipid<4; ipid++)
 480         {
 481             // Main histo for ConsiderPi0s().
 482             fOutputContainer->Add(new TH3F(Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
 483                                     Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
 484                                     ptMult, ptMin, ptMax,
 485                                     PhiMult, PhiMin, PhiMax,
 486                                     EtaMult, EtaMin, EtaMax ) );
 487             TH3F * h = static_cast<TH3F*>(fOutputContainer->Last()) ;
 488             h->Sumw2();
 489             h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
 490             h->GetYaxis()->SetTitle("#phi [rad]");
 491             h->GetZaxis()->SetTitle("#eta");
 492
 493             // For ConsiderPi0s_MBSelection().
 494             fOutputContainer->Add(new TH3F(Form("%s_ptphieta_ptAssoc_%3.1f_MB", spid[ipid].Data(), fAssocBins.At(i+1)),
 495                                     Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
 496                                     ptMult, ptMin, ptMax,
 497                                     PhiMult, PhiMin, PhiMax,
 498                                     EtaMult, EtaMin, EtaMax ) );
 499             h = static_cast<TH3F*>(fOutputContainer->Last()) ;
 500             h->Sumw2();
 501             h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
 502             h->GetYaxis()->SetTitle("#phi [rad]");
 503             h->GetZaxis()->SetTitle("#eta");
 504
 505             // For Mixed events in ConsiderTracksMix()
 506             fOutputContainer->Add(new TH3F(Form("mix_%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
 507                                     Form("Mixed %s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), fAssocBins.At(i+1)),
 508                                     ptMult, ptMin, ptMax,
 509                                     PhiMult, PhiMin, PhiMax,
 510                                     EtaMult, EtaMin, EtaMax ) );
 511             h = static_cast<TH3F*>(fOutputContainer->Last()) ;
 512             h->Sumw2();
 513             h->GetXaxis()->SetTitle("Pt_{triger} [GEV]");
 514             h->GetYaxis()->SetTitle("#phi [rad]");
 515             h->GetZaxis()->SetTitle("#eta");
 516         }
 517     }
 518 }
 519
 520 //_______________________________________________________________________________
 521 void AliPHOSCorrelations::SetHistPHOSClusterMap()
 522 {
 523     //  Cluster X/Z/E distribution.
 524     for(int i =  0; i<5; i++)
 525     {
 526         fOutputContainer->Add(new TH3F( Form("QA_cluXZE_mod%i", i), Form("PHOS Clusters XZE distribution of module %i", i),
 527                                 70, 0, 70,
 528                                 60, 0, 60,
 529                                 200, 0, 20 ) );
 530         TH3F *h = static_cast<TH3F*>(fOutputContainer->Last()) ;
 531         h->GetXaxis()->SetTitle("X");
 532         h->GetYaxis()->SetTitle("Z");
 533         h->GetZaxis()->SetTitle("E");
 534     }
 535 }
 536
 537 //_______________________________________________________________________________
 538 void AliPHOSCorrelations::UserExec(Option_t *)
 539 {
 540     // Main loop, called for each event analyze ESD/AOD
 541     // Step 0: Event Objects
 542     LogProgress(0);
 543
 544     fEvent = InputEvent();
 545     if( ! fEvent )
 546     {
 547         AliError("Event could not be retrieved");
 548         PostData(1, fOutputContainer);
 549         return ;
 550     }
 551     LogProgress(1);
 552
 553     // Step 1(done once):
 554     if( fRunNumber != fEvent->GetRunNumber() )
 555     {
 556         fRunNumber = fEvent->GetRunNumber();
 557         fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
 558         SetGeometry();
 559         SetESDTrackCuts();
 560     }
 561     LogProgress(2);
 562
 563     if(GetPeriod().Contains("0x0"))
 564     {
 565         AliWarning("Undefined period!");
 566     }
 567
 568     // Step 2: Preparation variables for new event
 569     ZeroingVariables();
 570
 571     fEventESD = dynamic_cast<AliESDEvent*>(fEvent);
 572     fEventAOD = dynamic_cast<AliAODEvent*>(fEvent);
 573
 574     // Get Event-Handler for the trigger information
 575     fEventHandler= dynamic_cast<AliInputEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 576     if (!fEventHandler)
 577     {
 578         AliError("Could not get InputHandler");
 579         PostData(1, fOutputContainer);
 580         return; // Reject!
 581     }
 582     LogProgress(3);
 583
 584     // Step 3: Event trigger selection
 585     // fPHOSEvent, fMBEvent
 586     if( RejectTriggerMaskSelection() )
 587     {
 588         PostData(1, fOutputContainer);
 589         return; // Reject!
 590     }
 591     LogProgress(4);
 592
 593     // Step 4: Vertex
 594     // fVertex, fVertexVector, fVtxBin
 595     SetVertex();
 596     if( RejectEventVertex() )
 597     {
 598         PostData(1, fOutputContainer);
 599         return; // Reject!
 600     }
 601     LogProgress(5);
 602
 603     // Step 5: Centrality
 604     // fCentrality, fCentBin
 605     SetCentrality();
 606     if( RejectEventCentrality() )
 607     {
 608         PostData(1, fOutputContainer);
 609         return; // Reject!
 610     }
 611     LogProgress(6);
 612     if(fPHOSEvent)     FillHistogram( "hCentralityTriggerEvent",    fCentrality, fInternalRunNumber-0.5 ) ;
 613     if(fMBEvent)     FillHistogram( "hCentralityMBEvent",        fCentrality, fInternalRunNumber-0.5 ) ;
 614     FillHistogram( "hCentrality", fCentrality, fInternalRunNumber-0.5 ) ;
 615
 616     // Step 6: Reaction Plane
 617     // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
 618     EvalReactionPlane();
 619     fEMRPBin = GetRPBin();
 620
 621     // Step 7: Event Photons (PHOS Clusters) selection
 622     SelectPhotonClusters();
 623     if( ! fCaloPhotonsPHOS->GetEntriesFast() )
 624         LogSelection(kHasPHOSClusters, fInternalRunNumber);
 625
 626     // Step 8: Event Associated particles (TPC Tracks) selection
 627     SelectAccosiatedTracks();
 628     if( ! fTracksTPC->GetEntriesFast() ) LogSelection(kHasTPCTracks, fInternalRunNumber);
 629     LogSelection(kTotalSelected, fInternalRunNumber);
 630
 631     // Step 9: Fill TPC's track mask and control bining hists.
 632     FillTrackEtaPhi();
 633     FillEventBiningProperties();
 634     LogProgress(7);
 635
 636     // Step 10: Extract one most energetic pi0 candidate in this event.
 637     SelectTriggerPi0ME();
 638
 639     // Step 11: Start correlation analysis.
 640     if (fPHOSEvent)
 641     {
 642         ConsiderPi0s(); // Consider the most energetic Pi0 in this event with all tracks of this event.
 643         LogProgress(8);
 644     }
 645
 646     if(GetPeriod().Contains("13") && fMBEvent)
 647     {
 648         ConsiderPi0s_MBSelection();
 649         LogProgress(9);
 650     }
 651
 652     // Filling mixing histograms:
 653     if (fMBEvent)
 654     {
 655         ConsiderPi0sMix();      // Make background for extracting pi0 mass.
 656         ConsiderTracksMix();    // Compare only one most energetic pi0 candidate with all tracks from previous MB events.
 657
 658         // Update pull using MB events only!
 659         UpdatePhotonLists();    // Updating pull of photons.
 660         UpdateTrackLists();     // Updating pull of tracks.
 661         LogProgress(10);
 662     }
 663
 664     LogProgress(11);
 665     // Post output data.
 666     PostData(1, fOutputContainer);
 667 }
 668
 669 //_______________________________________________________________________________
 670 void AliPHOSCorrelations::SetESDTrackCuts()
 671 {
 672     if( fEventESD )
 673     {
 674         // Create ESD track cut
 675         fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
 676         fESDtrackCuts->SetRequireTPCRefit(kTRUE) ;
 677     }
 678 }
 679
 680 //_______________________________________________________________________________
 681 Int_t AliPHOSCorrelations::ConvertToInternalRunNumber(const Int_t& run) const
 682 {
 683     // Manual setup using data from logbook.
 684     if(GetPeriod().Contains("11h"))
 685     {
 686         switch(run)
 687         {
 688             case  170593 : return 179 ;
 689             case  170572 : return 178 ;
 690             case  170556 : return 177 ;
 691             case  170552 : return 176 ;
 692             case  170546 : return 175 ;
 693             case  170390 : return 174 ;
 694             case  170389 : return 173 ;
 695             case  170388 : return 172 ;
 696             case  170387 : return 171 ;
 697             case  170315 : return 170 ;
 698             case  170313 : return 169 ;
 699             case  170312 : return 168 ;
 700             case  170311 : return 167 ;
 701             case  170309 : return 166 ;
 702             case  170308 : return 165 ;
 703             case  170306 : return 164 ;
 704             case  170270 : return 163 ;
 705             case  170269 : return 162 ;
 706             case  170268 : return 161 ;
 707             case  170267 : return 160 ;
 708             case  170264 : return 159 ;
 709             case  170230 : return 158 ;
 710             case  170228 : return 157 ;
 711             case  170208 : return 156 ;
 712             case  170207 : return 155 ;
 713             case  170205 : return 154 ;
 714             case  170204 : return 153 ;
 715             case  170203 : return 152 ;
 716             case  170195 : return 151 ;
 717             case  170193 : return 150 ;
 718             case  170163 : return 149 ;
 719             case  170162 : return 148 ;
 720             case  170159 : return 147 ;
 721             case  170155 : return 146 ;
 722             case  170152 : return 145 ;
 723             case  170091 : return 144 ;
 724             case  170089 : return 143 ;
 725             case  170088 : return 142 ;
 726             case  170085 : return 141 ;
 727             case  170084 : return 140 ;
 728             case  170083 : return 139 ;
 729             case  170081 : return 138 ;
 730             case  170040 : return 137 ;
 731             case  170038 : return 136 ;
 732             case  170036 : return 135 ;
 733             case  170027 : return 134 ;
 734             case  169981 : return 133 ;
 735             case  169975 : return 132 ;
 736             case  169969 : return 131 ;
 737             case  169965 : return 130 ;
 738             case  169961 : return 129 ;
 739             case  169956 : return 128 ;
 740             case  169926 : return 127 ;
 741             case  169924 : return 126 ;
 742             case  169923 : return 125 ;
 743             case  169922 : return 124 ;
 744             case  169919 : return 123 ;
 745             case  169918 : return 122 ;
 746             case  169914 : return 121 ;
 747             case  169859 : return 120 ;
 748             case  169858 : return 119 ;
 749             case  169855 : return 118 ;
 750             case  169846 : return 117 ;
 751             case  169838 : return 116 ;
 752             case  169837 : return 115 ;
 753             case  169835 : return 114 ;
 754             case  169683 : return 113 ;
 755             case  169628 : return 112 ;
 756             case  169591 : return 111 ;
 757             case  169590 : return 110 ;
 758             case  169588 : return 109 ;
 759             case  169587 : return 108 ;
 760             case  169586 : return 107 ;
 761             case  169584 : return 106 ;
 762             case  169557 : return 105 ;
 763             case  169555 : return 104 ;
 764             case  169554 : return 103 ;
 765             case  169553 : return 102 ;
 766             case  169550 : return 101 ;
 767             case  169515 : return 100 ;
 768             case  169512 : return 99 ;
 769             case  169506 : return 98 ;
 770             case  169504 : return 97 ;
 771             case  169498 : return 96 ;
 772             case  169475 : return 95 ;
 773             case  169420 : return 94 ;
 774             case  169419 : return 93 ;
 775             case  169418 : return 92 ;
 776             case  169417 : return 91 ;
 777             case  169415 : return 90 ;
 778             case  169411 : return 89 ;
 779             case  169238 : return 88 ;
 780             case  169236 : return 87 ;
 781             case  169167 : return 86 ;
 782             case  169160 : return 85 ;
 783             case  169156 : return 84 ;
 784             case  169148 : return 83 ;
 785             case  169145 : return 82 ;
 786             case  169144 : return 81 ;
 787             case  169143 : return 80 ;
 788             case  169138 : return 79 ;
 789             case  169099 : return 78 ;
 790             case  169094 : return 77 ;
 791             case  169091 : return 76 ;
 792             case  169045 : return 75 ;
 793             case  169044 : return 74 ;
 794             case  169040 : return 73 ;
 795             case  169035 : return 72 ;
 796             case  168992 : return 71 ;
 797             case  168988 : return 70 ;
 798             case  168984 : return 69 ;
 799             case  168826 : return 68 ;
 800             case  168777 : return 67 ;
 801             case  168514 : return 66 ;
 802             case  168512 : return 65 ;
 803             case  168511 : return 64 ;
 804             case  168467 : return 63 ;
 805             case  168464 : return 62 ;
 806             case  168461 : return 61 ;
 807             case  168460 : return 60 ;
 808             case  168458 : return 59 ;
 809             case  168362 : return 58 ;
 810             case  168361 : return 57 ;
 811             case  168356 : return 56 ;
 812             case  168342 : return 55 ;
 813             case  168341 : return 54 ;
 814             case  168325 : return 53 ;
 815             case  168322 : return 52 ;
 816             case  168318 : return 51 ;
 817             case  168311 : return 50 ;
 818             case  168310 : return 49 ;
 819             case  168213 : return 48 ;
 820             case  168212 : return 47 ;
 821             case  168208 : return 46 ;
 822             case  168207 : return 45 ;
 823             case  168206 : return 44 ;
 824             case  168205 : return 43 ;
 825             case  168204 : return 42 ;
 826             case  168203 : return 41 ;
 827             case  168181 : return 40 ;
 828             case  168177 : return 39 ;
 829             case  168175 : return 38 ;
 830             case  168173 : return 37 ;
 831             case  168172 : return 36 ;
 832             case  168171 : return 35 ;
 833             case  168115 : return 34 ;
 834             case  168108 : return 33 ;
 835             case  168107 : return 32 ;
 836             case  168105 : return 31 ;
 837             case  168104 : return 30 ;
 838             case  168103 : return 29 ;
 839             case  168076 : return 28 ;
 840             case  168069 : return 27 ;
 841             case  168068 : return 26 ;
 842             case  168066 : return 25 ;
 843             case  167988 : return 24 ;
 844             case  167987 : return 23 ;
 845             case  167986 : return 22 ;
 846             case  167985 : return 21 ;
 847             case  167921 : return 20 ;
 848             case  167920 : return 19 ;
 849             case  167915 : return 18 ;
 850             case  167909 : return 17 ;
 851             case  167903 : return 16 ;
 852             case  167902 : return 15 ;
 853             case  167818 : return 14 ;
 854             case  167814 : return 13 ;
 855             case  167813 : return 12 ;
 856             case  167808 : return 11 ;
 857             case  167807 : return 10 ;
 858             case  167806 : return 9 ;
 859             case  167713 : return 8 ;
 860             case  167712 : return 7 ;
 861             case  167711 : return 6 ;
 862             case  167706 : return 5 ;
 863             case  167693 : return 4 ;
 864             case  166532 : return 3 ;
 865             case  166530 : return 2 ;
 866             case  166529 : return 1 ;
 867
 868             default : return 199;
 869         }
 870     }
 871
 872     if(GetPeriod().Contains("10h"))
 873     {
 874         switch(run)
 875         {
 876             case  139517 : return 137;
 877             case  139514 : return 136;
 878             case  139513 : return 135;
 879             case  139511 : return 134;
 880             case  139510 : return 133;
 881             case  139507 : return 132;
 882             case  139505 : return 131;
 883             case  139504 : return 130;
 884             case  139503 : return 129;
 885             case  139470 : return 128;
 886             case  139467 : return 127;
 887             case  139466 : return 126;
 888             case  139465 : return 125;
 889             case  139440 : return 124;
 890             case  139439 : return 123;
 891             case  139438 : return 122;
 892             case  139437 : return 121;
 893             case  139360 : return 120;
 894             case  139329 : return 119;
 895             case  139328 : return 118;
 896             case  139314 : return 117;
 897             case  139311 : return 116;
 898             case  139310 : return 115;
 899             case  139309 : return 114;
 900             case  139308 : return 113;
 901             case  139173 : return 112;
 902             case  139172 : return 111;
 903             case  139110 : return 110;
 904             case  139107 : return 109;
 905             case  139105 : return 108;
 906             case  139104 : return 107;
 907             case  139042 : return 106;
 908             case  139038 : return 105;
 909             case  139037 : return 104;
 910             case  139036 : return 103;
 911             case  139029 : return 102;
 912             case  139028 : return 101;
 913             case  138983 : return 100;
 914             case  138982 : return 99;
 915             case  138980 : return 98;
 916             case  138979 : return 97;
 917             case  138978 : return 96;
 918             case  138977 : return 95;
 919             case  138976 : return 94;
 920             case  138973 : return 93;
 921             case  138972 : return 92;
 922             case  138965 : return 91;
 923             case  138924 : return 90;
 924             case  138872 : return 89;
 925             case  138871 : return 88;
 926             case  138870 : return 87;
 927             case  138837 : return 86;
 928             case  138830 : return 85;
 929             case  138828 : return 84;
 930             case  138826 : return 83;
 931             case  138796 : return 82;
 932             case  138795 : return 81;
 933             case  138742 : return 80;
 934             case  138732 : return 79;
 935             case  138730 : return 78;
 936             case  138666 : return 77;
 937             case  138662 : return 76;
 938             case  138653 : return 75;
 939             case  138652 : return 74;
 940             case  138638 : return 73;
 941             case  138624 : return 72;
 942             case  138621 : return 71;
 943             case  138583 : return 70;
 944             case  138582 : return 69;
 945             case  138579 : return 68;
 946             case  138578 : return 67;
 947             case  138534 : return 66;
 948             case  138469 : return 65;
 949             case  138442 : return 64;
 950             case  138439 : return 63;
 951             case  138438 : return 62;
 952             case  138396 : return 61;
 953             case  138364 : return 60;
 954             case  138359 : return 59;
 955             case  138275 : return 58;
 956             case  138225 : return 57;
 957             case  138201 : return 56;
 958             case  138200 : return 55;
 959             case  138197 : return 54;
 960             case  138192 : return 53;
 961             case  138190 : return 52;
 962             case  138154 : return 51;
 963             case  138153 : return 50;
 964             case  138151 : return 49;
 965             case  138150 : return 48;
 966             case  138126 : return 47;
 967             case  138125 : return 46;
 968             case  137848 : return 45;
 969             case  137847 : return 44;
 970             case  137844 : return 43;
 971             case  137843 : return 42;
 972             case  137752 : return 41;
 973             case  137751 : return 40;
 974             case  137748 : return 39;
 975             case  137724 : return 38;
 976             case  137722 : return 37;
 977             case  137718 : return 36;
 978             case  137704 : return 35;
 979             case  137693 : return 34;
 980             case  137692 : return 33;
 981             case  137691 : return 32;
 982             case  137689 : return 31;
 983             case  137686 : return 30;
 984             case  137685 : return 29;
 985             case  137639 : return 28;
 986             case  137638 : return 27;
 987             case  137608 : return 26;
 988             case  137595 : return 25;
 989             case  137549 : return 24;
 990             case  137546 : return 23;
 991             case  137544 : return 22;
 992             case  137541 : return 21;
 993             case  137539 : return 20;
 994             case  137531 : return 19;
 995             case  137530 : return 18;
 996             case  137443 : return 17;
 997             case  137441 : return 16;
 998             case  137440 : return 15;
 999             case  137439 : return 14;
1000             case  137434 : return 13;
1001             case  137432 : return 12;
1002             case  137431 : return 11;
1003             case  137430 : return 10;
1004             case  137366 : return 9;
1005             case  137243 : return 8;
1006             case  137236 : return 7;
1007             case  137235 : return 6;
1008             case  137232 : return 5;
1009             case  137231 : return 4;
1010             case  137165 : return 3;
1011             case  137162 : return 2;
1012             case  137161 : return 1;
1013             default : return 199;
1014         }
1015     }
1016
1017     if(GetPeriod().Contains("13"))
1018     {
1019         switch(run)
1020         {
1021             case  195344 : return 1;
1022             case  195346 : return 2;
1023             case  195351 : return 3;
1024             case  195389 : return 4;
1025             case  195390 : return 5;
1026             case  195391 : return 6;
1027             case  195478 : return 7;
1028             case  195479 : return 8;
1029             case  195480 : return 9;
1030             case  195481 : return 10;
1031             case  195482 : return 11;
1032             case  195483 : return 12;
1033             case  195529 : return 13;
1034             case  195531 : return 14;
1035             case  195532 : return 15;
1036             case  195566 : return 16;
1037             case  195567 : return 17;
1038             case  195568 : return 18;
1039             case  195592 : return 19;
1040             case  195593 : return 20;
1041             case  195596 : return 21;
1042             case  195633 : return 22;
1043             case  195635 : return 23;
1044             case  195644 : return 24;
1045             case  195673 : return 25;
1046             case  195675 : return 26;
1047             case  195676 : return 27;
1048             case  195677 : return 28;
1049             case  195681 : return 29;
1050             case  195682 : return 30;
1051             case  195720 : return 31;
1052             case  195721 : return 32;
1053             case  195722 : return 33;
1054             case  195724 : return 34;
1055             case  195725 : return 34;
1056             case  195726 : return 35;
1057             case  195727 : return 36;
1058             case  195760 : return 37;
1059             case  195761 : return 38;
1060             case  195765 : return 39;
1061             case  195767 : return 40;
1062             case  195783 : return 41;
1063             case  195787 : return 42;
1064             case  195826 : return 43;
1065             case  195827 : return 44;
1066             case  195829 : return 45;
1067             case  195830 : return 46;
1068             case  195831 : return 47;
1069             case  195867 : return 48;
1070             case  195869 : return 49;
1071             case  195871 : return 50;
1072             case  195872 : return 51;
1073             case  195873 : return 52;
1074             case  195935 : return 53;
1075             case  195949 : return 54;
1076             case  195950 : return 55;
1077             case  195954 : return 56;
1078             case  195955 : return 57;
1079             case  195958 : return 58;
1080             case  195989 : return 59;
1081             case  195994 : return 60;
1082             case  195998 : return 61;
1083             case  196000 : return 62;
1084             case  196006 : return 63;
1085             case  196085 : return 64;
1086             case  196089 : return 65;
1087             case  196090 : return 66;
1088             case  196091 : return 67;
1089             case  196099 : return 68;
1090             case  196105 : return 69;
1091             case  196107 : return 70;
1092             case  196185 : return 71;
1093             case  196187 : return 72;
1094             case  196194 : return 73;
1095             case  196197 : return 74;
1096             case  196199 : return 75;
1097             case  196200 : return 76;
1098             case  196201 : return 77;
1099             case  196203 : return 78;
1100             case  196208 : return 79;
1101             case  196214 : return 80;
1102             case  196308 : return 81;
1103             case  196309 : return 82;
1104             case  196310 : return 83;
1105             case  196311 : return 84;
1106             case  196433 : return 85;
1107             case  196474 : return 86;
1108             case  196475 : return 87;
1109             case  196477 : return 88;
1110             case  196528 : return 89;
1111             case  196533 : return 90;
1112             case  196535 : return 91;
1113             case  196563 : return 92;
1114             case  196564 : return 93;
1115             case  196566 : return 94;
1116             case  196568 : return 95;
1117             case  196601 : return 96;
1118             case  196605 : return 97;
1119             case  196608 : return 98;
1120             case  196646 : return 99;
1121             case  196648 : return 100;
1122             case  196701 : return 101;
1123             case  196702 : return 102;
1124             case  196703 : return 103;
1125             case  196706 : return 104;
1126             case  196714 : return 105;
1127             case  196720 : return 106;
1128             case  196721 : return 107;
1129             case  196722 : return 108;
1130             case  196772 : return 109;
1131             case  196773 : return 110;
1132             case  196774 : return 111;
1133             case  196869 : return 112;
1134             case  196870 : return 113;
1135             case  196874 : return 114;
1136             case  196876 : return 115;
1137             case  196965 : return 116;
1138             case  196967 : return 117;
1139             case  196972 : return 118;
1140             case  196973 : return 119;
1141             case  196974 : return 120;
1142             case  197003 : return 121;
1143             case  197011 : return 122;
1144             case  197012 : return 123;
1145             case  197015 : return 124;
1146             case  197027 : return 125;
1147             case  197031 : return 126;
1148             case  197089 : return 127;
1149             case  197090 : return 128;
1150             case  197091 : return 129;
1151             case  197092 : return 130;
1152             case  197094 : return 131;
1153             case  197098 : return 132;
1154             case  197099 : return 133;
1155             case  197138 : return 134;
1156             case  197139 : return 135;
1157             case  197142 : return 136;
1158             case  197143 : return 137;
1159             case  197144 : return 138;
1160             case  197145 : return 139;
1161             case  197146 : return 140;
1162             case  197147 : return 140;
1163             case  197148 : return 141;
1164             case  197149 : return 142;
1165             case  197150 : return 143;
1166             case  197152 : return 144;
1167             case  197153 : return 145;
1168             case  197184 : return 146;
1169             case  197189 : return 147;
1170             case  197247 : return 148;
1171             case  197248 : return 149;
1172             case  197254 : return 150;
1173             case  197255 : return 151;
1174             case  197256 : return 152;
1175             case  197258 : return 153;
1176             case  197260 : return 154;
1177             case  197296 : return 155;
1178             case  197297 : return 156;
1179             case  197298 : return 157;
1180             case  197299 : return 158;
1181             case  197300 : return 159;
1182             case  197302 : return 160;
1183             case  197341 : return 161;
1184             case  197342 : return 162;
1185             case  197348 : return 163;
1186             case  197349 : return 164;
1187             case  197351 : return 165;
1188             case  197386 : return 166;
1189             case  197387 : return 167;
1190             case  197388 : return 168;
1191             default : return 199;
1192         }
1193     }
1194
1195     if(GetPeriod().Contains("0x0") && fDebug >= 1)
1196     {
1197         AliWarning("Period not defined");
1198     }
1199     return 1;
1200 }
1201
1202 //_______________________________________________________________________________
1203 Bool_t AliPHOSCorrelations::RejectTriggerMaskSelection()
1204 {
1205     // Analyse trigger event and reject it if it not intresting.
1206     const Bool_t REJECT = true;
1207     const Bool_t ACCEPT = false;
1208
1209     if( fDebug >= 2 )
1210         AliInfo( Form("Event passed offline phos trigger test: %s ", fEvent->GetFiredTriggerClasses().Data() ) );
1211
1212     const Int_t physSelMask = fEventHandler->IsEventSelected();
1213
1214     Bool_t isAny            = physSelMask & AliVEvent::kAny;            // to accept any trigger
1215
1216     Bool_t isPHI1           = physSelMask & AliVEvent::kPHI1;           // PHOS trigger, CINT1 suite
1217     Bool_t isPHI7           = physSelMask & AliVEvent::kPHI7;           // PHOS trigger, CINT7 suite
1218     Bool_t isPHI8           = physSelMask & AliVEvent::kPHI8;           // PHOS trigger, CINT8 suite
1219     Bool_t isCentral        = physSelMask & AliVEvent::kCentral;        // PbPb central collision trigger
1220     Bool_t isSemiCentral    = physSelMask & AliVEvent::kSemiCentral;    // PbPb semicentral collision trigger
1221     Bool_t isPHOSPb         = physSelMask & AliVEvent::kPHOSPb;         // PHOS trigger, CINT8 suite for PbPb
1222
1223     Bool_t isMB         = physSelMask & AliVEvent::kMB;                 // Minimum bias trigger, i.e. interaction trigger, offline SPD or V0 selection
1224     Bool_t isINT7       = physSelMask & AliVEvent::kINT7;               // V0AND trigger, offline V0 selection
1225     Bool_t isAnyINT     = physSelMask & AliVEvent::kAnyINT;             // to accept any interaction (aka minimum bias) trigger
1226
1227     // Other triggers
1228     Bool_t isMUON           = physSelMask & AliVEvent::kMUON;           // Muon trigger, offline SPD or V0 selection
1229     Bool_t isHighMult       = physSelMask & AliVEvent::kHighMult;       // High-multiplicity trigger (threshold defined online), offline SPD or V0 selection
1230     Bool_t isEMC1           = physSelMask & AliVEvent::kEMC1;           // EMCAL trigger
1231     Bool_t isCINT5          = physSelMask & AliVEvent::kCINT5;          // Minimum bias trigger without SPD. i.e. interaction trigger, offline V0 selection
1232     Bool_t isCMUS5          = physSelMask & AliVEvent::kCMUS5;          // Muon trigger, offline V0 selection
1233     Bool_t isMUSPB          = physSelMask & AliVEvent::kMUSPB;          // idem for PbPb
1234     Bool_t isMUSH7          = physSelMask & AliVEvent::kMUSH7;          // Muon trigger: high pt, single muon, offline V0 selection, CINT7 suite
1235     Bool_t isMUSHPB         = physSelMask & AliVEvent::kMUSHPB;         // idem for PbPb
1236     Bool_t isMUL7           = physSelMask & AliVEvent::kMUL7;           // Muon trigger: like sign dimuon, offline V0 selection, CINT7 suite
1237     Bool_t isMuonLikePB     = physSelMask & AliVEvent::kMuonLikePB;     // idem for PbPb
1238     Bool_t isMUU7           = physSelMask & AliVEvent::kMUU7;           // Muon trigger, unlike sign dimuon, offline V0 selection, CINT7 suite
1239     Bool_t isMuonUnlikePB   = physSelMask & AliVEvent::kMuonUnlikePB;   // idem for PbPb
1240     Bool_t isEMC7           = physSelMask & AliVEvent::kEMC7;           // EMCAL trigger, CINT7 suite
1241     Bool_t isEMC8           = physSelMask & AliVEvent::kEMC8;           // EMCAL trigger, CINT8 suite
1242     Bool_t isMUS7           = physSelMask & AliVEvent::kMUS7;           // Muon trigger: low pt, single muon, offline V0 selection, CINT7 suite
1243
1244     Bool_t isEMCEJE         = physSelMask & AliVEvent::kEMCEJE;                 // EMCAL jet patch trigger
1245     Bool_t isEMCEGA         = physSelMask & AliVEvent::kEMCEGA;                 // EMCAL gamma trigger
1246
1247     Bool_t isDG5                = physSelMask & AliVEvent::kDG5;                // Double gap diffractive
1248     Bool_t isZED                = physSelMask & AliVEvent::kZED;                // ZDC electromagnetic dissociation
1249     Bool_t isSPI7               = physSelMask & AliVEvent::kSPI7;               // Power interaction trigger
1250     Bool_t isSPI                = physSelMask & AliVEvent::kSPI;                // Power interaction trigger
1251     Bool_t isINT8               = physSelMask & AliVEvent::kINT8;               // CINT8 trigger: 0TVX (T0 vertex) triger
1252     Bool_t isMuonSingleLowPt8   = physSelMask & AliVEvent::kMuonSingleLowPt8;   // Muon trigger : single muon, low pt, T0 selection, CINT8 suite
1253     Bool_t isMuonSingleHighPt8  = physSelMask & AliVEvent::kMuonSingleHighPt8;  // Muon trigger : single muon, high pt, T0 selection, CINT8 suite
1254     Bool_t isMuonLikeLowPt8     = physSelMask & AliVEvent::kMuonLikeLowPt8;     // Muon trigger : like sign muon, low pt, T0 selection, CINT8 suite
1255     Bool_t isMuonUnlikeLowPt8   = physSelMask & AliVEvent::kMuonUnlikeLowPt8;   // Muon trigger : unlike sign muon, low pt, T0 selection, CINT8 suite
1256     Bool_t isMuonUnlikeLowPt0   = physSelMask & AliVEvent::kMuonUnlikeLowPt0;   // Muon trigger : unlike sign muon, low pt, no additional L0 requirement
1257     Bool_t isUserDefined        = physSelMask & AliVEvent::kUserDefined;        // Set when custom trigger classes are set in AliPhysicsSelection, offline SPD or V0 selection
1258     Bool_t isTRD                = physSelMask & AliVEvent::kTRD;                // TRD trigger
1259     Bool_t isFastOnly           = physSelMask & AliVEvent::kFastOnly;           // The fast cluster fired. This bit is set in to addition another trigger bit, e.g. kMB
1260
1261     // All input events
1262     FillHistogram("hTriggerPassedEvents", 0 );
1263     if ( isAny )        FillHistogram("hTriggerPassedEvents",  1.) ;
1264
1265     // PHOS events.
1266     if ( isPHI1 )       FillHistogram("hTriggerPassedEvents",  2. );
1267     if ( isPHI7 )       FillHistogram("hTriggerPassedEvents",  3. );
1268     if ( isPHI8 )       FillHistogram("hTriggerPassedEvents",  4. );
1269     if ( isCentral )
1270     {
1271         FillHistogram("hTriggerPassedEvents",   5. );
1272         FillHistogram("hCentralTriggerperRunNumber", fInternalRunNumber );
1273     }
1274     if ( isSemiCentral )
1275     {
1276         FillHistogram("hTriggerPassedEvents",     6. );
1277         FillHistogram("hSemiCentralTriggerperRunNumber", fInternalRunNumber );
1278     }
1279     if ( isPHOSPb )     FillHistogram("hTriggerPassedEvents",     7. );
1280
1281     // MB events.
1282     if ( isMB )
1283     {
1284         FillHistogram("hTriggerPassedEvents",     8. );
1285         FillHistogram("hMBTriggerperRunNumber", fInternalRunNumber );
1286     }
1287     if ( isINT7 )       FillHistogram("hTriggerPassedEvents",     9. );
1288     if ( isAnyINT )     FillHistogram("hTriggerPassedEvents", 10. );
1289
1290     Bool_t isTriggerEvent   = false;
1291     Bool_t isMIXEvent       = false;
1292
1293     // TODO: Set false by default.
1294     fPHOSEvent  = false ;
1295     fMBEvent    = false ;
1296
1297     // Choosing triggers for different periods.
1298     if(GetPeriod().Contains("13"))
1299     {
1300         // Working: MB + TriggerPHOS events in real events; MB only in mixed events.
1301         isTriggerEvent  =  isPHI7 || isINT7 ;
1302         isMIXEvent      =  isINT7 ;
1303     }
1304
1305     if(GetPeriod().Contains("11h"))
1306     {
1307         // Working: The same trigger in real and mixed events.
1308         isTriggerEvent = isCentral || isSemiCentral ;
1309         isMIXEvent     = isCentral || isSemiCentral ;
1310     }
1311
1312     // Select events.
1313     if(isTriggerEvent || isMIXEvent)
1314     {
1315         if ( isTriggerEvent )
1316         {
1317             FillHistogram("hTriggerPassedEvents", 17.);
1318             fPHOSEvent = true;
1319         }
1320
1321         if ( isMIXEvent )
1322         {
1323             FillHistogram("hTriggerPassedEvents", 18.);
1324             fMBEvent = true;
1325         }
1326
1327         return ACCEPT;
1328     }
1329
1330     // other events
1331     FillHistogram("hTriggerPassedEvents",  19.);
1332
1333     //if ( isAny )      FillHistogram("hTriggerPassedEvents",  1.+20) ; // Not necessary
1334     // PHOS events.
1335     if ( isPHI1 )       FillHistogram("hTriggerPassedEvents",   2.+20 );
1336     if ( isPHI7 )       FillHistogram("hTriggerPassedEvents",   3.+20 );
1337     if ( isPHI8 )       FillHistogram("hTriggerPassedEvents",   4.+20 );
1338     if ( isCentral )    FillHistogram("hTriggerPassedEvents",   5.+20 );
1339     if ( isSemiCentral )FillHistogram("hTriggerPassedEvents",   6.+20 );
1340     if ( isPHOSPb )     FillHistogram("hTriggerPassedEvents",   7.+20 );
1341     // MB events.
1342     if ( isMB )         FillHistogram("hTriggerPassedEvents",   8.+20 );
1343     if ( isINT7 )       FillHistogram("hTriggerPassedEvents",   9.+20 );
1344     if ( isAnyINT )     FillHistogram("hTriggerPassedEvents",   10.+20 );
1345     // Other rejected events
1346     if ( isMUON )               FillHistogram("hTriggerPassedEvents",     11.+20 );
1347     if ( isHighMult )           FillHistogram("hTriggerPassedEvents",     12.+20 );
1348     if ( isEMC1 )               FillHistogram("hTriggerPassedEvents",     13.+20 );
1349     if ( isCINT5 )              FillHistogram("hTriggerPassedEvents",     14.+20 );
1350     if ( isCMUS5 )              FillHistogram("hTriggerPassedEvents",     15.+20 );
1351     if ( isMUSPB )              FillHistogram("hTriggerPassedEvents",     16.+20 );
1352     if ( isMUSH7 )              FillHistogram("hTriggerPassedEvents",     17.+20 );
1353     if ( isMUSHPB )             FillHistogram("hTriggerPassedEvents",     18.+20 );
1354     if ( isMUL7 )               FillHistogram("hTriggerPassedEvents",     19.+20 );
1355     if ( isMuonLikePB )         FillHistogram("hTriggerPassedEvents",     20.+20 );
1356     if ( isMUU7 )               FillHistogram("hTriggerPassedEvents",     21.+20 );
1357     if ( isMuonUnlikePB )       FillHistogram("hTriggerPassedEvents",     22.+20 );
1358     if ( isEMC7 )               FillHistogram("hTriggerPassedEvents",     23.+20 );
1359     if ( isEMC8 )               FillHistogram("hTriggerPassedEvents",     24.+20 );
1360     if ( isMUS7 )               FillHistogram("hTriggerPassedEvents",     25.+20 );
1361     if ( isEMCEJE )             FillHistogram("hTriggerPassedEvents",     26.+20 );
1362     if ( isEMCEGA )             FillHistogram("hTriggerPassedEvents",     27.+20 );
1363     if ( isDG5 )                FillHistogram("hTriggerPassedEvents",     28.+20 );
1364     if ( isZED )                FillHistogram("hTriggerPassedEvents",     29.+20 );
1365     if ( isSPI7 )               FillHistogram("hTriggerPassedEvents",     30.+20 );
1366     if ( isSPI )                FillHistogram("hTriggerPassedEvents",     31.+20 );
1367     if ( isINT8 )               FillHistogram("hTriggerPassedEvents",     32.+20 );
1368     if ( isMuonSingleLowPt8 )   FillHistogram("hTriggerPassedEvents",     33.+20 );
1369     if ( isMuonSingleHighPt8 )  FillHistogram("hTriggerPassedEvents",     34.+20 );
1370     if ( isMuonLikeLowPt8 )     FillHistogram("hTriggerPassedEvents",     35.+20 );
1371     if ( isMuonUnlikeLowPt8 )   FillHistogram("hTriggerPassedEvents",     36.+20 );
1372     if ( isMuonUnlikeLowPt0 )   FillHistogram("hTriggerPassedEvents",     37.+20 );
1373     if ( isUserDefined )        FillHistogram("hTriggerPassedEvents",     38.+20 );
1374     if ( isTRD )                FillHistogram("hTriggerPassedEvents",     39.+20 );
1375     if ( isFastOnly )           FillHistogram("hTriggerPassedEvents",     40.+20 );
1376
1377     return REJECT;
1378 }
1379
1380 //_______________________________________________________________________________
1381 void AliPHOSCorrelations::SetVertex()
1382 {
1383     const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
1384     if( primaryVertex )
1385     {
1386         fVertex[0] = primaryVertex->GetX();
1387         fVertex[1] = primaryVertex->GetY();
1388         fVertex[2] = primaryVertex->GetZ();
1389     }
1390     else
1391     {
1392         //AliError("Event has 0x0 Primary Vertex, defaulting to origo");
1393         fVertex[0] = 0;
1394         fVertex[1] = 0;
1395         fVertex[2] = 0;
1396     }
1397     fVertexVector = TVector3(fVertex);
1398
1399     fVtxBin=GetVertexBin(fVertexVector) ;
1400 }
1401
1402 //_______________________________________________________________________________
1403 Bool_t AliPHOSCorrelations::RejectEventVertex() const
1404 {
1405     if( ! fEvent->GetPrimaryVertex() ) return true; // reject
1406     LogSelection(kHasVertex, fInternalRunNumber);
1407
1408     if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ ) return true; // reject
1409     LogSelection(kHasAbsVertex, fInternalRunNumber);
1410
1411     return false; // accept event.
1412 }
1413
1414 //_______________________________________________________________________________
1415 void AliPHOSCorrelations::SetVertexBinning()
1416 {
1417     // Define vertex bins by their edges
1418     const int nbins = fNVtxZBins+1;
1419     const double binWidth = 2*fMaxAbsVertexZ/fNVtxZBins;
1420     Double_t edges[nbins+1];
1421     for (int i = 0; i < nbins; ++i)
1422     {
1423         edges[i] = -1.*fNVtxZBins + binWidth*(double)i;
1424     }
1425
1426     TArrayD vtxEdges(nbins, edges);
1427
1428     for(int i=0; i<vtxEdges.GetSize()-1; ++i)
1429     {
1430         if(vtxEdges.At(i) > vtxEdges.At(i+1))
1431             AliFatal("edges are not sorted");
1432
1433         fVtxEdges = vtxEdges;
1434     }
1435 }
1436
1437 //_______________________________________________________________________________
1438 Int_t AliPHOSCorrelations::GetVertexBin(const TVector3&  vertexVector)
1439 {
1440     int lastBinUpperIndex = fVtxEdges.GetSize() -1;
1441     if( vertexVector.z() > fVtxEdges[lastBinUpperIndex] )
1442     {
1443         if( fDebug >= 1 )
1444             AliWarning( Form("vertex (%f) larger then upper edge of last vertex bin (%f)!", vertexVector.z(), fVtxEdges[lastBinUpperIndex]) );
1445         return lastBinUpperIndex-1;
1446     }
1447     if( vertexVector.z() < fVtxEdges[0] )
1448     {
1449         if( fDebug >= 1 )
1450         AliWarning( Form("vertex (%f) smaller then lower edge of first bin (%f)!", vertexVector.z(), fVtxEdges[0]) );
1451         return 0;
1452     }
1453
1454     fVtxBin = TMath::BinarySearch<Double_t> ( GetNumberOfVertexBins(), fVtxEdges.GetArray(), vertexVector.z() );
1455     return fVtxBin;
1456 }
1457
1458 //_______________________________________________________________________________
1459 void AliPHOSCorrelations::SetCentrality()
1460 {
1461     AliCentrality *centrality = fEvent->GetCentrality();
1462     if( centrality )
1463         fCentrality=centrality->GetCentralityPercentile(fCentralityEstimator);
1464     else
1465     {
1466         AliError("Event has 0x0 centrality");
1467         fCentrality = -1.;
1468     }
1469
1470     fCentBin = GetCentralityBin(fCentrality);
1471 }
1472
1473 //_______________________________________________________________________________
1474 Bool_t AliPHOSCorrelations::RejectEventCentrality() const
1475 {
1476     if (fCentrality<fCentralityLowLimit)
1477         return true; //reject
1478     if(fCentrality>fCentralityHightLimit)
1479         return true; //reject
1480
1481     return false;  // accept event.
1482 }
1483
1484 //_______________________________________________________________________________
1485 void AliPHOSCorrelations::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed)
1486 {
1487     // Define centrality bins by their edges
1488     for(int i=0; i<edges.GetSize()-1; ++i)
1489     {
1490         if(edges.At(i) > edges.At(i+1)) AliFatal("edges are not sorted");
1491         if( edges.GetSize() != nMixed.GetSize()+1) AliFatal("edges and nMixed don't have appropriate relative sizes");
1492
1493         fCentEdges = edges;
1494         fCentNMixed = nMixed;
1495     }
1496 }
1497
1498 //_______________________________________________________________________________
1499 Int_t AliPHOSCorrelations::GetCentralityBin(const Float_t& centralityV0M)
1500 {
1501     int lastBinUpperIndex = fCentEdges.GetSize() -1;
1502     if( centralityV0M > fCentEdges[lastBinUpperIndex] )
1503     {
1504         if( fDebug >= 1 )
1505             AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1506         return lastBinUpperIndex-1;
1507     }
1508     if( centralityV0M < fCentEdges[0] )
1509     {
1510         if( fDebug >= 1 )
1511         AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1512         return 0;
1513     }
1514
1515     fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1516     return fCentBin;
1517 }
1518
1519 //_______________________________________________________________________________
1520 void AliPHOSCorrelations::SetCentralityBorders(const Double_t& downLimit , const Double_t& upLimit )
1521 {
1522     if (downLimit < 0. || upLimit > 100 || upLimit<=downLimit)
1523         AliError( Form("Warning. Bad value of centrality borders. Setting as default: fCentralityLowLimit=%2.f, fCentralityHightLimit=%2.f", fCentralityLowLimit, fCentralityHightLimit) );
1524     else
1525     {
1526         fCentralityLowLimit     = downLimit;
1527         fCentralityHightLimit     = upLimit;
1528         AliInfo( Form("Centrality border was set as fCentralityLowLimit=%2.f, fCentralityHightLimit=%2.f", fCentralityLowLimit, fCentralityHightLimit ) );
1529     }
1530 }
1531
1532 //_______________________________________________________________________________
1533 void AliPHOSCorrelations::EvalReactionPlane()
1534 {
1535     // assigns: fHaveTPCRP and fRP
1536     // also does RP histogram fill
1537
1538     AliEventplane *eventPlane = fEvent->GetEventplane();
1539     if( ! eventPlane )
1540         { AliError("Event has no event plane"); return; }
1541
1542     Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
1543
1544     if(reactionPlaneQ>=999 || reactionPlaneQ < 0.)
1545     {
1546         //reaction plain was not defined
1547         fHaveTPCRP = kFALSE;
1548     }
1549     else
1550     {
1551         //reaction plain defined
1552         fHaveTPCRP = kTRUE;
1553     }
1554
1555     if(fHaveTPCRP)
1556         fRP = reactionPlaneQ;
1557     else
1558         fRP = 0.;
1559
1560     FillHistogram("phiRPflat",fRP,fCentrality) ;
1561 }
1562
1563 //_______________________________________________________________________________
1564 Int_t AliPHOSCorrelations::GetRPBin()
1565 {
1566     Double_t averageRP;
1567     averageRP = fRP ;         // If possible, it is better to have EP bin from TPC
1568                             // to have similar events for miximng (including jets etc)   (fRPV0A+fRPV0C+fRP) /3.;
1569     fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1570     if(fEMRPBin > (Int_t)fNEMRPBins-1)
1571         fEMRPBin = fNEMRPBins-1 ;
1572     else
1573         if(fEMRPBin < 0) fEMRPBin=0;
1574
1575     return fEMRPBin;
1576 }
1577
1578 //_______________________________________________________________________________
1579 void AliPHOSCorrelations::SelectPhotonClusters()
1580 {
1581     //Selects PHOS clusters
1582
1583     // clear (or create) array for holding events photons/clusters
1584     if(fCaloPhotonsPHOS)
1585         fCaloPhotonsPHOS->Clear();
1586     else
1587     {
1588         fCaloPhotonsPHOS = new TClonesArray("AliCaloPhoton",200);
1589         fCaloPhotonsPHOS->SetOwner();
1590     }
1591
1592     Int_t inPHOS = 0 ;
1593
1594     for (Int_t i = 0;  i < fEvent->GetNumberOfCaloClusters();  i++)
1595     {
1596         AliVCluster *clu = fEvent->GetCaloCluster(i);
1597         if (!clu->IsPHOS() || clu->E()< fMinClusterEnergy) continue; // reject cluster
1598
1599         Float_t  position[3];
1600         clu->GetPosition(position);
1601         TVector3 global(position) ;
1602         Int_t relId[4] ;
1603         fPHOSGeo->GlobalPos2RelId(global,relId) ;
1604         Int_t modPHOS  = relId[0] ;
1605         Int_t cellXPHOS = relId[2];
1606         Int_t cellZPHOS = relId[3] ;
1607
1608         Double_t distBC=clu->GetDistanceToBadChannel();
1609         if(distBC<fMinBCDistance)             continue ; // reject cluster
1610         if(clu->GetNCells() < fMinNCells)     continue ; // reject cluster
1611         if(clu->GetM02() < fMinM02)           continue ; // reject cluster
1612
1613         FillHistogram("hTOFcut", clu->GetTOF()*1.e9 );
1614
1615         if(fTOFCutEnabled)
1616         {
1617             Double_t tof = clu->GetTOF();
1618             if(TMath::Abs(tof) > fTOFCut ) continue ;
1619         }
1620         TLorentzVector lorentzMomentum;
1621         Double_t ecore = clu->GetCoreEnergy();
1622         //Double_t ecore = clu->E();
1623
1624         FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
1625
1626         Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0) ?
1627         clu->GetMomentum(lorentzMomentum, origo);
1628
1629         if(inPHOS>=fCaloPhotonsPHOS->GetSize())
1630             fCaloPhotonsPHOS->Expand(inPHOS+50) ;
1631
1632         AliCaloPhoton * ph =new((*fCaloPhotonsPHOS)[inPHOS]) AliCaloPhoton(lorentzMomentum.X(), lorentzMomentum.Py(), lorentzMomentum.Z(), lorentzMomentum.E());
1633         inPHOS++ ;
1634         ph->SetCluster(clu);
1635
1636         // Manual PHOS module number calculation
1637         /*Float_t cellId=clu->GetCellAbsId(0) ;
1638         Int_t mod = (Int_t)TMath:: Ceil(cellId/(56*64) ) ; */
1639         ph->SetModule(modPHOS) ;
1640
1641         lorentzMomentum*=ecore/lorentzMomentum.E() ;
1642
1643         //ph->SetNCells(clu->GetNCells());
1644         ph->SetMomV2(&lorentzMomentum) ;
1645         ph->SetDispBit(clu->GetDispersion() < 2.5) ;
1646         ph->SetCPVBit(clu->GetEmcCpvDistance() > 2.) ;
1647
1648         FillHistogram(Form("QA_cluXZE_mod%i", modPHOS), cellXPHOS, cellZPHOS, lorentzMomentum.E() ) ;
1649     }
1650 }
1651
1652 //_______________________________________________________________________________
1653 void AliPHOSCorrelations::SelectAccosiatedTracks()
1654 {
1655     // clear (or create) array for holding events tracks
1656     if(fTracksTPC)
1657         fTracksTPC->Clear();
1658     else
1659     {
1660         fTracksTPC = new TClonesArray("TLorentzVector",12000);
1661     }
1662     Int_t iTracks = 0 ;
1663     for (Int_t i = 0; i < fEvent->GetNumberOfTracks(); i++)
1664     {
1665
1666         AliVTrack * track = (AliVTrack*)fEvent->GetTrack(i);
1667
1668         //Select tracks under certain conditions, TPCrefit, ITSrefit ...
1669         ULong_t status = track->GetStatus();
1670         if (fTrackStatus && !((status & fTrackStatus) == fTrackStatus))
1671         {
1672             if( fDebug > 2 ) printf("\t Reject track, status != fTrackStatus\n" );
1673             continue ;
1674         }
1675
1676         if(fEventESD)
1677         {
1678             if(!SelectESDTrack((AliESDtrack*)track)) continue ; // reject track
1679         }
1680         else
1681         {
1682             if(!SelectAODTrack((AliAODTrack*)track)) continue ;    // reject track
1683         }
1684
1685         Double_t px = track->Px();
1686         Double_t py = track->Py();
1687         Double_t pz = track->Pz();
1688         Double_t e  = track->E() ;
1689
1690         if(iTracks >= fTracksTPC->GetSize())
1691             fTracksTPC->Expand(iTracks+50) ;
1692
1693         new((*fTracksTPC)[iTracks]) TLorentzVector(px, py, pz, e);
1694         iTracks++ ;
1695     }
1696 }
1697
1698 //_______________________________________________________________________________
1699 void AliPHOSCorrelations::SelectTriggerPi0ME()
1700 {
1701     const Int_t nPHOS = fCaloPhotonsPHOS->GetEntriesFast() ;
1702     for(Int_t i1 = 0; i1 < nPHOS-1; i1++)
1703     {
1704         AliCaloPhoton * ph1 = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1705         for (Int_t i2 = i1+1; i2 < nPHOS; i2++)
1706         {
1707             AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
1708             TLorentzVector p12 = *ph1 + *ph2;
1709
1710             Double_t phiTrigger = p12.Phi() ;
1711             Double_t etaTrigger = p12.Eta() ;
1712
1713             Double_t m      = p12.M() ;
1714             Double_t pt  = p12.Pt();
1715             Double_t eff = 1./GetEfficiency(pt);
1716             int mod1 = ph1->Module() ;
1717             int mod2 = ph2->Module() ;
1718
1719             FillHistogram("clu_phieta",       phiTrigger, etaTrigger );
1720             FillHistogram("clusingle_phieta", ph1->Phi(), ph1->Eta() );
1721             FillHistogram("clusingle_phieta", ph2->Phi(), ph2->Eta() );
1722
1723             FillHistogram("all_mpt", m, pt, eff );
1724
1725             if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1726             {
1727                 FillHistogram("cpv_mpt", m, pt, eff );
1728             }
1729
1730             if ( ph1->IsDispOK() && ph2->IsDispOK() )
1731             {
1732                 FillHistogram("disp_mpt", m, pt, eff );
1733                 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1734                 {
1735                     FillHistogram("both_mpt", m, pt, eff );
1736                     if(mod1 == mod2) // for each module
1737                     {
1738                         FillHistogram(Form("both%d_mpt", mod1), m, pt, eff );
1739                     }
1740                 }
1741             }
1742
1743             if(!TestMass(m,pt)) continue; //reject this pair
1744
1745             Int_t modCase = GetModCase(mod1, mod2);
1746
1747             //Now we choosing most energetic pi0.
1748             TestPi0ME(kPidAll, p12, modCase);
1749             if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1750                 TestPi0ME(kPidCPV, p12, modCase);
1751             if ( ph1->IsDispOK() && ph2->IsDispOK() )
1752             {
1753                 TestPi0ME(kPidDisp, p12, modCase);
1754                 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1755                     TestPi0ME(kPidBoth, p12, modCase);
1756             }
1757         }
1758     }
1759 }
1760
1761 //_______________________________________________________________________________
1762 void AliPHOSCorrelations::ConsiderPi0s()
1763 {
1764     TString spid[4] = {"all","cpv","disp","both"} ;
1765     // Counting number of trigger particles.
1766     for (int ipid = 0; ipid < 4; ipid++)
1767     {
1768         if (fMEExists[ipid])
1769             FillHistogram( Form("nTrigger_%s", spid[ipid].Data()), GetMEPt(ipid), 1./GetEfficiency(GetMEPt(ipid)) );
1770     }
1771
1772     // Take track's angles and compare with trigger's angles.
1773     for(Int_t i3 = 0; i3 < fTracksTPC->GetEntriesFast(); i3++)
1774     {
1775         TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i3);
1776
1777         Double_t phiAssoc = track->Phi();
1778         Double_t etaAssoc = track->Eta();
1779         Double_t ptAssoc  = track->Pt() ;
1780
1781         Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
1782         Double_t dPhi(0.), dEta(0.);
1783
1784         for (int ipid = 0; ipid < 4; ipid++)
1785         {
1786             if (GetMEExists(ipid))
1787             {
1788                 dPhi = GetMEPhi(ipid) - phiAssoc;
1789                 while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1790                 while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1791                 dEta = GetMEEta(ipid) - etaAssoc;
1792                 FillHistogram( Form("%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), ptAssocBin), GetMEPt(ipid), dPhi, dEta, 1./GetEfficiency(GetMEPt(ipid)) );
1793             }
1794         }
1795     }
1796 }
1797
1798 //_______________________________________________________________________________
1799 void AliPHOSCorrelations::ConsiderPi0s_MBSelection()
1800 {
1801     TString spid[4] = {"all","cpv","disp","both"} ;
1802     // Counting number of trigger particles.
1803     for (int ipid = 0; ipid < 4; ipid++)
1804     {
1805         if (GetMEExists(ipid))
1806         {
1807             FillHistogram( Form("nTrigger_%s_MB", spid[ipid].Data()), GetMEPt(ipid), 1./GetEfficiency(GetMEPt(ipid)) );
1808         }
1809     }
1810
1811     // Take track's angles and compare with trigger's angles.
1812     for(Int_t i3 = 0; i3 < fTracksTPC->GetEntriesFast(); i3++)
1813     {
1814         TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i3);
1815
1816         Double_t phiAssoc = track->Phi();
1817         Double_t etaAssoc = track->Eta();
1818         Double_t ptAssoc  = track->Pt();
1819
1820         Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
1821         Double_t dPhi(0.), dEta(0.);
1822
1823         for (int ipid = 0; ipid < 4; ipid++)
1824         {
1825             if (GetMEExists(ipid))
1826             {
1827                 dPhi = GetMEPhi(ipid) - phiAssoc;
1828                 while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1829                 while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1830                 dEta = GetMEEta(ipid) - etaAssoc;
1831                 FillHistogram(Form("%s_ptphieta_ptAssoc_%3.1f_MB", spid[ipid].Data(), ptAssocBin),  GetMEPt(ipid), dPhi, dEta, 1./GetEfficiency(GetMEPt(ipid)) );
1832             }
1833         }
1834     }
1835 }
1836
1837 //_______________________________________________________________________________
1838 void AliPHOSCorrelations::ConsiderPi0sMix()
1839 {
1840     TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1841     FillHistogram("hCentralityMixingProgress", fCentBin, arrayList->GetEntries() );
1842     FillHistogram("hVertexMixingProgress", fVtxBin, arrayList->GetEntries() );
1843     FillHistogram("hRPMixingProgress", fEMRPBin, arrayList->GetEntries() );
1844     for(Int_t evi = 0; evi < arrayList->GetEntries(); evi++)
1845     {
1846         TClonesArray * mixPHOS = static_cast<TClonesArray*>(arrayList->At(evi));
1847         for (Int_t i1 = 0; i1 < fCaloPhotonsPHOS->GetEntriesFast(); i1++)
1848         {
1849             AliCaloPhoton * ph1 = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1850             for(Int_t i2 = 0; i2 < mixPHOS->GetEntriesFast(); i2++)
1851             {
1852                 AliCaloPhoton * ph2 = (AliCaloPhoton*)mixPHOS->At(i2) ;
1853                 TLorentzVector p12 = *ph1 + *ph2;
1854                 Double_t m      = p12.M() ;
1855                 Double_t pt  = p12.Pt() ;
1856                 Double_t eff = 1./GetEfficiency(pt);
1857
1858                 int mod1 = ph1->Module() ;
1859                 int mod2 = ph2->Module() ;
1860
1861                 FillHistogram("mix_all_mpt", m, pt, eff);
1862                 if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1863                 {
1864                     FillHistogram("mix_cpv_mpt",m, pt, eff);
1865                 }
1866                 if ( ph1->IsDispOK() && ph2->IsDispOK() )
1867                 {
1868                     FillHistogram("mix_disp_mpt",m, pt, eff);
1869                     if ( ph1->IsCPVOK() && ph2->IsCPVOK() )
1870                     {
1871                         FillHistogram("mix_both_mpt",m, pt, eff);
1872                         if (mod1 == mod2) // for each module
1873                         {
1874                             FillHistogram(Form("mix_both%d_mpt",mod1),m, pt, eff);
1875                         }
1876                     }
1877                 }
1878             }
1879         }
1880     }
1881 }
1882
1883 //_______________________________________________________________________________
1884 void AliPHOSCorrelations::ConsiderTracksMix()
1885 {
1886     TString spid[4] = {"all","cpv","disp","both"} ;
1887
1888     TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
1889
1890     for(Int_t evi = 0; evi < arrayList->GetEntries();evi++)
1891     {
1892         TClonesArray * mixTracks = static_cast<TClonesArray*>(arrayList->At(evi));
1893         for(Int_t i3 = 0; i3 < mixTracks->GetEntriesFast(); i3++)
1894         {
1895             TLorentzVector * track = (TLorentzVector*)mixTracks->At(i3);
1896
1897             Double_t phiAssoc = track->Phi();
1898             Double_t etaAssoc = track->Eta();
1899             Double_t ptAssoc  = track->Pt();
1900
1901             Double_t ptAssocBin = GetAssocBin(ptAssoc) ;
1902
1903             Double_t ptTrigger(0.);
1904
1905             Double_t dPhi(0.), dEta(0.);
1906
1907             for (int ipid = 0; ipid < 4; ipid++)
1908             {
1909                 if (GetMEExists(ipid))
1910                 {
1911                     dPhi = GetMEPhi(ipid) - phiAssoc;
1912                     while (dPhi > 1.5*TMath::Pi()) dPhi -= 2*TMath::Pi();
1913                     while (dPhi < -.5*TMath::Pi()) dPhi += 2*TMath::Pi();
1914                     dEta = GetMEEta(ipid) - etaAssoc;
1915                     ptTrigger = GetMEPt(ipid);
1916
1917                     FillHistogram(Form("mix_%s_ptphieta_ptAssoc_%3.1f", spid[ipid].Data(), ptAssocBin), ptTrigger, dPhi, dEta, 1./GetEfficiency(ptTrigger));
1918                 }
1919             }
1920         }
1921     }
1922 }
1923
1924 //_______________________________________________________________________________
1925 TList* AliPHOSCorrelations::GetCaloPhotonsPHOSList(const UInt_t vtxBin, const UInt_t centBin, const UInt_t rpBin)
1926 {
1927     int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1928     if( fCaloPhotonsPHOSLists->At(offset) )
1929     {
1930         // list exists
1931         TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1932         return list;
1933     }
1934     else
1935     {
1936         // no list for this bin has been created, yet
1937         TList* list = new TList();
1938         fCaloPhotonsPHOSLists->AddAt(list, offset);
1939         return list;
1940     }
1941 }
1942
1943 //_______________________________________________________________________________
1944 TList* AliPHOSCorrelations::GetTracksTPCList(const UInt_t vtxBin, const UInt_t centBin, const UInt_t rpBin)
1945 {
1946     int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins + centBin * fNEMRPBins + rpBin;
1947     if( fTracksTPCLists->At(offset) )
1948     {
1949         // list exists
1950         TList* list = dynamic_cast<TList*> (fTracksTPCLists->At(offset));
1951         return list;
1952     }
1953     else
1954     {
1955         // no list for this bin has been created, yet
1956         TList* list = new TList();
1957         fTracksTPCLists->AddAt(list, offset);
1958         return list;
1959     }
1960 }
1961
1962 //_______________________________________________________________________________
1963 Double_t AliPHOSCorrelations::GetAssocBin(const Double_t& pt) const
1964 {
1965     //Calculates bin of associated particle pt.
1966     for(Int_t i=1; i<fAssocBins.GetSize(); i++)
1967     {
1968         if(pt>fAssocBins.At(i-1) && pt<fAssocBins.At(i))
1969             return fAssocBins.At(i) ;
1970     }
1971
1972     return fAssocBins.At(fAssocBins.GetSize()-1) ;
1973 }
1974
1975 //_______________________________________________________________________________
1976 void AliPHOSCorrelations::FillTrackEtaPhi() const
1977 {
1978     // Distribution TPC's tracks by angles.
1979     for (Int_t i1=0; i1<fTracksTPC->GetEntriesFast(); i1++)
1980     {
1981         TLorentzVector * track = (TLorentzVector*)fTracksTPC->At(i1);
1982         Double_t phiAssoc = track->Phi();
1983         Double_t etaAssoc = track->Eta();
1984         Double_t ptAssoc  = track->Pt() ;
1985
1986         FillHistogram( "track_phieta", phiAssoc, etaAssoc, ptAssoc );
1987     }
1988 }
1989
1990 //_______________________________________________________________________________
1991 void AliPHOSCorrelations::UpdatePhotonLists()
1992 {
1993     //Now we either add current events to stack or remove
1994     //If no photons in current event - no need to add it to mixed
1995
1996     TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1997     if( fDebug >= 2 )
1998         AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1999     if(fCaloPhotonsPHOS->GetEntriesFast()>0)
2000     {
2001         arrayList->AddFirst(fCaloPhotonsPHOS) ;
2002         fCaloPhotonsPHOS=0x0;
2003         if(arrayList->GetEntries() > fCentNMixed[fCentBin])
2004         {
2005             // Remove redundant events
2006             TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
2007             arrayList->RemoveLast() ;
2008             delete tmp;
2009         }
2010     }
2011 }
2012
2013 //_______________________________________________________________________________
2014 void AliPHOSCorrelations::UpdateTrackLists()
2015 {
2016     //Now we either add current events to stack or remove
2017     //If no photons in current event - no need to add it to mixed
2018
2019     TList * arrayList = GetTracksTPCList(fVtxBin, fCentBin, fEMRPBin);
2020
2021     if( fDebug >= 2 )
2022         AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
2023     if(fTracksTPC->GetEntriesFast()>0)
2024     {
2025         arrayList->AddFirst(fTracksTPC) ;
2026         fTracksTPC=0x0;
2027         if(arrayList->GetEntries() > fCentNMixed[fCentBin])
2028         {
2029             // Remove redundant events
2030             TClonesArray * tmp = static_cast<TClonesArray*>(arrayList->Last()) ;
2031             arrayList->RemoveLast() ;
2032             delete tmp;
2033         }
2034     }
2035 }
2036
2037 //_______________________________________________________________________________
2038 Bool_t AliPHOSCorrelations::SelectESDTrack(const AliESDtrack * t) const
2039 {
2040     // Estimate if this track can be used for correlation analisys. If all right - return "TRUE".
2041     Float_t pt = t->Pt();
2042     if( pt<0.5 || pt>20. ) return kFALSE ;
2043     if( TMath::Abs( t->Eta() ) > 0.8 ) return kFALSE;
2044     if( !fESDtrackCuts->AcceptTrack(t) ) return kFALSE ;
2045
2046     return kTRUE ;
2047 }
2048
2049 //_______________________________________________________________________________
2050 Bool_t AliPHOSCorrelations::SelectAODTrack(const AliAODTrack * t) const
2051 {
2052     // Estimate if this track can be used for correlation analisys. If all right - return "TRUE".
2053     Float_t pt = t->Pt() ;
2054     if( pt<0.5 || pt>20.     ) return kFALSE ;
2055     if( TMath::Abs(t->Eta()) > 0.8 ) return kFALSE ;
2056
2057     if(fSelectHybridTracks)
2058     {
2059         if(!t->IsHybridGlobalConstrainedGlobal())
2060         {
2061             if( fDebug > 2 ) printf("\t Reject track, IsHybridGlobalConstrainedGlobal() is FALSE\n ");
2062             return kFALSE ;
2063         }
2064     }
2065
2066     if(fSelectSPDHitTracks) //Not much sense to use with TPC only or Hybrid tracks
2067     {
2068         if(!t->HasPointOnITSLayer(0) && !t->HasPointOnITSLayer(1))
2069         {
2070             if( fDebug > 2 ) printf("\t Reject track, HasPointOnITSLayer(0) is %i and HasPointOnITSLayer(1) is %i\n", t->HasPointOnITSLayer(0), t->HasPointOnITSLayer(1) );
2071             return kFALSE ; ;
2072         }
2073     }
2074
2075     if(fSelectFractionTPCSharedClusters)
2076     {
2077         Double_t frac = Double_t(t->GetTPCnclsS()) / Double_t(t->GetTPCncls());
2078         if (frac > fCutTPCSharedClustersFraction)
2079         {
2080             if( fDebug > 2 ) printf("\t Reject track, shared cluster fraction %f > %f\n",frac, fCutTPCSharedClustersFraction);
2081             return kFALSE ;
2082         }
2083     }
2084
2085     return kTRUE ;
2086 }
2087
2088 //_______________________________________________________________________________
2089 void AliPHOSCorrelations::LogProgress(int step)
2090 {
2091     // Fill "step by step" hist
2092     FillHistogram("hTotSelEvents", step+0.5);
2093 }
2094
2095 //_______________________________________________________________________________
2096 void AliPHOSCorrelations::LogSelection(const int& step, const int& internalRunNumber) const
2097 {
2098     // the +0.5 is not realy neccisarry, but oh well... -henrik
2099     FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
2100 }
2101
2102 //_______________________________________________________________________________
2103 Bool_t AliPHOSCorrelations::TestMass(const Double_t& m, const Double_t& pt) const
2104 {
2105     // If pi0 candidate outside of mass peak, then return false.
2106     if (fUseMassWindowParametrisation && fNSigmaWidth)
2107     {
2108         // Parametrization
2109         FillHistogram("massWindow", MassMeanFunction(pt), MassSigmaFunction(pt)*fNSigmaWidth);
2110         if ( MassMeanFunction(pt)-MassSigmaFunction(pt)*fNSigmaWidth<=m && m<=MassMeanFunction(pt)+MassSigmaFunction(pt)*fNSigmaWidth )
2111         {
2112             FillHistogram("massWindowPass", 1);
2113             return true;
2114         }
2115         else
2116         {
2117             FillHistogram("massWindowPass", 2);
2118             return false;
2119         }
2120     }
2121     else
2122     {
2123         if(!fNSigmaWidth)
2124         {
2125             AliInfo( Form("fNSigmaWidth equel 0. Class will use default mass window: from %f to %f GeV.", fMassInvMeanMin, fMassInvMeanMax ) );
2126         }
2127
2128         if(fMassInvMeanMin<=m && m<=fMassInvMeanMax)
2129         {
2130             FillHistogram("massWindowPass", 3);
2131             FillHistogram("massWindow", m, 0.025); // 0.025 just for filling.
2132             return true;
2133         }
2134         else
2135         {
2136             FillHistogram("massWindowPass", 4);
2137             return false;
2138         }
2139     }
2140 }
2141
2142 //_______________________________________________________________________________
2143 Double_t AliPHOSCorrelations::MassMeanFunction(const Double_t &pt) const
2144 {
2145     // Parametrization mean of mass window
2146     return ( fMassMean[0]*pt + fMassMean[1] );
2147 }
2148
2149 //_______________________________________________________________________________
2150 Double_t AliPHOSCorrelations::MassSigmaFunction(const Double_t &pt) const
2151 {
2152     // Parametrization sigma of mass window
2153     return ( TMath::Sqrt(fMassSigma[0]*fMassSigma[0]/pt + fMassSigma[1]*fMassSigma[1]/pt/pt + fMassSigma[2]*fMassSigma[2]));
2154 }
2155
2156 //_____________________________________________________________________________
2157 void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x) const
2158 {
2159     //FillHistogram
2160     TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
2161     if(hist)
2162         hist->Fill(x) ;
2163     else
2164         AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
2165 }
2166
2167 //_____________________________________________________________________________
2168 void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x, Double_t y)const
2169 {
2170     //FillHistogram
2171     TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
2172     if(th1)
2173         th1->Fill(x, y) ;
2174     else
2175         AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
2176 }
2177
2178 //_____________________________________________________________________________
2179 void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x, Double_t y, Double_t z) const
2180 {
2181     //Fills 1D histograms with key
2182     TObject * obj = fOutputContainer->FindObject(key);
2183
2184     TH2 * th2 = dynamic_cast<TH2*> (obj);
2185     if(th2)
2186     {
2187         th2->Fill(x, y, z) ;
2188         return;
2189     }
2190
2191     TH3 * th3 = dynamic_cast<TH3*> (obj);
2192     if(th3)
2193     {
2194         th3->Fill(x, y, z) ;
2195         return;
2196     }
2197
2198     AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
2199 }
2200
2201 //_____________________________________________________________________________
2202 void AliPHOSCorrelations::FillHistogram(const char * key, Double_t x,Double_t y, Double_t z, Double_t w) const
2203 {
2204     //Fills 1D histograms with key
2205     TObject * obj = fOutputContainer->FindObject(key);
2206
2207     TH3 * th3 = dynamic_cast<TH3*> (obj);
2208     if(th3)
2209     {
2210         th3->Fill(x, y, z, w) ;
2211         return;
2212     }
2213
2214     AliError(Form("can not find histogram (of instance TH3) <%s> ",key)) ;
2215 }
2216
2217 //_____________________________________________________________________________
2218 void AliPHOSCorrelations::SetGeometry()
2219 {
2220     // Initialize the PHOS geometry
2221     //Init geometry
2222     if(!fPHOSGeo)
2223     {
2224         AliOADBContainer geomContainer("phosGeo");
2225         geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
2226         TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
2227         fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
2228         for(Int_t mod=0; mod<5; mod++)
2229         {
2230             if(!matrixes->At(mod))
2231             {
2232                 if( fDebug )
2233                 AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2234                 continue;
2235             }
2236             else
2237             {
2238                 fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
2239                 if( fDebug >1 )
2240                     AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2241             }
2242         }
2243     }
2244 }
2245
2246 //_____________________________________________________________________________
2247 Double_t AliPHOSCorrelations::GetEfficiency(const Double_t& pt) const
2248 {
2249     Double_t x = pt;
2250     //Efficiency for Both2core only!
2251     if (!fUseEfficiency)
2252         return 1.;
2253
2254     Double_t e =1.;
2255      // From 0 to 5 - 11h for different centrality.
2256      /*0: 0-5%
2257     1: 5-10%
2258     2: 10-20%
2259     3: 20-40%
2260     4: 40-60%
2261     5: 60-80%
2262     6: 0-20%
2263     7: 0-10%*/
2264     Double_t par0[9] = { -798863,       339.714,    6407.1,     -457.778,   1283.65,    -117.075,   -19.3764,   0,          0       };
2265     Double_t par1[9] = { -799344,       -1852.1,    3326.29,    -384.229,   504.046,    562.608,    130.518,    0,          0       };
2266     Double_t par2[9] = { -858904,       -1923.28,   5350.74,    -568.946,   945.497,    419.647,    101.911,    0,          0       };
2267     Double_t par3[9] = { -795652,       -1495.97,   2926.46,    -357.804,   478.961,    551.127,    128.86,     0,          0       };
2268     Double_t par4[9] = { -891951,       279626,     -123110,    -5464.75,   27470.8,    283264,     15355.1,    192762,     44828.6 };
2269     Double_t par5[9] = { -1.1094e+06,   -986.915,   2127.71,    -268.908,   375.594,    380.791,    89.4053,    0,          0       };
2270     // Double_t par6[7] = {4.86106e+09, 4.47013e+08, -1.48079e+09, 1.47233e+08, -2.62356e+08, -1.00639e+08, -2.45629e+07, 0, 0};
2271     // Double_t par7[7] = {-1.36243e+06, -26011.1, 135838, -12161.3, 24956.8, 4985.4, 1285.57, 0, 0};
2272
2273      // 8 bin for pPb13 and 0-100%
2274     Double_t par8[9] = { 6.87095e+06,   8.36553e+06,    -3.29572e+06,   2.18688e+06,    -739490,    521666,     106661,     0,  0   };
2275
2276     Double_t* pFitPoint;
2277
2278     // TODO: fix functions value below 1 GeV/c.
2279     if(x < 1.) x = 1.;
2280
2281     if(GetPeriod().Contains("11h"))
2282     {
2283         if (fCentrality <= 5)                       pFitPoint = &par0[0];
2284         if (fCentrality > 5 && fCentrality <= 10)   pFitPoint = &par1[0];
2285         if (fCentrality > 10 && fCentrality <= 20)  pFitPoint = &par2[0];
2286         if (fCentrality > 20 && fCentrality <= 40)  pFitPoint = &par3[0];
2287         if (fCentrality > 40 && fCentrality <= 60)  pFitPoint = &par4[0];
2288         if (fCentrality > 60)                       pFitPoint = &par5[0];
2289
2290         Double_t pFit[9];
2291         for (int i = 0; i < 9; ++i)
2292          {
2293              pFit[i] = *(pFitPoint+i);
2294          }
2295
2296         if (fCentrality > 40 && fCentrality <= 60)
2297             e = TMath::Exp(-(((((1.+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))+(pFit[7]*(x*(x*(x*x)))))/((((pFit[3]*x)+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x))))+(pFit[8]*(x*(x*(x*x))))))) ;
2298         else
2299             e = TMath::Exp(-((((1.+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))/(((pFit[3]*x)+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x)))))) ;
2300     }
2301     else
2302     if(GetPeriod().Contains("13"))
2303     {
2304         pFitPoint = &par8[0];
2305         Double_t pFit[9];
2306         for( int i = 0; i < 9; i++ )
2307          {
2308              pFit[i] = *(pFitPoint+i);
2309          }
2310
2311         e = TMath::Exp(-((((pFit[0]+(pFit[1]*x))+(pFit[2]*(x*x)))+(pFit[5]*(x*(x*x))))/(((1.+(pFit[3]*x))+(pFit[4]*(x*x)))+(pFit[6]*(x*(x*x)))))) ;
2312     }
2313     else
2314     {
2315         // No case
2316         AliWarning(Form("No efficiensy choise. Return 1"));
2317         e = 1.;
2318     }
2319
2320     return e;
2321 }
2322
2323 //_____________________________________________________________________________
2324 Int_t AliPHOSCorrelations::GetModCase(const Int_t &mod1, const Int_t &mod2) const
2325 {
2326     // Return modules pair namber.
2327     if(mod1 == mod2)
2328     {
2329         if(mod1 == 1) return 1;
2330         if(mod1 == 2) return 2;
2331         if(mod1 == 3) return 3;
2332     }
2333     else
2334     {
2335         if(mod1 == 1 || mod2 == 1)
2336             if(mod1 == 2 || mod2 == 2)
2337                 return 12;
2338
2339         if(mod1 == 1 || mod2 == 1)
2340             if(mod1 == 3 || mod2 == 3)
2341                 return 13;
2342         if(mod1 == 2 || mod2 == 2)
2343             if(mod1 == 3 || mod2 == 3)
2344                 return 23;
2345     }
2346
2347     AliError(Form("No choise for mod1 = %i, mod2 = %i", mod1, mod2));
2348     return 1;
2349 }
2350
2351 //_____________________________________________________________________________
2352 void AliPHOSCorrelations::TestPi0ME(const Int_t& ipid, const TLorentzVector& p12, const Int_t& modCase)
2353 {
2354     Double_t phiTrigger = p12.Phi() ;
2355     Double_t etaTrigger = p12.Eta() ;
2356     Double_t pt         = p12.Pt() ;
2357
2358     if ( GetMEExists(ipid) )
2359     {
2360         if ( pt >= GetMEPt(ipid) )
2361         {
2362             SetMEPt(ipid,pt);
2363             SetMEPhi(ipid, phiTrigger);
2364             SetMEEta(ipid, etaTrigger);
2365             SetMEModCase(ipid, modCase);
2366         }
2367     }
2368     else
2369     {
2370         SetMEPt(ipid,pt);
2371         SetMEPhi(ipid, phiTrigger);
2372         SetMEEta(ipid, etaTrigger);
2373         SetMEModCase(ipid, modCase);
2374         SetMEExists(ipid);
2375     }
2376 }
2377
2378 //_____________________________________________________________________________
2379 void AliPHOSCorrelations::ZeroingVariables()
2380 {
2381     // Set Phi, Eta, pT, modNumber andtrigger variable of moust energetic trigger particle to zero.
2382     for (int i = 0; i < 4; ++i)
2383     {
2384         fMEExists[i] = false;
2385         fMEPhi[i] = fMEEta[i] = fMEPt[i] = -99;
2386         fMEModCase[i] = 1;
2387     }
2388 }
2389
2390 //_____________________________________________________________________________
2391 void AliPHOSCorrelations::SetMassMeanParametrs(const Double_t par[2])
2392 {
2393     for (int i = 0; i < 2; ++i)
2394     {
2395         fMassMean[i] = par[i] ;
2396     }
2397 }
2398
2399 //_____________________________________________________________________________
2400 void AliPHOSCorrelations::SetMassSigmaParametrs(const Double_t par[3])
2401 {
2402     for (int i = 0; i < 3; ++i)
2403     {
2404         fMassSigma[i] = par[i] ;
2405     }
2406 }
2407
2408 //_____________________________________________________________________________
2409 void AliPHOSCorrelations::FillEventBiningProperties() const
2410 {
2411     // Fill fCentBin, fEMRPBin, fVtxBin.
2412     if(fPHOSEvent || fMBEvent)
2413     {
2414         FillHistogram( "hCentralityBining", fCentBin) ;
2415         FillHistogram( "phiRPflatBining",   fEMRPBin) ;
2416         FillHistogram( "hVertexZBining",    fVtxBin)  ;
2417         if(fPHOSEvent)
2418         {
2419             FillHistogram( "hCentralityBiningTrigger", fCentBin) ;
2420             FillHistogram( "phiRPflatBiningTrigger",   fEMRPBin) ;
2421             FillHistogram( "hVertexZBiningTrigger",    fVtxBin) ;
2422         }
2423         if(fMBEvent)
2424         {
2425             FillHistogram( "hCentralityBiningMB", fCentBin) ;
2426             FillHistogram( "phiRPflatBiningMB",   fEMRPBin) ;
2427             FillHistogram( "hVertexZBiningMB",    fVtxBin) ;
2428         }
2429     }
2430 }
2431
2432