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