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