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