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