PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0Flow.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 #include "TChain.h"
  17 #include "TTree.h"
  18 #include "TObjArray.h"
  19 #include "TF1.h"
  20 #include "TFile.h"
  21 #include "TH1F.h"
  22 #include "TH2F.h"
  23 #include "TH2I.h"
  24 #include "TH3F.h"
  25 #include "TParticle.h"
  26 #include "TCanvas.h"
  27 #include "TStyle.h"
  28 #include "TRandom.h"
  29 #include "THashList.h"
  30
  31 #include "AliAnalysisManager.h"
  32 #include "AliMCEventHandler.h"
  33 #include "AliMCEvent.h"
  34 #include "AliStack.h"
  35 #include "AliAnalysisTaskSE.h"
  36 #include "AliAnalysisTaskPi0Flow.h"
  37 #include "AliCaloPhoton.h"
  38 #include "AliPHOSGeometry.h"
  39 #include "TGeoManager.h"
  40 #include "AliPHOSEsdCluster.h"
  41 #include "AliPHOSCalibData.h"
  42 #include "AliESDEvent.h"
  43 #include "AliESDCaloCells.h"
  44 #include "AliESDVertex.h"
  45 #include "AliESDtrackCuts.h"
  46 #include "AliAODEvent.h"
  47 #include "AliLog.h"
  48 #include "AliPID.h"
  49 #include "AliCDBManager.h"
  50 #include <AliAODCaloCluster.h>
  51 #include "AliCentrality.h"
  52 #include "AliESDtrackCuts.h"
  53 #include "AliEventplane.h"
  54 #include "TProfile.h"
  55 #include "AliOADBContainer.h"
  56 #include "AliEPFlattener.h"
  57
  58 // Analysis task to fill histograms with PHOS ESD or AOD clusters and cells
  59 // Authors : Dmitri Peressounko
  60 // Date    : 28.05.2011
  61 // Modified: 03.08.2012 Henrik Qvigstad
  62 /* $Id$ */
  63
  64 ClassImp(AliAnalysisTaskPi0Flow);
  65
  66 //________________________________________________________________________
  67 Double_t rnlin(Double_t *x, Double_t * /*par*/)
  68 {
  69   //a = par[0], b = par[1].
  70   //1+a*exp(-e/b)
  71
  72 // return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0] ;
  73  return 1.015*(0.0241+1.0504*x[0]+0.000249*x[0]*x[0]) ;
  74
  75 }
  76
  77 //________________________________________________________________________
  78 AliAnalysisTaskPi0Flow::AliAnalysisTaskPi0Flow(const char *name, Period period)
  79 : AliAnalysisTaskSE(name),
  80   fCentEdges(10),
  81   fCentNMixed(10),
  82   fNEMRPBins(9),
  83   fPeriod(period),
  84   fMaxAbsVertexZ(10.),
  85   fManualV0EPCalc(false),
  86   fOutputContainer(0x0),
  87   fNonLinCorr(0),
  88   fEvent(0x0),
  89   fEventESD(0x0),
  90   fEventAOD(0x0),
  91   fMCStack(0x0),
  92   fRunNumber(-999),
  93   fInternalRunNumber(0),
  94   fPHOSGeo(0),
  95   fMultV0(0x0),
  96   fV0Cpol(0.),fV0Apol(0.),
  97   fESDtrackCuts(0x0),
  98   fPHOSCalibData(0x0),
  99   fEPcalibFileName("$ALICE_ROOT/OADB/PHOS/PHOSflat.root"),
 100   fTPCFlat(0x0),
 101   fV0AFlat(0x0),
 102   fV0CFlat(0x0),
 103   fVertexVector(),
 104   fVtxBin(0),
 105   fCentralityV0M(0.),
 106   fCentBin(0),
 107   fHaveTPCRP(0),
 108   fRP(0),
 109   fRPV0A(0),
 110   fRPV0C(0),
 111   fEMRPBin(0),
 112   fCaloPhotonsPHOS(0x0),
 113   fCaloPhotonsPHOSLists(0x0)
 114 {
 115   const int nbins = 9;
 116   Double_t edges[nbins+1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
 117   TArrayD centEdges(nbins+1, edges);
 118   Int_t nMixed[nbins] = {4,4,6,10,20,30,50,100,100};
 119   TArrayI centNMixed(nbins, nMixed);
 120   SetCentralityBinning(centEdges, centNMixed);
 121
 122   for(int mod=1; mod <= kNMod; ++mod)
 123     fModuleEnabled[mod-1] = kTRUE;
 124
 125   for(Int_t i=0;i<kNCenBins;i++){
 126     for(Int_t j=0;j<2; j++)
 127       for(Int_t k=0; k<2; k++) {
 128         fMeanQ[i][j][k]=0.;
 129         fWidthQ[i][j][k]=0.;
 130       }
 131   }
 132
 133   fVertex[0]=0; fVertex[1]=0; fVertex[2]=0;
 134
 135   // Output slots #0 write into a TH1 container
 136   DefineOutput(1,TList::Class());
 137
 138
 139
 140   // Set bad channel map
 141   char key[55] ;
 142   for(Int_t i=0; i<6; i++){
 143     snprintf(key,55,"PHOS_BadMap_mod%d",i) ;
 144     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
 145   }
 146
 147
 148   // Initialize non-linrarity correction
 149   fNonLinCorr = new TF1("nonlib",rnlin,0.,40.,0);
 150
 151
 152 }
 153 //___________________________________________________________________________
 154 AliAnalysisTaskPi0Flow::~AliAnalysisTaskPi0Flow()
 155 {
 156   delete fNonLinCorr;
 157   delete fESDtrackCuts;
 158   delete fPHOSCalibData;
 159   delete fCaloPhotonsPHOSLists;
 160   if(fTPCFlat)delete fTPCFlat;  fTPCFlat=0x0;
 161   if(fV0AFlat)delete fV0AFlat;  fV0AFlat=0x0;
 162   if(fV0CFlat)delete fV0CFlat;  fV0CFlat=0x0;
 163
 164 }
 165 //________________________________________________________________________
 166 void AliAnalysisTaskPi0Flow::UserCreateOutputObjects()
 167 {
 168   // Create histograms
 169   // Called once
 170   const Int_t nRuns=200 ;
 171
 172   // histograms
 173   if(fOutputContainer != NULL){
 174     delete fOutputContainer;
 175   }
 176   fOutputContainer = new THashList();
 177   fOutputContainer->SetOwner(kTRUE);
 178
 179   //========QA histograms=======
 180
 181   //Event selection
 182   fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 12,0.,13.,nRuns,0.,float(nRuns))) ;
 183   fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 12,0.,12.)) ;
 184
 185   //vertex distribution
 186   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
 187
 188   //Centrality
 189   fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
 190   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
 191   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
 192   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
 193   fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,10.,100,0.,100.)) ;
 194
 195
 196   //Reaction plane
 197   fOutputContainer->Add(new TH3F("hPHOSphi","cos" ,10,0.,100.,20,0.,10.,100,-TMath::Pi(),TMath::Pi()));
 198
 199   fOutputContainer->Add(new TH2F("cos2AC","RP correlation between TPC subs", 100,-1.,1.,20,0.,100.)) ;
 200   fOutputContainer->Add(new TH2F("cos2V0AC","RP correlation between VO A and C sides", 100,-1.,1.,20,0.,100.)) ;
 201   fOutputContainer->Add(new TH2F("cos2V0ATPC","RP correlation between TPC and V0A", 100,-1.,1.,20,0.,100.)) ;
 202   fOutputContainer->Add(new TH2F("cos2V0CTPC","RP correlation between TPC and V0C", 100,-1.,1.,20,0.,100.)) ;
 203
 204   fOutputContainer->Add(new TH2F("phiRP","RP distribution with TPC", 100,0.,TMath::Pi(),20,0.,100.)) ;
 205   fOutputContainer->Add(new TH2F("phiRPflat","RP distribution with TPC flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
 206   fOutputContainer->Add(new TH2F("phiRPV0A","RP distribution with V0A", 100,0.,TMath::Pi(),20,0.,100.)) ;
 207   fOutputContainer->Add(new TH2F("phiRPV0C","RP distribution with V0C", 100,0.,TMath::Pi(),20,0.,100.)) ;
 208   fOutputContainer->Add(new TH3F("phiRPV0AC","RP distribution with V0A and V0C", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
 209   fOutputContainer->Add(new TH2F("phiRPV0Aflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
 210   fOutputContainer->Add(new TH2F("phiRPV0Cflat","RP distribution with V0 flat", 100,0.,TMath::Pi(),20,0.,100.)) ;
 211   fOutputContainer->Add(new TH3F("phiRPV0ATPC","RP distribution with V0A + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
 212   fOutputContainer->Add(new TH3F("phiRPV0CTPC","RP distribution with V0C + TPC", 100,0.,TMath::Pi(),100,0.,TMath::Pi(),20,0.,100.)) ;
 213
 214
 215   //PHOS QA
 216   fOutputContainer->Add(new TH1I("hCellMultEvent"  ,"PHOS cell multiplicity per event"    ,2000,0,2000));
 217   fOutputContainer->Add(new TH1I("hCellMultEventM1","PHOS cell multiplicity per event, M1",2000,0,2000));
 218   fOutputContainer->Add(new TH1I("hCellMultEventM2","PHOS cell multiplicity per event, M2",2000,0,2000));
 219   fOutputContainer->Add(new TH1I("hCellMultEventM3","PHOS cell multiplicity per event, M3",2000,0,2000));
 220
 221   fOutputContainer->Add(new TH1F("hCellEnergy"  ,"Cell energy"            ,3000,0.,30.));
 222   fOutputContainer->Add(new TH1F("hCellEnergyM1","Cell energy in module 1",3000,0.,30.));
 223   fOutputContainer->Add(new TH1F("hCellEnergyM2","Cell energy in module 2",3000,0.,30.));
 224   fOutputContainer->Add(new TH1F("hCellEnergyM3","Cell energy in module 3",3000,0.,30.));
 225
 226   fOutputContainer->Add(new TH2F("hCellNXZM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 227   fOutputContainer->Add(new TH2F("hCellNXZM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 228   fOutputContainer->Add(new TH2F("hCellNXZM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 229   fOutputContainer->Add(new TH2F("hCellEXZM1","Cell E(X,Z), M1",64,0.5,64.5, 56,0.5,56.5));
 230   fOutputContainer->Add(new TH2F("hCellEXZM2","Cell E(X,Z), M2",64,0.5,64.5, 56,0.5,56.5));
 231   fOutputContainer->Add(new TH2F("hCellEXZM3","Cell E(X,Z), M3",64,0.5,64.5, 56,0.5,56.5));
 232
 233   //Bad Map
 234   fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 235   fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 236   fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 237
 238   fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 239   fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 240   fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 241
 242   fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 243   fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 244   fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 245
 246   fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 247   fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 248   fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 249
 250
 251   //Single photon and pi0 spectrum
 252   const Int_t nPtPhot = 400 ;
 253   const Double_t ptPhotMax = 40 ;
 254   const Int_t nM       = 500;
 255   const Double_t mMin  = 0.0;
 256   const Double_t mMax  = 1.0;
 257
 258   //PHOS calibration QA
 259   fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 260   fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 261   fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 262   fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 263   fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 264   fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 265
 266   // Histograms for different centralities
 267   const int kNPID = 16;
 268   const char* pidNames[kNPID] = {"All", "Allcore", "Allwou", "Disp", "Disp2", "Dispcore",  "Disp2core", "Dispwou", "CPV", "CPVcore", "CPV2", "CPV2core", "Both", "Bothcore", "Both2", "Both2core"};
 269   char key[55];
 270   TString name, title;
 271   for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
 272     for(Int_t ipid=0; ipid < kNPID; ipid++){
 273       name = Form("hPhot%s_cen%i", pidNames[ipid], cent );
 274       title = Form("%s clusters", pidNames[ipid]);
 275       fOutputContainer->Add(new TH1F(name.Data(), title.Data(), nPtPhot,0.,ptPhotMax));
 276
 277       name = Form("hPi0%s_cen%i", pidNames[ipid], cent );
 278       title = Form("%s clusters", pidNames[ipid]);
 279       fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 280
 281       name = Form("hSingle%s_cen%i", pidNames[ipid], cent );
 282       title = Form("%s clusters", pidNames[ipid]);
 283       fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 284
 285       name = Form("hMiPi0%s_cen%i", pidNames[ipid], cent );
 286       title = Form("%s clusters", pidNames[ipid]);
 287       fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 288
 289       name = Form("hMiSingle%s_cen%i", pidNames[ipid], cent );
 290       title = Form("%s clusters", pidNames[ipid]);
 291       fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 292
 293       // PhotPhi histograms
 294       const Int_t nPt      = 20;
 295       const Double_t xPt[21]={0.6,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.,5.5,6.,7.,8.,9.,10.,12.,14.,16.,18.,20.} ;
 296       const Int_t nPhi=10 ;
 297       Double_t xPhi[nPhi+1] ;
 298       for(Int_t i=0;i<=nPhi;i++)
 299         xPhi[i]=i*TMath::Pi() /nPhi ;
 300       const Int_t nMm=200 ;
 301       Double_t xM[nMm+1] ;
 302       for(Int_t i=0;i<=nMm;i++)
 303         xM[i]=i*0.5 /nMm;
 304       const Int_t kNPhiTitles = 3;
 305       const char* phiTitles[kNPhiTitles] = {"TPC", "V0A", "V0C"};
 306       for(Int_t iRP=0; iRP<3; iRP++){
 307         name = Form("hPhotPhi%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
 308         title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
 309         fOutputContainer->Add(new TH2F(name.Data(), title.Data(), nPt,xPt,nPhi,xPhi));
 310
 311         name = Form("hMassPt%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
 312         title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
 313         fOutputContainer->Add(new TH3F(name.Data(), title.Data(), nMm,xM,nPt,xPt,nPhi,xPhi));
 314
 315         name = Form("hMiMassPt%s%s_cen%i", phiTitles[iRP], pidNames[ipid], cent );
 316         title = Form("(M,p_{T},d#phi)_{#gamma#gamma}");
 317         fOutputContainer->Add(new TH3F(name.Data(), title.Data(), nMm,xM,nPt,xPt,nPhi,xPhi));
 318       }
 319     }
 320
 321     // a07 histograms
 322     snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
 323     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 324     snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
 325     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 326     snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
 327     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 328     snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
 329     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 330     snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
 331     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 332
 333     snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
 334     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 335     snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
 336     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 337     snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
 338     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 339     snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
 340     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 341     snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
 342     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 343   }
 344
 345   // Setup photon lists
 346   Int_t kapacity = kNVtxZBins * GetNumberOfCentralityBins() * fNEMRPBins;
 347   fCaloPhotonsPHOSLists = new TObjArray(kapacity);
 348   fCaloPhotonsPHOSLists->SetOwner();
 349
 350   PostData(1, fOutputContainer);
 351 }
 352
 353 //________________________________________________________________________
 354 void AliAnalysisTaskPi0Flow::UserExec(Option_t *)
 355 {
 356   // Main loop, called for each event
 357   // Analyze ESD/AOD
 358
 359
 360   // Step 0: Event Objects
 361   fEvent = GetEvent();
 362   fEventESD = dynamic_cast<AliESDEvent*> (fEvent);
 363   fEventAOD = dynamic_cast<AliAODEvent*> (fEvent);
 364   fMCStack = GetMCStack();
 365   LogProgress(0);
 366
 367
 368   // Step 1: Run Number, Misalignment Matrix, and Calibration
 369   // fRunNumber, fInternalRunNumber, fMultV0, fV0Cpol, fV0Apol, fMeanQ, fWidthQ
 370   if( fRunNumber != fEvent->GetRunNumber()) { // Check run number
 371     // this should run only at first call of UserExec(),
 372     // or if task runs over multiple runs, which should not occur in normal use.
 373
 374     // if run number has changed, set run variables
 375     fRunNumber = fEvent->GetRunNumber();
 376     fInternalRunNumber = ConvertToInternalRunNumber(fRunNumber);
 377     // then set misalignment and V0 calibration
 378     SetGeometry();
 379     SetMisalignment();
 380     SetV0Calibration();
 381     SetESDTrackCuts();
 382     SetPHOSCalibData();
 383     SetFlatteningData();
 384   }
 385   LogProgress(1);
 386   LogSelection(0, fInternalRunNumber);
 387
 388
 389   // Step 2: Vertex
 390   // fVertex, fVertexVector, fVtxBin
 391   SetVertex();
 392   if( RejectEventVertex() ) {
 393     PostData(1, fOutputContainer);
 394     return; // Reject!
 395   }
 396   LogProgress(2);
 397
 398 // Step 3:
 399 //   if(event->IsPileupFromSPD()){
 400 //     PostData(1, fOutputContainer);
 401 //     return; // Reject!
 402 //   }
 403   LogProgress(3);
 404
 405
 406   // Step 4: Centrality
 407   // fCentralityV0M, fCentBin
 408   SetCentrality();
 409   if( RejectCentrality() ){
 410     PostData(1, fOutputContainer);
 411     return; // Reject!
 412   }
 413   LogProgress(4);
 414
 415
 416   // Step 5: Reaction Plane
 417   // fHaveTPCRP, fRP, fRPV0A, fRPV0C, fRPBin
 418   EvalReactionPlane(); //TODO: uncomment this, or at least deal with it
 419   EvalV0ReactionPlane(); //TODO: uncomment this, or at least deal with it
 420   fEMRPBin = GetRPBin(); //TODO: uncomment this, or at least deal with it
 421   LogProgress(5);
 422
 423
 424   // Step 6: MC
 425   //  ProcessMC() ;
 426   LogProgress(6);
 427
 428
 429   // Step 7: QA PHOS cells
 430   FillPHOSCellQAHists();
 431   LogProgress(7);
 432
 433
 434   // Step 8: Event Photons (PHOS Clusters) selection
 435   SelectPhotonClusters();
 436   FillSelectedClusterHistograms();
 437   LogProgress(8);
 438
 439   if( ! fCaloPhotonsPHOS->GetEntriesFast() )
 440     return;
 441   else
 442     LogSelection(6, fInternalRunNumber);
 443
 444
 445   // Step 9: Consider pi0 (photon/cluster) pairs.
 446   ConsiderPi0s();
 447   LogProgress(9);
 448
 449   // Step 10; Mixing
 450   ConsiderPi0sMix();
 451   LogProgress(10);
 452
 453   // Step 11: Update lists
 454   UpdateLists();
 455   LogProgress(11);
 456
 457
 458   // Post output data.
 459   PostData(1, fOutputContainer);
 460 }
 461
 462 //________________________________________________________________________
 463 // void AliAnalysisTaskPi0Flow::Terminate(Option_t *)
 464 // {
 465 //   // Draw result to the screen
 466 //   // Called once at the end of the query
 467 //   // new TCanvas;
 468 //   // TH1 * hTotSelEvents = dynamic_cast<TH1*>(fOutputContainer->FindObject("hTotSelEvents"));
 469 //   // hTotSelEvents->Draw();
 470 // }
 471 //________________________________________________________________________
 472 void AliAnalysisTaskPi0Flow::SetCentralityBinning(const TArrayD& edges, const TArrayI& nMixed)
 473 {
 474   // Define centrality bins by their edges
 475   if( edges.At(0) < 0.) AliFatal("lower edge less then 0");
 476   if( 90. < edges.At(edges.GetSize()-1)  ) AliFatal("upper edge larger then 90.");
 477   for(int i=0; i<edges.GetSize()-1; ++i)
 478     if(edges.At(i) > edges.At(i+1)) AliFatal("edges are not sorted");
 479   if( edges.GetSize() != nMixed.GetSize()+1) AliFatal("edges and nMixed don't have appropriate relative sizes");
 480
 481   fCentEdges = edges;
 482   fCentNMixed = nMixed;
 483 }
 484
 485 //_____________________________________________________________________________
 486 void AliAnalysisTaskPi0Flow::SetEnablePHOSModule(int module, Bool_t enable)
 487 {
 488   if( module < 1 || kNMod < module )
 489     AliFatal(Form("PHOS Module must be between 1 and %i", kNMod));
 490   else
 491     fModuleEnabled[module-1] = enable;
 492 }
 493
 494
 495 //________________________________________________________________________
 496 void AliAnalysisTaskPi0Flow::SetPHOSBadMap(Int_t mod, TH2I* badMapHist)
 497   {
 498     if(fPHOSBadMap[mod])
 499       delete fPHOSBadMap[mod];
 500
 501     fPHOSBadMap[mod]=new TH2I(*badMapHist);
 502     if(fDebug)
 503       AliInfo(Form("Setting Bad Map Histogram  %s",fPHOSBadMap[mod]->GetName()));
 504   }
 505
 506 //________________________________________________________________________
 507 Bool_t AliAnalysisTaskPi0Flow::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
 508 {
 509   //Check if this channel belogs to the good ones
 510
 511   if(strcmp(det,"PHOS")==0){
 512     if(mod>5 || mod<1){
 513       AliError(Form("No bad map for PHOS module %d",mod)) ;
 514       return kTRUE ;
 515     }
 516     if(!fPHOSBadMap[mod]){
 517       AliError(Form("No Bad map for PHOS module %d, !fPHOSBadMap[mod]",mod)) ;
 518       return kTRUE ;
 519     }
 520     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
 521       return kFALSE ;
 522     else
 523       return kTRUE ;
 524   }
 525   else{
 526     AliError(Form("Can not find bad channels for detector %s ",det)) ;
 527   }
 528
 529   //Remove 6 noisy channels in run 139036, LHC10h
 530   if( 139036 == fRunNumber
 531     && mod==1
 532     && (ix==9||ix==10||ix==11)
 533     && (iz==45 || iz==46))
 534     return kFALSE;
 535
 536   return kTRUE ;
 537 }
 538 //_____________________________________________________________________________
 539 void AliAnalysisTaskPi0Flow::FillPHOSCellQAHists()
 540 {
 541   // Fill cell occupancy per module
 542
 543   AliVCaloCells * cells = fEvent->GetPHOSCells();
 544
 545   FillHistogram("hCenPHOSCells",fCentralityV0M,cells->GetNumberOfCells()) ;
 546   FillHistogram("hCenTrack",fCentralityV0M,fEvent->GetNumberOfTracks()) ;
 547
 548
 549   Int_t nCellModule[3] = {0,0,0};
 550   for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
 551     Int_t cellAbsId = cells->GetCellNumber(iCell);
 552     Int_t relId[4] = {0,0,0,0};
 553     fPHOSGeo->AbsToRelNumbering(cellAbsId,relId);
 554     Int_t mod1  = relId[0];
 555     Int_t cellX = relId[2];
 556     Int_t cellZ = relId[3] ;
 557     Float_t energy = cells->GetAmplitude(iCell);
 558     FillHistogram("hCellEnergy",energy);
 559     if(mod1==1) {
 560       nCellModule[0]++;
 561       FillHistogram("hCellEnergyM1",cells->GetAmplitude(iCell));
 562       FillHistogram("hCellNXZM1",cellX,cellZ,1.);
 563       FillHistogram("hCellEXZM1",cellX,cellZ,energy);
 564     }
 565     else if (mod1==2) {
 566       nCellModule[1]++;
 567       FillHistogram("hCellEnergyM2",cells->GetAmplitude(iCell));
 568       FillHistogram("hCellNXZM2",cellX,cellZ,1.);
 569       FillHistogram("hCellEXZM2",cellX,cellZ,energy);
 570     }
 571     else if (mod1==3) {
 572       nCellModule[2]++;
 573       FillHistogram("hCellEnergyM3",cells->GetAmplitude(iCell));
 574       FillHistogram("hCellNXZM3",cellX,cellZ,1.);
 575       FillHistogram("hCellEXZM3",cellX,cellZ,energy);
 576     }
 577   }
 578   FillHistogram("hCellMultEvent",nCellModule[0]+nCellModule[1]+nCellModule[2]);
 579   FillHistogram("hCellMultEventM1",nCellModule[0]);
 580   FillHistogram("hCellMultEventM2",nCellModule[1]);
 581   FillHistogram("hCellMultEventM3",nCellModule[2]);
 582
 583 }
 584 //_____________________________________________________________________________
 585 void AliAnalysisTaskPi0Flow::SelectPhotonClusters()
 586 {
 587   // clear (or create) array for holding events photons/clusters
 588   if(fCaloPhotonsPHOS)
 589     fCaloPhotonsPHOS->Clear();
 590   else {
 591     fCaloPhotonsPHOS = new TObjArray(200);
 592     fCaloPhotonsPHOS->SetOwner();
 593   }
 594
 595
 596   AliVCaloCells* cells = dynamic_cast<AliVCaloCells*> (fEvent->GetPHOSCells());
 597   for (Int_t i=0; i<fEvent->GetNumberOfCaloClusters(); i++) {
 598     AliVCluster *clu = fEvent->GetCaloCluster(i);
 599
 600     if ( !clu->IsPHOS() || clu->E()< kMinClusterEnergy) continue; // reject cluster
 601
 602
 603
 604     // check if cell/channel is good.
 605     Float_t  position[3];
 606     clu->GetPosition(position);
 607     TVector3 global(position) ;
 608     Int_t relId[4] ;
 609     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 610     Int_t mod  = relId[0] ;
 611     Int_t cellX = relId[2];
 612     Int_t cellZ = relId[3] ;
 613     if ( ! fModuleEnabled[mod-1] )
 614       continue;
 615     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 616       continue ; // reject if not.
 617
 618     Double_t distBC=clu->GetDistanceToBadChannel();
 619     if(distBC<kMinBCDistance)
 620       continue ;
 621
 622     FillHistogram("hCluEvsClu", clu->E(), clu->GetNCells()) ;
 623
 624     if(clu->GetNCells() < kMinNCells) continue ;
 625     if(clu->GetM02() < kMinM02)   continue ;
 626
 627     TLorentzVector lorentzMomentum;
 628     Double_t ecore;
 629     ecore = CoreEnergy(clu,cells);
 630
 631     //if ESD, Apply re-Calibreation
 632     Double_t origo[3] = {0,0,0}; // don't rely on event vertex, assume (0,0,0)
 633     if( fEventESD ) {
 634       AliPHOSEsdCluster cluPHOS1( *(AliESDCaloCluster*) (clu) );
 635       cluPHOS1.Recalibrate(fPHOSCalibData, static_cast<AliESDCaloCells*> (cells)); // modify the cell energies
 636       Reclusterize(&cluPHOS1) ;
 637       cluPHOS1.EvalAll(kLogWeight, fVertexVector);         // recalculate the cluster parameters
 638       cluPHOS1.SetE(fNonLinCorr->Eval(cluPHOS1.E()));// Users's nonlinearity
 639
 640       if(cluPHOS1.E()<0.3) continue; // check energy again
 641
 642       //correct misalignment
 643       TVector3 localPos;
 644       const Float_t shiftX[6]={0.,-2.3,-2.11,-1.53,0.,0.} ;
 645       const Float_t shiftZ[6]={0.,-0.4, 0.52, 0.8,0.,0.} ;
 646       fPHOSGeo->Global2Local(localPos,global,mod) ;
 647       fPHOSGeo->Local2Global(mod,localPos.X()+shiftX[mod],localPos.Z()+shiftZ[mod],global);
 648       position[0]=global.X() ;
 649       position[1]=global.Y() ;
 650       position[2]=global.Z() ;
 651       cluPHOS1.SetPosition(position);
 652
 653       cluPHOS1.GetMomentum(lorentzMomentum ,origo);
 654
 655       //TODO: Check, this may be LHC10h specific:
 656       if(mod==2) lorentzMomentum*=135.5/134.0 ;
 657       if(mod==3) lorentzMomentum*=135.5/137.2 ;
 658       if(mod==2) ecore*=135.5/134.0 ;
 659       if(mod==3) ecore*=135.5/137.2 ;
 660
 661     }
 662     else if (fEventAOD) { // is ! ESD, AOD.
 663       AliESDCaloCluster* aodCluster = (AliESDCaloCluster*) (clu);
 664       aodCluster->GetMomentum(lorentzMomentum ,origo);
 665     }
 666     else {
 667       AliError("(Calo)Cluster is neither ESD nor AOD");
 668       continue;
 669     }
 670
 671
 672     FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.);
 673     if(lorentzMomentum.E()>1.5){
 674       FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.);
 675     }
 676
 677     fCaloPhotonsPHOS->Add(new  AliCaloPhoton(lorentzMomentum.X(),lorentzMomentum.Py(),lorentzMomentum.Z(),lorentzMomentum.E()) );
 678     AliCaloPhoton * ph = (AliCaloPhoton*) fCaloPhotonsPHOS->At( fCaloPhotonsPHOS->GetLast() );
 679
 680     ph->SetModule(mod) ;
 681     lorentzMomentum*= ecore/lorentzMomentum.E() ;
 682     ph->SetMomV2(&lorentzMomentum) ;
 683     ph->SetNCells(clu->GetNCells());
 684     ph->SetDispBit(TestLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 685     //Evaluate CoreDispersion
 686     Double_t m02=0.,m20=0. ;
 687     EvalCoreLambdas(clu, cells, m02, m20) ;
 688     ph->SetDisp2Bit(TestCoreLambda(clu->E(),m20,m02)) ; //Correct order m20,m02
 689 //    ph->SetDisp2Bit(TestCoreLambda(clu->E(),clu->GetM20(),clu->GetM02())) ;
 690     if(ph->IsDispOK()){
 691       FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
 692     }
 693
 694     // Track Matching
 695     Double_t dx=clu->GetTrackDx() ;
 696     Double_t dz=clu->GetTrackDz() ;
 697     Bool_t cpvBit=kTRUE ; //No track matched by default. True means: not from charged, according to veto.
 698     Bool_t cpvBit2=kTRUE ; //More Strict criterion
 699     if( fEventESD ) {
 700
 701       TArrayI * itracks = static_cast<AliESDCaloCluster*> (clu)->GetTracksMatched() ;
 702       if(itracks->GetSize()>0){
 703         Int_t iTr = itracks->At(0);
 704         if(iTr>=0 && iTr<fEvent->GetNumberOfTracks()){
 705           AliVParticle* track = fEvent->GetTrack(iTr);
 706           Double_t pt = track->Pt() ;
 707           Short_t charge = track->Charge() ;
 708           Double_t r=TestCPV(dx, dz, pt, charge) ;
 709           cpvBit=(r>2.) ;
 710           cpvBit2=(r>4.) ;
 711         }
 712       }
 713     }
 714     else if ( fEventAOD ) {
 715       int nTracksMatched = clu->GetNTracksMatched();
 716       if(nTracksMatched > 0) {
 717         AliVTrack* track = dynamic_cast<AliVTrack*> (clu->GetTrackMatched(0));
 718         if ( track ) {
 719           Double_t pt = track->Pt();
 720           Short_t charge = track->Charge();
 721           Double_t r = TestCPV(dx, dz, pt, charge) ;
 722           cpvBit=(r>2.) ;
 723           cpvBit2=(r>4.) ;
 724         }
 725       }
 726     }
 727     ph->SetCPVBit(cpvBit) ;
 728     ph->SetCPV2Bit(cpvBit2) ;
 729     if(cpvBit){
 730       FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.);
 731     }
 732     ph->SetEMCx(float(cellX)) ;
 733     ph->SetEMCz(float(cellZ)) ;
 734     //    ph->SetLambdas(clu->GetM20(),clu->GetM02()) ;
 735     ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolde
 736
 737   }
 738   FillHistogram("hCenPHOS",fCentralityV0M, fCaloPhotonsPHOS->GetEntriesFast()) ;
 739 }
 740 //_____________________________________________________________________________
 741 void AliAnalysisTaskPi0Flow::FillSelectedClusterHistograms()
 742 {
 743   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
 744     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
 745
 746     Double_t dphiA=ph1->Phi()-fRPV0A ;
 747     while(dphiA<0)dphiA+=TMath::Pi() ;
 748     while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
 749
 750     Double_t dphiC=ph1->Phi()-fRPV0C ;
 751     while(dphiC<0)dphiC+=TMath::Pi() ;
 752     while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
 753
 754     Double_t dphiT=ph1->Phi()-fRP ;
 755     while(dphiT<0)dphiT+=TMath::Pi() ;
 756     while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
 757
 758     Double_t pt = ph1->Pt() ;
 759     Double_t ptcore = ph1->GetMomV2()->Pt() ;
 760
 761     FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),pt,dphiA) ;
 762     FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),pt,dphiC) ;
 763     if(fHaveTPCRP)
 764       FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),pt,dphiT) ;
 765     FillHistogram(Form("hPhotPhiV0AAllcore_cen%d",fCentBin),ptcore,dphiA) ;
 766     FillHistogram(Form("hPhotPhiV0CAllcore_cen%d",fCentBin),ptcore,dphiC) ;
 767     if(fHaveTPCRP)
 768       FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCentBin),ptcore,dphiT) ;
 769
 770     FillHistogram(Form("hPhotAll_cen%d",fCentBin),pt) ;
 771     FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ptcore) ;
 772     if(ph1->IsntUnfolded()){
 773       FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),pt) ;
 774       FillHistogram(Form("hPhotPhiV0AAllwou_cen%d",fCentBin),pt,dphiA) ;
 775       FillHistogram(Form("hPhotPhiV0CAllwou_cen%d",fCentBin),pt,dphiC) ;
 776       if(fHaveTPCRP)
 777         FillHistogram(Form("hPhotPhiTPCAllwou_cen%d",fCentBin),pt,dphiT) ;
 778     }
 779     if(ph1->IsCPVOK()){
 780       FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),pt,dphiA) ;
 781       FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),pt,dphiC) ;
 782       if(fHaveTPCRP)
 783         FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),pt,dphiT) ;
 784
 785       FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ptcore,dphiA) ;
 786       FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ptcore,dphiC) ;
 787       if(fHaveTPCRP)
 788         FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ptcore,dphiT) ;
 789
 790       FillHistogram(Form("hPhotCPV_cen%d",fCentBin),pt) ;
 791       FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ptcore) ;
 792     }
 793     if(ph1->IsCPV2OK()){
 794       FillHistogram(Form("hPhotPhiV0ACPV2_cen%d",fCentBin),pt,dphiA) ;
 795       FillHistogram(Form("hPhotPhiV0CCPV2_cen%d",fCentBin),pt,dphiC) ;
 796       if(fHaveTPCRP)
 797         FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCentBin),pt,dphiT) ;
 798
 799       FillHistogram(Form("hPhotPhiV0ACPV2core_cen%d",fCentBin),ptcore,dphiA) ;
 800       FillHistogram(Form("hPhotPhiV0CCPV2core_cen%d",fCentBin),ptcore,dphiC) ;
 801       if(fHaveTPCRP)
 802         FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCentBin),ptcore,dphiT) ;
 803       FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),pt) ;
 804       FillHistogram(Form("hPhotCPV2core_cen%d",fCentBin),ptcore) ;
 805     }
 806     if(ph1->IsDispOK()){
 807       FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),pt,dphiA) ;
 808       FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),pt,dphiC) ;
 809       if(fHaveTPCRP)
 810         FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),pt,dphiT) ;
 811
 812       FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ptcore,dphiA) ;
 813       FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ptcore,dphiC) ;
 814       if(fHaveTPCRP)
 815         FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ptcore,dphiT) ;
 816
 817       if(ph1->IsntUnfolded()){
 818         FillHistogram(Form("hPhotPhiV0ADispwou_cen%d",fCentBin),pt,dphiA) ;
 819         FillHistogram(Form("hPhotPhiV0CDispwou_cen%d",fCentBin),pt,dphiC) ;
 820         if(fHaveTPCRP)
 821           FillHistogram(Form("hPhotPhiTPCDispwou_cen%d",fCentBin),pt,dphiT) ;
 822
 823       }
 824       FillHistogram(Form("hPhotDisp_cen%d",fCentBin),pt) ;
 825       FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ptcore) ;
 826       if(ph1->IsntUnfolded()){
 827         FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),pt) ;
 828       }
 829       if(ph1->IsCPVOK()){
 830         FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),pt,dphiA) ;
 831         FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),pt,dphiC) ;
 832         if(fHaveTPCRP)
 833           FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),pt,dphiT) ;
 834
 835         FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ptcore,dphiA) ;
 836         FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ptcore,dphiC) ;
 837         if(fHaveTPCRP)
 838           FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ptcore,dphiT) ;
 839
 840         FillHistogram(Form("hPhotBoth_cen%d",fCentBin),pt) ;
 841         FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ptcore) ;
 842       }
 843     }
 844     if(ph1->IsDisp2OK()){
 845       FillHistogram(Form("hPhotPhiV0ADisp2_cen%d",fCentBin),pt,dphiA) ;
 846       FillHistogram(Form("hPhotPhiV0CDisp2_cen%d",fCentBin),pt,dphiC) ;
 847       if(fHaveTPCRP)
 848         FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCentBin),pt,dphiT) ;
 849       FillHistogram(Form("hPhotPhiV0ADisp2core_cen%d",fCentBin),ptcore,dphiA) ;
 850       FillHistogram(Form("hPhotPhiV0CDisp2core_cen%d",fCentBin),ptcore,dphiC) ;
 851       if(fHaveTPCRP)
 852         FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCentBin),ptcore,dphiT) ;
 853
 854       FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),pt) ;
 855       FillHistogram(Form("hPhotDisp2core_cen%d",fCentBin),ptcore) ;
 856       if(ph1->IsCPVOK()){
 857         FillHistogram(Form("hPhotPhiV0ABoth2_cen%d",fCentBin),pt,dphiA) ;
 858         FillHistogram(Form("hPhotPhiV0CBoth2_cen%d",fCentBin),pt,dphiC) ;
 859         if(fHaveTPCRP)
 860           FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCentBin),pt,dphiT) ;
 861
 862         FillHistogram(Form("hPhotPhiV0ABoth2core_cen%d",fCentBin),ptcore,dphiA) ;
 863         FillHistogram(Form("hPhotPhiV0CBoth2core_cen%d",fCentBin),ptcore,dphiC) ;
 864         if(fHaveTPCRP)
 865           FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCentBin),ptcore,dphiT) ;
 866
 867         FillHistogram(Form("hPhotBoth2_cen%d",fCentBin),pt) ;
 868         FillHistogram(Form("hPhotBoth2core_cen%d",fCentBin),ptcore) ;
 869       }
 870     }
 871   }
 872 }
 873 //_____________________________________________________________________________
 874 void AliAnalysisTaskPi0Flow::ConsiderPi0s()
 875 {
 876   char key[55];
 877   for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
 878     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
 879     for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
 880       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
 881       TLorentzVector p12  = *ph1  + *ph2;
 882       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 883       FillHistogram("hPHOSphi",fCentralityV0M,p12.Pt(),p12.Phi());
 884       Double_t dphiA=p12.Phi()-fRPV0A ;
 885       while(dphiA<0)dphiA+=TMath::Pi() ;
 886       while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
 887
 888       Double_t dphiC=p12.Phi()-fRPV0C ;
 889       while(dphiC<0)dphiC+=TMath::Pi() ;
 890       while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
 891
 892       Double_t dphiT=p12.Phi()-fRP ;
 893       while(dphiT<0)dphiT+=TMath::Pi() ;
 894       while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
 895
 896       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 897       Double_t m=p12.M() ;
 898       Double_t mcore=pv12.M() ;
 899       Double_t pt=p12.Pt() ;
 900       Double_t ptcore=pv12.Pt() ;
 901       Double_t pt1=ph1->Pt() ;
 902       Double_t pt2=ph2->Pt() ;
 903       Double_t ptcore1=ph1->GetMomV2()->Pt() ;
 904       Double_t ptcore2=ph2->GetMomV2()->Pt() ;
 905
 906       FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA) ;
 907       FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC) ;
 908       if(fHaveTPCRP)
 909         FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
 910
 911       FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
 912       FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
 913       if(fHaveTPCRP)
 914         FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
 915
 916
 917       FillHistogram(Form("hPi0All_cen%d",fCentBin),m,pt) ;
 918       FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
 919       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 920         FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),m,pt) ;
 921         FillHistogram(Form("hMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA) ;
 922         FillHistogram(Form("hMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC) ;
 923         if(fHaveTPCRP)
 924           FillHistogram(Form("hMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT) ;
 925       }
 926
 927       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt1) ;
 928       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt2) ;
 929       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore1) ;
 930       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore2) ;
 931       if(ph1->IsntUnfolded())
 932         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt1) ;
 933       if(ph2->IsntUnfolded())
 934         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt2) ;
 935       if(ph1->IsCPVOK()){
 936         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt1) ;
 937         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
 938       }
 939       if(ph2->IsCPVOK()){
 940         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt2) ;
 941         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
 942       }
 943       if(ph1->IsCPV2OK()){
 944         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt1) ;
 945         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
 946       }
 947       if(ph2->IsCPV2OK()){
 948         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt2) ;
 949         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
 950       }
 951       if(ph1->IsDispOK()){
 952         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt1) ;
 953         if(ph1->IsntUnfolded()){
 954           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt1) ;
 955         }
 956         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
 957       }
 958       if(ph2->IsDispOK()){
 959         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt2) ;
 960         if(ph1->IsntUnfolded()){
 961           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt2) ;
 962         }
 963         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
 964       }
 965       if(ph1->IsDisp2OK()){
 966         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt1) ;
 967         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
 968       }
 969       if(ph2->IsDisp2OK()){
 970         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt2) ;
 971         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
 972       }
 973       if(ph1->IsDispOK() && ph1->IsCPVOK()){
 974         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt1) ;
 975         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore1) ;
 976       }
 977       if(ph2->IsDispOK() && ph2->IsCPVOK()){
 978         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt2) ;
 979         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore2) ;
 980       }
 981       if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
 982         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt1) ;
 983         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore1) ;
 984       }
 985       if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
 986         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt2) ;
 987         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore2) ;
 988       }
 989
 990
 991       if(a<kAlphaCut){
 992         FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),m,pt) ;
 993       }
 994
 995       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 996         snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
 997         FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
 998         FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
 999         if(fHaveTPCRP)
1000           FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
1001
1002         FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1003         FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1004         if(fHaveTPCRP)
1005           FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1006
1007         FillHistogram(Form("hPi0CPV_cen%d",fCentBin),m,pt) ;
1008         FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
1009
1010         if(a<kAlphaCut){
1011           FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),m,pt) ;
1012         }
1013       }
1014       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1015         FillHistogram(Form("hMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA) ;
1016         FillHistogram(Form("hMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC) ;
1017         if(fHaveTPCRP)
1018           FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT) ;
1019         FillHistogram(Form("hMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1020         FillHistogram(Form("hMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1021         if(fHaveTPCRP)
1022           FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1023
1024         FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),m,pt) ;
1025         FillHistogram(Form("hPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
1026         if(a<kAlphaCut){
1027           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
1028         }
1029       }
1030       if(ph1->IsDispOK() && ph2->IsDispOK()){
1031         snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
1032         FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1033         FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
1034         if(fHaveTPCRP)
1035           FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
1036
1037         FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1038         FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1039         if(fHaveTPCRP)
1040           FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1041
1042         FillHistogram(Form("hPi0Disp_cen%d",fCentBin),m,pt) ;
1043         FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),mcore, ptcore) ;
1044
1045         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1046           FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),m,pt) ;
1047
1048           FillHistogram(Form("hMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA) ;
1049           FillHistogram(Form("hMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC) ;
1050           if(fHaveTPCRP)
1051             FillHistogram(Form("hMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT) ;
1052         }
1053
1054         if(a<kAlphaCut){
1055           FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),m,pt) ;
1056         }
1057         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1058           FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1059           FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
1060           if(fHaveTPCRP)
1061             FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
1062
1063           FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1064           FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1065           if(fHaveTPCRP)
1066             FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1067
1068           FillHistogram(Form("hPi0Both_cen%d",fCentBin),m,pt) ;
1069           FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),mcore,ptcore) ;
1070
1071           if(a<kAlphaCut){
1072             snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
1073             FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),m,pt) ;
1074           }
1075           if(ph1->Module()==1 && ph2->Module()==1)
1076             FillHistogram("hPi0M11",m,pt );
1077           else if(ph1->Module()==2 && ph2->Module()==2)
1078             FillHistogram("hPi0M22",m,pt );
1079           else if(ph1->Module()==3 && ph2->Module()==3)
1080             FillHistogram("hPi0M33",m,pt );
1081           else if(ph1->Module()==1 && ph2->Module()==2)
1082             FillHistogram("hPi0M12",m,pt );
1083           else if(ph1->Module()==1 && ph2->Module()==3)
1084             FillHistogram("hPi0M13",m,pt );
1085           else if(ph1->Module()==2 && ph2->Module()==3)
1086             FillHistogram("hPi0M23",m,pt );
1087
1088         }
1089
1090       }
1091
1092
1093       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1094         FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),m,pt) ;
1095         FillHistogram(Form("hPi0Disp2core_cen%d",fCentBin),mcore, ptcore) ;
1096
1097         FillHistogram(Form("hMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA) ;
1098         FillHistogram(Form("hMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC) ;
1099         if(fHaveTPCRP)
1100           FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT) ;
1101
1102         FillHistogram(Form("hMassPtV0ADisp2core_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1103         FillHistogram(Form("hMassPtV0CDisp2core_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1104         if(fHaveTPCRP)
1105           FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1106
1107         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1108           FillHistogram(Form("hMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA) ;
1109           FillHistogram(Form("hMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC) ;
1110           if(fHaveTPCRP)
1111             FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT) ;
1112
1113           FillHistogram(Form("hMassPtV0ABoth2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1114           FillHistogram(Form("hMassPtV0CBoth2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1115           if(fHaveTPCRP)
1116             FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1117
1118           FillHistogram(Form("hPi0Both2_cen%d",fCentBin),m,pt) ;
1119           FillHistogram(Form("hPi0Both2core_cen%d",fCentBin),mcore,ptcore) ;
1120         }
1121
1122       }
1123     } // end of loop i2
1124   } // end of loop i1
1125 }
1126 //_____________________________________________________________________________
1127 void AliAnalysisTaskPi0Flow::ConsiderPi0sMix()
1128 {
1129   char key[55];
1130
1131   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1132
1133   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
1134     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1135     for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
1136       TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
1137       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
1138         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
1139         TLorentzVector p12  = *ph1  + *ph2;
1140         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
1141
1142         Double_t dphiA=p12.Phi()-fRPV0A ;
1143         while(dphiA<0)dphiA+=TMath::Pi() ;
1144         while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
1145
1146         Double_t dphiC=p12.Phi()-fRPV0C ;
1147         while(dphiC<0)dphiC+=TMath::Pi() ;
1148         while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
1149
1150         Double_t dphiT=p12.Phi()-fRP ;
1151         while(dphiT<0)dphiT+=TMath::Pi() ;
1152         while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
1153
1154
1155         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
1156         Double_t m=p12.M() ;
1157         Double_t mcore=pv12.M() ;
1158         Double_t pt=p12.Pt() ;
1159         Double_t ptcore=pv12.Pt() ;
1160         Double_t pt1=ph1->Pt() ;
1161         Double_t pt2=ph2->Pt() ;
1162         Double_t ptcore1=ph1->GetMomV2()->Pt() ;
1163         Double_t ptcore2=ph2->GetMomV2()->Pt() ;
1164
1165
1166         snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
1167         FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA) ;
1168         FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC) ;
1169         if(fHaveTPCRP)
1170           FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT) ;
1171
1172         FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),mcore, ptcore, dphiA) ;
1173         FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),mcore, ptcore, dphiC) ;
1174         if(fHaveTPCRP)
1175           FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),mcore, ptcore, dphiT) ;
1176
1177         FillHistogram(Form("hMiPi0All_cen%d",fCentBin),m,pt) ;
1178         FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),mcore,ptcore) ;
1179         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1180           FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),m,pt) ;
1181           FillHistogram(Form("hMiMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA) ;
1182           FillHistogram(Form("hMiMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC) ;
1183           if(fHaveTPCRP)
1184             FillHistogram(Form("hMiMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT) ;
1185         }
1186
1187         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt1) ;
1188         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt2) ;
1189         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore1) ;
1190         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore2) ;
1191         if(ph1->IsntUnfolded())
1192           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt1) ;
1193         if(ph2->IsntUnfolded())
1194           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt2) ;
1195         if(ph1->IsCPVOK()){
1196           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt1) ;
1197           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore1) ;
1198         }
1199         if(ph2->IsCPVOK()){
1200           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt2) ;
1201           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore2) ;
1202         }
1203         if(ph1->IsCPV2OK()){
1204           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt1) ;
1205           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore1) ;
1206         }
1207         if(ph2->IsCPV2OK()){
1208           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt2) ;
1209           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore2) ;
1210         }
1211         if(ph1->IsDispOK()){
1212           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt1) ;
1213           if(ph1->IsntUnfolded()){
1214             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt1) ;
1215           }
1216           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore1) ;
1217         }
1218         if(ph2->IsDispOK()){
1219           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt2) ;
1220           if(ph1->IsntUnfolded()){
1221             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt2) ;
1222           }
1223           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore2) ;
1224         }
1225         if(ph1->IsDisp2OK()){
1226           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt1) ;
1227           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore1) ;
1228         }
1229         if(ph2->IsDisp2OK()){
1230           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt2) ;
1231           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore2) ;
1232         }
1233         if(ph1->IsDispOK() && ph1->IsCPVOK()){
1234           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1235           FillHistogram(key,m,pt1) ;
1236           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1237           FillHistogram(key,mcore,ptcore1) ;
1238         }
1239         if(ph2->IsDispOK() && ph2->IsCPVOK()){
1240           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
1241           FillHistogram(key,m,pt2) ;
1242           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin) ;
1243           FillHistogram(key,mcore,ptcore2) ;
1244         }
1245         if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1246           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt1) ;
1247           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore1) ;
1248         }
1249         if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1250           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt2) ;
1251           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore2) ;
1252         }
1253
1254
1255
1256         if(a<kAlphaCut){
1257           FillHistogram(Form("hMiPi0All_a07_cen%d",fCentBin),m,pt) ;
1258         }
1259         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1260           FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA) ;
1261           FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC) ;
1262           if(fHaveTPCRP)
1263             FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT) ;
1264
1265           FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),mcore, ptcore,dphiA) ;
1266           FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),mcore, ptcore,dphiC) ;
1267           if(fHaveTPCRP)
1268             FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),mcore, ptcore,dphiT) ;
1269
1270           FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),m,pt) ;
1271           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),mcore, ptcore) ;
1272
1273           if(a<kAlphaCut){
1274             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),m,pt) ;
1275           }
1276         }
1277         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1278           FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),m,pt) ;
1279           FillHistogram(Form("hMiPi0CPV2core_cen%d",fCentBin),mcore, ptcore) ;
1280
1281           FillHistogram(Form("hMiMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA) ;
1282           FillHistogram(Form("hMiMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC) ;
1283           if(fHaveTPCRP)
1284             FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT) ;
1285           FillHistogram(Form("hMiMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA) ;
1286           FillHistogram(Form("hMiMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC) ;
1287           if(fHaveTPCRP)
1288             FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT) ;
1289
1290           if(a<kAlphaCut){
1291             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),m,pt) ;
1292           }
1293         }
1294         if(ph1->IsDispOK() && ph2->IsDispOK()){
1295           FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA) ;
1296           FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC) ;
1297           if(fHaveTPCRP)
1298             FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT) ;
1299
1300           FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1301           FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1302           if(fHaveTPCRP)
1303             FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1304
1305
1306           FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),m,pt) ;
1307           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1308           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1309             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),m,pt) ;
1310             FillHistogram(Form("hMiMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA) ;
1311             FillHistogram(Form("hMiMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC) ;
1312             if(fHaveTPCRP)
1313               FillHistogram(Form("hMiMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT) ;
1314           }
1315
1316           if(a<kAlphaCut){
1317             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),m,pt) ;
1318           }
1319           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1320             FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA) ;
1321             FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC) ;
1322             if(fHaveTPCRP)
1323               FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT) ;
1324
1325             FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1326             FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1327             if(fHaveTPCRP)
1328               FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1329
1330             FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),m,pt) ;
1331             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1332
1333             if(a<kAlphaCut){
1334               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),m,pt) ;
1335             }
1336           }
1337         }
1338
1339         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1340           FillHistogram(Form("hMiMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA) ;
1341           FillHistogram(Form("hMiMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC) ;
1342           if(fHaveTPCRP)
1343             FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT) ;
1344
1345           FillHistogram(Form("hMiMassPtV0ADisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1346           FillHistogram(Form("hMiMassPtV0CDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1347           if(fHaveTPCRP)
1348             FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1349
1350
1351           FillHistogram(Form("hMiPi0Disp2_cen%d",fCentBin),m,pt) ;
1352           FillHistogram(Form("hMiPi0Disp2core_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1353
1354           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1355             FillHistogram(Form("hMiMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA) ;
1356             FillHistogram(Form("hMiMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC) ;
1357             if(fHaveTPCRP)
1358               FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT) ;
1359
1360             FillHistogram(Form("hMiMassPtV0ABoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA) ;
1361             FillHistogram(Form("hMiMassPtV0CBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC) ;
1362             if(fHaveTPCRP)
1363               FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT) ;
1364
1365             FillHistogram(Form("hMiPi0Both2_cen%d",fCentBin),m,pt) ;
1366             FillHistogram(Form("hMiPi0Both2core_cen%d",fCentBin),pv12.M(),pv12.Pt()) ;
1367
1368           }
1369         }
1370       } // end of loop i2
1371     }
1372   } // end of loop i1
1373 }
1374 //_____________________________________________________________________________
1375 void AliAnalysisTaskPi0Flow::UpdateLists()
1376 {
1377   //Now we either add current events to stack or remove
1378   //If no photons in current event - no need to add it to mixed
1379
1380   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
1381   if( fDebug >= 2 )
1382     AliInfo( Form("fCentBin=%d, fCentNMixed[]=%d",fCentBin,fCentNMixed[fCentBin]) );
1383   if(fCaloPhotonsPHOS->GetEntriesFast()>0){
1384     arrayList->AddFirst(fCaloPhotonsPHOS) ;
1385     fCaloPhotonsPHOS=0;
1386     if(arrayList->GetEntries() > fCentNMixed[fCentBin]){ // Remove redundant events
1387       TObjArray * tmp = static_cast<TObjArray*>(arrayList->Last()) ;
1388       arrayList->RemoveLast() ;
1389       delete tmp ; // TODO: may conflict with delete done by list being owner.
1390     }
1391   }
1392   else
1393     fCaloPhotonsPHOS->Clear(); // TODO: redundant???
1394 }
1395 //_____________________________________________________________________________
1396 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x)const{
1397   //FillHistogram
1398   TH1 * hist = dynamic_cast<TH1*>(fOutputContainer->FindObject(key)) ;
1399   if(hist)
1400     hist->Fill(x) ;
1401   else
1402     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
1403 }
1404 //_____________________________________________________________________________
1405 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y)const{
1406   //FillHistogram
1407   TH1 * th1 = dynamic_cast<TH1*> (fOutputContainer->FindObject(key));
1408   if(th1)
1409     th1->Fill(x, y) ;
1410   else
1411     AliError(Form("can not find histogram (of instance TH1) <%s> ",key)) ;
1412 }
1413
1414 //_____________________________________________________________________________
1415 void AliAnalysisTaskPi0Flow::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1416   //Fills 1D histograms with key
1417   TObject * obj = fOutputContainer->FindObject(key);
1418
1419   TH2 * th2 = dynamic_cast<TH2*> (obj);
1420   if(th2) {
1421     th2->Fill(x, y, z) ;
1422     return;
1423   }
1424
1425   TH3 * th3 = dynamic_cast<TH3*> (obj);
1426   if(th3) {
1427     th3->Fill(x, y, z) ;
1428     return;
1429   }
1430
1431   AliError(Form("can not find histogram (of instance TH2) <%s> ",key)) ;
1432 }
1433
1434 //_____________________________________________________________________________
1435 AliVEvent* AliAnalysisTaskPi0Flow::GetEvent()
1436 {
1437   fEvent = InputEvent();
1438   if( ! fEvent ) {
1439     AliError("Event could not be retrieved");
1440     PostData(1, fOutputContainer);
1441   }
1442   return fEvent;
1443 }
1444
1445
1446 //___________________________________________________________________________
1447 AliStack* AliAnalysisTaskPi0Flow::GetMCStack()
1448 {
1449   fMCStack = 0;
1450   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
1451   if(eventHandler){
1452     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
1453     if( mcEventHandler)
1454       fMCStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
1455   }
1456   return fMCStack;
1457 }
1458
1459 //___________________________________________________________________________
1460 Int_t AliAnalysisTaskPi0Flow::GetCentralityBin(Float_t centralityV0M)
1461 {
1462  /* fCentBin=1+Int_t(centralityV0M/100. *kNCenBins) ;
1463   if(centralityV0M < 5. || fCentBin < 0)
1464    fCentBin=0 ;
1465   if(fCentBin > kNCenBins-1)
1466     fCentBin = kNCenBins-1 ;
1467  */
1468   int lastBinUpperIndex = fCentEdges.GetSize() -1;
1469   if( centralityV0M > fCentEdges[lastBinUpperIndex] ) {
1470     if( fDebug >= 1 )
1471       AliWarning( Form("centrality (%f) larger then upper edge of last centrality bin (%f)!", centralityV0M, fCentEdges[lastBinUpperIndex]) );
1472     return lastBinUpperIndex-1;
1473   }
1474   if( centralityV0M < fCentEdges[0] ) {
1475     if( fDebug >= 1 )
1476       AliWarning( Form("centrality (%f) smaller then lower edge of first bin (%f)!", centralityV0M, fCentEdges[0]) );
1477     return 0;
1478   }
1479
1480   fCentBin = TMath::BinarySearch<Double_t> ( GetNumberOfCentralityBins(), fCentEdges.GetArray(), centralityV0M );
1481   return fCentBin;
1482 }
1483
1484 //___________________________________________________________________________
1485 Int_t AliAnalysisTaskPi0Flow::GetRPBin()
1486 {
1487   Double_t averageRP;
1488   if(fHaveTPCRP)
1489     averageRP = fRP ;// If possible, it is better to have EP bin from TPC
1490                      // to have similar events for miximng (including jets etc)   (fRPV0A+fRPV0C+fRP) /3.;
1491   else
1492     averageRP = (fRPV0A+fRPV0C) /2.;
1493
1494   fEMRPBin = Int_t(fNEMRPBins*(averageRP)/TMath::Pi());
1495
1496   if(fEMRPBin> (Int_t) fNEMRPBins-1)
1497     fEMRPBin=fNEMRPBins-1 ;
1498   else if(fEMRPBin<0)
1499     fEMRPBin=0;
1500
1501   if ( fDebug >= 2 )
1502     AliInfo(Form("Event Mixing Reaction Plane bin is: %d", fEMRPBin));
1503
1504   return fEMRPBin;
1505 }
1506
1507
1508 //_____________________________________________________________________________
1509 void AliAnalysisTaskPi0Flow::LogProgress(int step)
1510 {
1511   if(fDebug >= 2) {
1512     AliInfo(Form("step %d completed", step));
1513   }
1514   // the +0.5 is not realy neccisarry, but oh well... -henrik
1515   //FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1516   //FillHistogram("hTotSelEvents", step+0.5);
1517 }
1518
1519 void AliAnalysisTaskPi0Flow::LogSelection(int step, int internalRunNumber)
1520 {
1521   // if(fDebug > 1) {
1522   //   AliInfo(Form("step %d completed", step));
1523   // }
1524   // the +0.5 is not realy neccisarry, but oh well... -henrik
1525   FillHistogram("hSelEvents", step+0.5, internalRunNumber-0.5);
1526   FillHistogram("hTotSelEvents", step+0.5);
1527 }
1528
1529
1530 //___________________________________________________________________________
1531 Int_t AliAnalysisTaskPi0Flow::ConvertToInternalRunNumber(Int_t run){
1532   if( kLHC11h == fPeriod ) {
1533     switch(run){
1534     case  170593 : return 179 ;
1535     case  170572 : return 178 ;
1536     case  170556 : return 177 ;
1537     case  170552 : return 176 ;
1538     case  170546 : return 175 ;
1539     case  170390 : return 174 ;
1540     case  170389 : return 173 ;
1541     case  170388 : return 172 ;
1542     case  170387 : return 171 ;
1543     case  170315 : return 170 ;
1544     case  170313 : return 169 ;
1545     case  170312 : return 168 ;
1546     case  170311 : return 167 ;
1547     case  170309 : return 166 ;
1548     case  170308 : return 165 ;
1549     case  170306 : return 164 ;
1550     case  170270 : return 163 ;
1551     case  170269 : return 162 ;
1552     case  170268 : return 161 ;
1553     case  170267 : return 160 ;
1554     case  170264 : return 159 ;
1555     case  170230 : return 158 ;
1556     case  170228 : return 157 ;
1557     case  170208 : return 156 ;
1558     case  170207 : return 155 ;
1559     case  170205 : return 154 ;
1560     case  170204 : return 153 ;
1561     case  170203 : return 152 ;
1562     case  170195 : return 151 ;
1563     case  170193 : return 150 ;
1564     case  170163 : return 149 ;
1565     case  170162 : return 148 ;
1566     case  170159 : return 147 ;
1567     case  170155 : return 146 ;
1568     case  170152 : return 145 ;
1569     case  170091 : return 144 ;
1570     case  170089 : return 143 ;
1571     case  170088 : return 142 ;
1572     case  170085 : return 141 ;
1573     case  170084 : return 140 ;
1574     case  170083 : return 139 ;
1575     case  170081 : return 138 ;
1576     case  170040 : return 137 ;
1577     case  170038 : return 136 ;
1578     case  170036 : return 135 ;
1579     case  170027 : return 134 ;
1580     case  169981 : return 133 ;
1581     case  169975 : return 132 ;
1582     case  169969 : return 131 ;
1583     case  169965 : return 130 ;
1584     case  169961 : return 129 ;
1585     case  169956 : return 128 ;
1586     case  169926 : return 127 ;
1587     case  169924 : return 126 ;
1588     case  169923 : return 125 ;
1589     case  169922 : return 124 ;
1590     case  169919 : return 123 ;
1591     case  169918 : return 122 ;
1592     case  169914 : return 121 ;
1593     case  169859 : return 120 ;
1594     case  169858 : return 119 ;
1595     case  169855 : return 118 ;
1596     case  169846 : return 117 ;
1597     case  169838 : return 116 ;
1598     case  169837 : return 115 ;
1599     case  169835 : return 114 ;
1600     case  169683 : return 113 ;
1601     case  169628 : return 112 ;
1602     case  169591 : return 111 ;
1603     case  169590 : return 110 ;
1604     case  169588 : return 109 ;
1605     case  169587 : return 108 ;
1606     case  169586 : return 107 ;
1607     case  169584 : return 106 ;
1608     case  169557 : return 105 ;
1609     case  169555 : return 104 ;
1610     case  169554 : return 103 ;
1611     case  169553 : return 102 ;
1612     case  169550 : return 101 ;
1613     case  169515 : return 100 ;
1614     case  169512 : return 99 ;
1615     case  169506 : return 98 ;
1616     case  169504 : return 97 ;
1617     case  169498 : return 96 ;
1618     case  169475 : return 95 ;
1619     case  169420 : return 94 ;
1620     case  169419 : return 93 ;
1621     case  169418 : return 92 ;
1622     case  169417 : return 91 ;
1623     case  169415 : return 90 ;
1624     case  169411 : return 89 ;
1625     case  169238 : return 88 ;
1626     case  169236 : return 87 ;
1627     case  169167 : return 86 ;
1628     case  169160 : return 85 ;
1629     case  169156 : return 84 ;
1630     case  169148 : return 83 ;
1631     case  169145 : return 82 ;
1632     case  169144 : return 81 ;
1633     case  169143 : return 80 ;
1634     case  169138 : return 79 ;
1635     case  169099 : return 78 ;
1636     case  169094 : return 77 ;
1637     case  169091 : return 76 ;
1638     case  169045 : return 75 ;
1639     case  169044 : return 74 ;
1640     case  169040 : return 73 ;
1641     case  169035 : return 72 ;
1642     case  168992 : return 71 ;
1643     case  168988 : return 70 ;
1644     case  168984 : return 69 ;
1645     case  168826 : return 68 ;
1646     case  168777 : return 67 ;
1647     case  168514 : return 66 ;
1648     case  168512 : return 65 ;
1649     case  168511 : return 64 ;
1650     case  168467 : return 63 ;
1651     case  168464 : return 62 ;
1652     case  168461 : return 61 ;
1653     case  168460 : return 60 ;
1654     case  168458 : return 59 ;
1655     case  168362 : return 58 ;
1656     case  168361 : return 57 ;
1657     case  168356 : return 56 ;
1658     case  168342 : return 55 ;
1659     case  168341 : return 54 ;
1660     case  168325 : return 53 ;
1661     case  168322 : return 52 ;
1662     case  168318 : return 51 ;
1663     case  168311 : return 50 ;
1664     case  168310 : return 49 ;
1665     case  168213 : return 48 ;
1666     case  168212 : return 47 ;
1667     case  168208 : return 46 ;
1668     case  168207 : return 45 ;
1669     case  168206 : return 44 ;
1670     case  168205 : return 43 ;
1671     case  168204 : return 42 ;
1672     case  168203 : return 41 ;
1673     case  168181 : return 40 ;
1674     case  168177 : return 39 ;
1675     case  168175 : return 38 ;
1676     case  168173 : return 37 ;
1677     case  168172 : return 36 ;
1678     case  168171 : return 35 ;
1679     case  168115 : return 34 ;
1680     case  168108 : return 33 ;
1681     case  168107 : return 32 ;
1682     case  168105 : return 31 ;
1683     case  168104 : return 30 ;
1684     case  168103 : return 29 ;
1685     case  168076 : return 28 ;
1686     case  168069 : return 27 ;
1687     case  168068 : return 26 ;
1688     case  168066 : return 25 ;
1689     case  167988 : return 24 ;
1690     case  167987 : return 23 ;
1691     case  167986 : return 22 ;
1692     case  167985 : return 21 ;
1693     case  167921 : return 20 ;
1694     case  167920 : return 19 ;
1695     case  167915 : return 18 ;
1696     case  167909 : return 17 ;
1697     case  167903 : return 16 ;
1698     case  167902 : return 15 ;
1699     case  167818 : return 14 ;
1700     case  167814 : return 13 ;
1701     case  167813 : return 12 ;
1702     case  167808 : return 11 ;
1703     case  167807 : return 10 ;
1704     case  167806 : return 9 ;
1705     case  167713 : return 8 ;
1706     case  167712 : return 7 ;
1707     case  167711 : return 6 ;
1708     case  167706 : return 5 ;
1709     case  167693 : return 4 ;
1710     case  166532 : return 3 ;
1711     case  166530 : return 2 ;
1712     case  166529 : return 1 ;
1713
1714     default : return 199;
1715     }
1716   }
1717   if( kLHC10h == fPeriod ) {
1718     switch(run){
1719     case  139517 : return 137;
1720     case  139514 : return 136;
1721     case  139513 : return 135;
1722     case  139511 : return 134;
1723     case  139510 : return 133;
1724     case  139507 : return 132;
1725     case  139505 : return 131;
1726     case  139504 : return 130;
1727     case  139503 : return 129;
1728     case  139470 : return 128;
1729     case  139467 : return 127;
1730     case  139466 : return 126;
1731     case  139465 : return 125;
1732     case  139440 : return 124;
1733     case  139439 : return 123;
1734     case  139438 : return 122;
1735     case  139437 : return 121;
1736     case  139360 : return 120;
1737     case  139329 : return 119;
1738     case  139328 : return 118;
1739     case  139314 : return 117;
1740     case  139311 : return 116;
1741     case  139310 : return 115;
1742     case  139309 : return 114;
1743     case  139308 : return 113;
1744     case  139173 : return 112;
1745     case  139172 : return 111;
1746     case  139110 : return 110;
1747     case  139107 : return 109;
1748     case  139105 : return 108;
1749     case  139104 : return 107;
1750     case  139042 : return 106;
1751     case  139038 : return 105;
1752     case  139037 : return 104;
1753     case  139036 : return 103;
1754     case  139029 : return 102;
1755     case  139028 : return 101;
1756     case  138983 : return 100;
1757     case  138982 : return 99;
1758     case  138980 : return 98;
1759     case  138979 : return 97;
1760     case  138978 : return 96;
1761     case  138977 : return 95;
1762     case  138976 : return 94;
1763     case  138973 : return 93;
1764     case  138972 : return 92;
1765     case  138965 : return 91;
1766     case  138924 : return 90;
1767     case  138872 : return 89;
1768     case  138871 : return 88;
1769     case  138870 : return 87;
1770     case  138837 : return 86;
1771     case  138830 : return 85;
1772     case  138828 : return 84;
1773     case  138826 : return 83;
1774     case  138796 : return 82;
1775     case  138795 : return 81;
1776     case  138742 : return 80;
1777     case  138732 : return 79;
1778     case  138730 : return 78;
1779     case  138666 : return 77;
1780     case  138662 : return 76;
1781     case  138653 : return 75;
1782     case  138652 : return 74;
1783     case  138638 : return 73;
1784     case  138624 : return 72;
1785     case  138621 : return 71;
1786     case  138583 : return 70;
1787     case  138582 : return 69;
1788     case  138579 : return 68;
1789     case  138578 : return 67;
1790     case  138534 : return 66;
1791     case  138469 : return 65;
1792     case  138442 : return 64;
1793     case  138439 : return 63;
1794     case  138438 : return 62;
1795     case  138396 : return 61;
1796     case  138364 : return 60;
1797     case  138359 : return 59;
1798     case  138275 : return 58;
1799     case  138225 : return 57;
1800     case  138201 : return 56;
1801     case  138200 : return 55;
1802     case  138197 : return 54;
1803     case  138192 : return 53;
1804     case  138190 : return 52;
1805     case  138154 : return 51;
1806     case  138153 : return 50;
1807     case  138151 : return 49;
1808     case  138150 : return 48;
1809     case  138126 : return 47;
1810     case  138125 : return 46;
1811     case  137848 : return 45;
1812     case  137847 : return 44;
1813     case  137844 : return 43;
1814     case  137843 : return 42;
1815     case  137752 : return 41;
1816     case  137751 : return 40;
1817     case  137748 : return 39;
1818     case  137724 : return 38;
1819     case  137722 : return 37;
1820     case  137718 : return 36;
1821     case  137704 : return 35;
1822     case  137693 : return 34;
1823     case  137692 : return 33;
1824     case  137691 : return 32;
1825     case  137689 : return 31;
1826     case  137686 : return 30;
1827     case  137685 : return 29;
1828     case  137639 : return 28;
1829     case  137638 : return 27;
1830     case  137608 : return 26;
1831     case  137595 : return 25;
1832     case  137549 : return 24;
1833     case  137546 : return 23;
1834     case  137544 : return 22;
1835     case  137541 : return 21;
1836     case  137539 : return 20;
1837     case  137531 : return 19;
1838     case  137530 : return 18;
1839     case  137443 : return 17;
1840     case  137441 : return 16;
1841     case  137440 : return 15;
1842     case  137439 : return 14;
1843     case  137434 : return 13;
1844     case  137432 : return 12;
1845     case  137431 : return 11;
1846     case  137430 : return 10;
1847     case  137366 : return 9;
1848     case  137243 : return 8;
1849     case  137236 : return 7;
1850     case  137235 : return 6;
1851     case  137232 : return 5;
1852     case  137231 : return 4;
1853     case  137165 : return 3;
1854     case  137162 : return 2;
1855     case  137161 : return 1;
1856     default : return 199;
1857     }
1858   }
1859   if((fPeriod == kUndefinedPeriod) && (fDebug >= 1) ) {
1860     AliWarning("Period not defined");
1861   }
1862   return 1;
1863 }
1864 //_____________________________________________________________________________
1865 Bool_t AliAnalysisTaskPi0Flow::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1866
1867   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1868   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1869   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1870   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1871   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1872   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1873     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1874     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1875   return (R2<2.5*2.5) ;
1876
1877 }
1878 //_____________________________________________________________________________
1879 Bool_t AliAnalysisTaskPi0Flow::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1880
1881   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1882   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1883   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1884   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1885   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1886   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1887     0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1888     0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1889   return (R2<1.5*1.5) ;
1890
1891 }
1892 //____________________________________________________________________________
1893 TList* AliAnalysisTaskPi0Flow::GetCaloPhotonsPHOSList(UInt_t vtxBin, UInt_t centBin, UInt_t rpBin)
1894 {
1895   int offset = vtxBin * GetNumberOfCentralityBins() * fNEMRPBins
1896               + centBin * fNEMRPBins
1897               + rpBin;
1898   if( fCaloPhotonsPHOSLists->At(offset) ) { // list exists
1899     TList* list = dynamic_cast<TList*> (fCaloPhotonsPHOSLists->At(offset));
1900     if( ! list )
1901       AliError("object in fCaloPhotonsPHOSLists at %i did not cast");
1902     return list;
1903   }
1904   else {// no list for this bin has been created, yet
1905     TList* list = new TList();
1906     list->SetOwner();
1907     fCaloPhotonsPHOSLists->AddAt(list, offset);
1908     return list;
1909   }
1910 }
1911
1912 //____________________________________________________________________________
1913 Double_t AliAnalysisTaskPi0Flow::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1914   //Parameterization of LHC10h period
1915   //_true if neutral_
1916
1917   Double_t meanX=0;
1918   Double_t meanZ=0.;
1919   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1920                          6.58365e-01*5.91917e-01*5.91917e-01/((pt-9.61306e-01)*(pt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
1921   Double_t sz=TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(pt*pt+1.91456e-02*1.91456e-02)+1.60) ;
1922   Double_t mf = 0.; //
1923   if(fEventAOD) mf = fEventAOD->GetMagneticField(); //Positive for ++ and negative for --
1924   else if(fEventESD) mf = fEventESD->GetMagneticField(); //Positive for ++ and negative for --
1925
1926
1927   if(mf<0.){ //field --
1928     meanZ = -0.468318 ;
1929     if(charge>0)
1930       meanX=TMath::Min(7.3, 3.89994*1.20679*1.20679/(pt*pt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-pt*3.33650e+01)) ;
1931     else
1932       meanX=-TMath::Min(7.7,3.86040*0.912499*0.912499/(pt*pt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-pt*2.57070e+01)) ;
1933   }
1934   else{ //Field ++
1935     meanZ= -0.468318;
1936     if(charge>0)
1937       meanX=-TMath::Min(8.0,3.86040*1.31357*1.31357/(pt*pt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-pt*3.08451e+01)) ;
1938     else
1939       meanX= TMath::Min(6.85, 3.89994*1.16240*1.16240/(pt*pt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-pt*2.40913e+01)) ;
1940   }
1941
1942   Double_t rz=(dz-meanZ)/sz ;
1943   Double_t rx=(dx-meanX)/sx ;
1944   return TMath::Sqrt(rx*rx+rz*rz) ;
1945 }
1946 //____________________________________________________________________________
1947 void AliAnalysisTaskPi0Flow::SetFlatteningData(){
1948   //Read objects with flattening parameters
1949   AliOADBContainer flatContainer("phosFlat");
1950   flatContainer.InitFromFile(fEPcalibFileName.Data(),"phosFlat");
1951   TObjArray *maps = (TObjArray*)flatContainer.GetObject(fRunNumber,"phosFlat");
1952   if(!maps){
1953       AliError(Form("Can not read Flattening for run %d. \n From file >%s<\n",fRunNumber,fEPcalibFileName.Data())) ;
1954   }
1955   else{
1956     AliInfo(Form("Setting PHOS flattening with name %s \n",maps->GetName())) ;
1957     AliEPFlattener * h = (AliEPFlattener*)maps->At(0) ;
1958     if(fTPCFlat) delete fTPCFlat ;
1959     fTPCFlat = new AliEPFlattener() ;
1960     fTPCFlat = h ;
1961     h = (AliEPFlattener*)maps->At(1) ;
1962     if(fV0AFlat) delete fV0AFlat ;
1963     fV0AFlat = new AliEPFlattener() ;
1964     fV0AFlat = h ;
1965     h = (AliEPFlattener*)maps->At(2) ;
1966     if(fV0CFlat) delete fV0CFlat ;
1967     fV0CFlat = new AliEPFlattener() ;
1968     fV0CFlat = h ;
1969   }
1970
1971 }
1972  //____________________________________________________________________________
1973 Double_t  AliAnalysisTaskPi0Flow::ApplyFlattening(Double_t phi, Double_t c){
1974
1975   if(fTPCFlat)
1976     return fTPCFlat->MakeFlat(phi,c);
1977   return phi ;
1978
1979 }
1980 //____________________________________________________________________________
1981 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0A(Double_t phi, Double_t c){
1982
1983   if(fV0AFlat)
1984     return fV0AFlat->MakeFlat(phi,c);
1985   return phi ;
1986
1987 }
1988 //____________________________________________________________________________
1989 Double_t  AliAnalysisTaskPi0Flow::ApplyFlatteningV0C(Double_t phi, Double_t c){
1990
1991   if(fV0CFlat)
1992     return fV0CFlat->MakeFlat(phi,c);
1993   return phi ;
1994
1995 }
1996 //____________________________________________________________________________
1997 Double_t  AliAnalysisTaskPi0Flow::CoreEnergy(AliVCluster * clu, AliVCaloCells * cells)
1998 {
1999   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
2000
2001   //Can not use already calculated coordinates?
2002   //They have incidence correction...
2003   const Double_t distanceCut =3.5 ;
2004   const Double_t logWeight=4.5 ;
2005
2006   const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2007 // Calculates the center of gravity in the local PHOS-module coordinates
2008   Float_t wtot = 0;
2009   const Int_t mulDigit=clu->GetNCells() ;
2010   Double_t xc[mulDigit] ;
2011   Double_t zc[mulDigit] ;
2012   Double_t ei[mulDigit] ;
2013   Double_t x = 0 ;
2014   Double_t z = 0 ;
2015   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2016     Int_t relid[4] ;
2017     Float_t xi ;
2018     Float_t zi ;
2019     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
2020     fPHOSGeo->RelPosInModule(relid, xi, zi);
2021     xc[iDigit]=xi ;
2022     zc[iDigit]=zi ;
2023     ei[iDigit]=elist[iDigit]*cells->GetCellAmplitude(clu->GetCellsAbsId()[iDigit]);
2024     if( fDebug >= 3 )
2025       printf("%f ",ei[iDigit]);
2026     if (clu->E()>0 && ei[iDigit]>0) {
2027       Float_t w = TMath::Max( 0., logWeight + TMath::Log( ei[iDigit] / clu->E() ) ) ;
2028       x    += xc[iDigit] * w ;
2029       z    += zc[iDigit] * w ;
2030       wtot += w ;
2031     }
2032   }
2033   if (wtot>0) {
2034     x /= wtot ;
2035     z /= wtot ;
2036   }
2037   Double_t coreE=0. ;
2038   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
2039     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
2040     if(distance < distanceCut)
2041       coreE += ei[iDigit] ;
2042   }
2043   //Apply non-linearity correction
2044   return fNonLinCorr->Eval(coreE) ;
2045 }
2046 //____________________________________________________________________________
2047 Bool_t  AliAnalysisTaskPi0Flow::AreNeibors(Int_t id1,Int_t id2){
2048   // return true if absId are "Neighbors" (adjacent, including diagornaly,)
2049   // false if not.
2050
2051   Int_t relid1[4] ;
2052   fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
2053
2054   Int_t relid2[4] ;
2055   fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
2056
2057   // if inside the same PHOS module
2058   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) {
2059     const Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
2060     const Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
2061
2062     // and if diff in both direction is 1 or less
2063     if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
2064       return true; // are neighbors
2065   }
2066
2067   // else false
2068   return false;
2069 }
2070 //____________________________________________________________________________
2071 void  AliAnalysisTaskPi0Flow::Reclusterize(AliVCluster * clu){
2072   //Re-clusterize to make continues cluster
2073
2074   const Int_t oldMulDigit=clu->GetNCells() ;
2075   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2076   UShort_t * dlist = clu->GetCellsAbsId();
2077
2078   Int_t index[oldMulDigit] ;
2079   Bool_t used[oldMulDigit] ;
2080   for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
2081   Int_t inClu=0 ;
2082   Double_t eMax=0. ;
2083   //find maximum
2084   for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
2085     if(eMax<elist[iDigit]){
2086       eMax=elist[iDigit];
2087       index[0]=iDigit ;
2088       inClu=1 ;
2089     }
2090   }
2091   if(inClu==0){ //empty cluster
2092     return ;
2093   }
2094   used[index[0]]=kTRUE ; //mark as used
2095   for(Int_t i=0; i<inClu; i++){
2096     for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
2097        if(used[iDigit]) //already used
2098          continue ;
2099        if(AreNeibors(dlist[index[i]],dlist[iDigit])){
2100          index[inClu]= iDigit ;
2101          inClu++ ;
2102          used[iDigit]=kTRUE ;
2103        }
2104     }
2105   }
2106
2107   if(inClu==oldMulDigit) //no need to modify
2108     return ;
2109
2110   clu->SetNCells(inClu);
2111   //copy
2112   UShort_t tmpD[oldMulDigit] ;
2113   Double_t tmpE[oldMulDigit] ;
2114   for(Int_t i=0; i<oldMulDigit; i++){
2115     tmpD[i]=dlist[i] ;
2116     tmpE[i]=elist[i] ;
2117   }
2118   //change order of digits in list so that
2119   //first inClu cells were true ones
2120   for(Int_t i=0; i<inClu; i++){
2121     dlist[i]=tmpD[index[i]] ;
2122     elist[i]=tmpE[index[i]] ;
2123   }
2124
2125
2126 }
2127
2128 //_____________________________________________________________________________
2129 void AliAnalysisTaskPi0Flow::SetMisalignment(){
2130   // sets the misalignment vertex if ESD
2131   if( fEventESD ) {
2132     for(Int_t mod=0; mod<5; mod++) {
2133       const TGeoHMatrix* modMatrix = fEvent->GetPHOSMatrix(mod);
2134       if( ! modMatrix) {
2135         if( fDebug )
2136           AliInfo(Form("no PHOS Geometric Misalignment Matrix for module %d", mod));
2137         continue;
2138       }
2139       else {
2140         fPHOSGeo->SetMisalMatrix(modMatrix, mod);
2141         if( fDebug )
2142           AliInfo(Form("PHOS Geometric Misalignment Matrix set for module %d", mod));
2143       }
2144     }
2145   }
2146 }
2147
2148 //_____________________________________________________________________________
2149 void AliAnalysisTaskPi0Flow::SetV0Calibration(){
2150     // assigns: fMultV0, fV0Cpol, fV0Apol, fMeanQ, and fWidthQ
2151
2152     if ( ! fManualV0EPCalc ) {
2153       if( fDebug >=2 )
2154         AliInfo("Not setting V0Calibration, only needed for manual V0 EP Calculation");
2155       return;
2156     }
2157
2158     int runNumber = this->fRunNumber;
2159
2160     TString oadbfilename = "$ALICE_ROOT/OADB/PWGCF/VZERO/VZEROcalibEP.root";
2161     TFile *foadb = TFile::Open(oadbfilename.Data());
2162
2163     if(!foadb){
2164         AliError(Form("OADB file %s cannot be opened\n", oadbfilename.Data()));
2165         AliError("V0 Calibration not set !\n");
2166         return;
2167     }
2168
2169     AliOADBContainer *cont = (AliOADBContainer*) foadb->Get("hMultV0BefCorr");
2170     if(!cont){
2171         AliError("OADB object hMultV0BefCorr is not available in the file");
2172         AliError("V0 Calibration not set!\n");
2173         return;
2174     }
2175
2176     if(!(cont->GetObject(runNumber))){
2177         AliError(Form("OADB object hMultV0BefCorr is not available for run %i, trying 137366)",runNumber));
2178         runNumber = 137366;
2179     }
2180     if(!(cont->GetObject(runNumber))){
2181         AliError(Form("OADB object hMultV0BefCorr is not available for run %i ",runNumber));
2182         AliError("V0 Calibration not set!\n");
2183         return;
2184     }
2185
2186     if( fDebug )  AliInfo("Setting V0 calibration") ;
2187     fMultV0 = ((TH2F *) cont->GetObject(runNumber))->ProfileX();
2188
2189     TF1 *fpol0 = new TF1("fpol0","pol0");
2190     fMultV0->Fit(fpol0,"Q0","",0,31);
2191     fV0Cpol = fpol0->GetParameter(0);
2192     fMultV0->Fit(fpol0,"Q0","",32,64);
2193     fV0Apol = fpol0->GetParameter(0);
2194
2195     for(Int_t iside=0;iside<2;iside++){
2196         for(Int_t icoord=0;icoord<2;icoord++){
2197             for(Int_t i=0;i  < kNCenBins;i++){
2198                 char namecont[100];
2199                 if(iside==0 && icoord==0)
2200                     snprintf(namecont,100,"hQxc2_%i",i);
2201                 else if(iside==1 && icoord==0)
2202                     snprintf(namecont,100,"hQxa2_%i",i);
2203                 else if(iside==0 && icoord==1)
2204                     snprintf(namecont,100,"hQyc2_%i",i);
2205                 else if(iside==1 && icoord==1)
2206                     snprintf(namecont,100,"hQya2_%i",i);
2207
2208                 cont = (AliOADBContainer*) foadb->Get(namecont);
2209                 if(!cont){
2210                     AliError(Form("OADB object %s is not available in the file %s", namecont, oadbfilename.Data()));
2211                     AliError("V0 Calibration not fully set!\n");
2212                     return;
2213                 }
2214
2215                 if(!(cont->GetObject(runNumber))){
2216                     AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2217                     runNumber = 137366;
2218                 }
2219                 if(!(cont->GetObject(runNumber))){
2220                   AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2221                   AliError("V0 Calibration not fully set!\n");
2222                   return;
2223                 }
2224                 fMeanQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetMean();
2225                 fWidthQ[i][iside][icoord] = ((TH1F *) cont->GetObject(runNumber))->GetRMS();
2226
2227                 //for v3
2228 //              if(iside==0 && icoord==0)
2229 //                  snprintf(namecont,100,"hQxc3_%i",i);
2230 //              else if(iside==1 && icoord==0)
2231 //                  snprintf(namecont,100,"hQxa3_%i",i);
2232 //              else if(iside==0 && icoord==1)
2233 //                  snprintf(namecont,100,"hQyc3_%i",i);
2234 //              else if(iside==1 && icoord==1)
2235 //                  snprintf(namecont,100,"hQya3_%i",i);
2236 //
2237 //              cont = (AliOADBContainer*) foadb->Get(namecont);
2238 //              if(!cont){
2239 //                  AliError(Form("OADB object %s is not available in the file",namecont));
2240 //                  AliError("V0 Calibration not fully set!\n");
2241 //                  return;
2242 //              }
2243 //
2244 //              if(!(cont->GetObject(runNumber))){
2245 //                  AliError(Form("OADB object %s is not available for run %i, trying run 137366",namecont,runNumber));
2246 //                  runNumber = 137366;
2247 //              }
2248 //              if(!(cont->GetObject(runNumber))){
2249 //                AliError(Form("OADB object %s is not available for run %i",namecont,runNumber));
2250 //                AliError("V0 Calibration not fully set!\n");
2251 //                return;
2252 //              }
2253 //              fMeanQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetMean();
2254 //              fWidthQv3[i][iside][icoord] = ((TH1F *) cont->GetObject(run))->GetRMS();
2255
2256             }
2257         }
2258     }
2259
2260     delete fpol0; fpol0=0;
2261 }
2262
2263 //_____________________________________________________________________________
2264 void AliAnalysisTaskPi0Flow::SetESDTrackCuts()
2265 {
2266   if( fEventESD ) {
2267     // Create ESD track cut
2268     fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts() ;
2269     //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010();
2270     fESDtrackCuts->SetRequireTPCRefit(kTRUE);
2271   }
2272 }
2273
2274 //_____________________________________________________________________________
2275 void AliAnalysisTaskPi0Flow::SetGeometry()
2276 {
2277   // Initialize the PHOS geometry
2278   if( kLHC10h == fPeriod && fEventESD ) {
2279     TGeoManager::Import("geometry.root"); //TODO: should perhaps not be done
2280     fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
2281     if( ! fPHOSGeo )
2282       AliError("geometry (fPHOSGeo) not initialised");
2283   }
2284
2285   //Init geometry
2286   if(!fPHOSGeo){
2287     AliOADBContainer geomContainer("phosGeo");
2288     geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
2289     TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
2290     fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP") ;
2291     for(Int_t mod=0; mod<5; mod++) {
2292       if(!matrixes->At(mod)) {
2293         if( fDebug )
2294           AliInfo(Form("No PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2295         continue;
2296       }
2297       else {
2298         fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod) ;
2299         if( fDebug >1 )
2300           AliInfo(Form("Adding PHOS Matrix for mod:%d, geo=%p\n", mod, fPHOSGeo));
2301       }
2302     }
2303   }
2304 }
2305
2306 //_____________________________________________________________________________
2307 void AliAnalysisTaskPi0Flow::SetPHOSCalibData()
2308 {
2309   if( fPHOSCalibData )
2310     delete fPHOSCalibData;
2311   fPHOSCalibData = 0;
2312
2313   // Calibration only needed for ESD
2314   if( fEventESD /*&& */ ) {
2315     if( kLHC10h == fPeriod && fEventESD ) {
2316       //We have to apply re-calibration for pass1 LCH10h
2317       // Initialize decalibration factors in the form of the OCDB object
2318       AliCDBManager * man = AliCDBManager::Instance();
2319       man->SetRun(140000) ; //TODO; revise, this should probably not b done.
2320       man->SetDefaultStorage("local://OCDB");
2321     }
2322     fPHOSCalibData = new AliPHOSCalibData();
2323   }
2324 }
2325
2326 //_____________________________________________________________________________
2327 void AliAnalysisTaskPi0Flow::SetVertex()
2328 {
2329   const AliVVertex *primaryVertex = fEvent->GetPrimaryVertex();
2330   if( primaryVertex ) {
2331     fVertex[0] = primaryVertex->GetX();
2332     fVertex[1] = primaryVertex->GetY();
2333     fVertex[2] = primaryVertex->GetZ();
2334   }
2335   else {
2336     AliError("Event has 0x0 Primary Vertex, defaulting to origo");
2337     fVertex[0] = 0;
2338     fVertex[1] = 0;
2339     fVertex[2] = 0;
2340   }
2341   fVertexVector = TVector3(fVertex);
2342   FillHistogram("hZvertex", fVertexVector.z(), fInternalRunNumber-0.5);
2343
2344   if( fDebug >= 2 )
2345     AliInfo(Form("Vertex is set to (%.1f,%.1f,%.1f)", fVertex[0], fVertex[1], fVertex[2]));
2346
2347   fVtxBin=0 ;// No support for vtx binning implemented.
2348 }
2349
2350 //_____________________________________________________________________________
2351 Bool_t AliAnalysisTaskPi0Flow::RejectEventVertex()
2352 {
2353   if( ! fEvent->GetPrimaryVertex() )
2354     return true; // reject
2355   LogSelection(1, fInternalRunNumber);
2356
2357   if ( TMath::Abs(fVertexVector.z()) > fMaxAbsVertexZ )
2358     return true; // reject
2359   LogSelection(2, fInternalRunNumber);
2360
2361   return false; // accept event.
2362 }
2363
2364 //_____________________________________________________________________________
2365 void AliAnalysisTaskPi0Flow::SetCentrality()
2366 {
2367   AliCentrality *centrality = fEvent->GetCentrality();
2368   if( centrality )
2369     fCentralityV0M=centrality->GetCentralityPercentile("V0M");
2370   else {
2371     AliError("Event has 0x0 centrality");
2372     fCentralityV0M = -1.;
2373   }
2374   FillHistogram("hCentrality",fCentralityV0M,fInternalRunNumber-0.5) ;
2375
2376   fCentBin = GetCentralityBin(fCentralityV0M);
2377
2378   if ( fDebug >= 2 )
2379     AliInfo(Form("Centrality (bin) is: %f (%d)", fCentralityV0M, fCentBin));
2380 }
2381
2382 //_____________________________________________________________________________
2383 Bool_t AliAnalysisTaskPi0Flow::RejectCentrality()
2384 {
2385   if( ! fEvent->GetCentrality() )
2386     return true; // reject
2387   LogSelection(3, fInternalRunNumber);
2388
2389 //   if( fCentralityV0M <= 0. || fCentralityV0M>80. )
2390 //     return true; // reject
2391
2392   int lastBinUpperIndex = fCentEdges.GetSize() -1;
2393   if( fCentralityV0M > fCentEdges[lastBinUpperIndex] ) {
2394     if( fDebug )
2395       AliInfo("Rejecting due to centrality outside of binning.");
2396     return true; // reject
2397   }
2398   LogSelection(4, fInternalRunNumber);
2399
2400   if( fCentralityV0M < fCentEdges[0] ) {
2401     if( fDebug )
2402       AliInfo("Rejecting due to centrality outside of binning.");
2403     return true; // reject
2404   }
2405   LogSelection(5, fInternalRunNumber);
2406
2407   return false;
2408 }
2409
2410
2411 //_____________________________________________________________________________
2412 void AliAnalysisTaskPi0Flow::EvalReactionPlane()
2413 {
2414   // assigns: fHaveTPCRP and fRP
2415   // also does a few histogram fills
2416
2417   AliEventplane *eventPlane = fEvent->GetEventplane();
2418   if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2419
2420   Double_t reactionPlaneQ = eventPlane->GetEventplane("Q");
2421   FillHistogram("phiRP",reactionPlaneQ,fCentralityV0M) ;
2422
2423   if(reactionPlaneQ==999 || reactionPlaneQ < 0.){ //reaction plain was not defined
2424     if( fDebug ) AliInfo(Form("No Q Reaction Plane, value is %f", reactionPlaneQ));
2425     fHaveTPCRP = kFALSE;
2426   }
2427   else{
2428     if( fDebug >= 2 ) AliInfo(Form("Q Reaction Plane is %f", reactionPlaneQ));
2429     fHaveTPCRP = kTRUE;
2430   }
2431
2432   if(fHaveTPCRP){
2433     fRP = ApplyFlattening(reactionPlaneQ, fCentralityV0M) ;
2434
2435     while(fRP<0)  fRP+=TMath::Pi();
2436     while(fRP>TMath::Pi())  fRP-=TMath::Pi();
2437     FillHistogram("phiRPflat",fRP,fCentralityV0M) ;
2438     Double_t dPsi = eventPlane->GetQsubRes() ;
2439     FillHistogram("cos2AC",TMath::Cos(2.*dPsi),fCentralityV0M) ;
2440   }
2441   else
2442     fRP=0.;
2443 }
2444
2445
2446 //____________________________________________________________________________
2447 void  AliAnalysisTaskPi0Flow::EvalV0ReactionPlane(){
2448   // set: fRPV0A and fRPV0C
2449
2450   // Do Manual V0 EP Calculation
2451   if ( fManualV0EPCalc )
2452     {
2453       //VZERO data
2454       AliVVZERO* v0 = fEvent->GetVZEROData();
2455
2456       //reset Q vector info
2457       Double_t Qxa2 = 0, Qya2 = 0;
2458       Double_t Qxc2 = 0, Qyc2 = 0;
2459
2460       for (Int_t iv0 = 0; iv0 < 64; iv0++) {
2461         Double_t phiV0 = TMath::PiOver4()*(0.5 + iv0 % 8);
2462         Float_t multv0 = v0->GetMultiplicity(iv0);
2463         if (iv0 < 32){ // V0C
2464           Qxc2 += TMath::Cos(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2465           Qyc2 += TMath::Sin(2*phiV0) * multv0*fV0Cpol/fMultV0->GetBinContent(iv0+1);
2466         } else {       // V0A
2467           Qxa2 += TMath::Cos(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2468           Qya2 += TMath::Sin(2*phiV0) * multv0*fV0Apol/fMultV0->GetBinContent(iv0+1);
2469         }
2470       }
2471
2472       Int_t iC = -1;
2473       // centrality bins
2474       if(fCentralityV0M < 5) iC = 0;
2475       else if(fCentralityV0M < 10) iC = 1;
2476       else if(fCentralityV0M < 20) iC = 2;
2477       else if(fCentralityV0M < 30) iC = 3;
2478       else if(fCentralityV0M < 40) iC = 4;
2479       else if(fCentralityV0M < 50) iC = 5;
2480       else if(fCentralityV0M < 60) iC = 6;
2481       else if(fCentralityV0M < 70) iC = 7;
2482       else iC = 8;
2483
2484       //grab for each centrality the proper histo with the Qx and Qy to do the recentering
2485       Double_t Qxamean2 = fMeanQ[iC][1][0];
2486       Double_t Qxarms2  = fWidthQ[iC][1][0];
2487       Double_t Qyamean2 = fMeanQ[iC][1][1];
2488       Double_t Qyarms2  = fWidthQ[iC][1][1];
2489
2490       Double_t Qxcmean2 = fMeanQ[iC][0][0];
2491       Double_t Qxcrms2  = fWidthQ[iC][0][0];
2492       Double_t Qycmean2 = fMeanQ[iC][0][1];
2493       Double_t Qycrms2  = fWidthQ[iC][0][1];
2494
2495       Double_t QxaCor2 = (Qxa2 - Qxamean2)/Qxarms2;
2496       Double_t QyaCor2 = (Qya2 - Qyamean2)/Qyarms2;
2497       Double_t QxcCor2 = (Qxc2 - Qxcmean2)/Qxcrms2;
2498       Double_t QycCor2 = (Qyc2 - Qycmean2)/Qycrms2;
2499
2500       fRPV0A = TMath::ATan2(QyaCor2, QxaCor2)/2.;
2501       fRPV0C = TMath::ATan2(QycCor2, QxcCor2)/2.;
2502     }
2503   else // Use Official V0 EP Calculation.
2504     {
2505       AliEventplane *eventPlane = fEvent->GetEventplane();
2506       if( ! eventPlane ) { AliError("Event has no event plane"); return; }
2507       fRPV0A = eventPlane->GetEventplane("V0A", fEvent);
2508       fRPV0C = eventPlane->GetEventplane("V0C", fEvent);
2509     }
2510
2511   // Check that the A&C RP are within allowed range.
2512   if( fDebug >= 3 && (fRPV0A<0 || fRPV0A>TMath::Pi() ) )
2513     AliInfo(Form("RPV0A outside of permited range [0,pi]: %f, correcting", fRPV0A));
2514   if( fDebug >= 3 && (fRPV0C<0 || fRPV0C>TMath::Pi() ) )
2515     AliInfo(Form("RPV0C outside of permited range [0,pi]: %f, correcting", fRPV0C));
2516   while (fRPV0A<0          ) fRPV0A+=TMath::Pi() ;
2517   while (fRPV0A>TMath::Pi()) fRPV0A-=TMath::Pi() ;
2518   while (fRPV0C<0          ) fRPV0C+=TMath::Pi() ;
2519   while (fRPV0C>TMath::Pi()) fRPV0C-=TMath::Pi() ;
2520
2521   // Reaction plane histograms before flattening
2522   if( fDebug >= 2 )
2523     AliInfo(Form("V0 Reaction Plane before flattening: A side: %f, C side: %f", fRPV0A, fRPV0C));
2524
2525   FillHistogram("phiRPV0A" ,fRPV0A,fCentralityV0M);
2526   FillHistogram("phiRPV0C" ,fRPV0C,fCentralityV0M);
2527   FillHistogram("phiRPV0AC",fRPV0A,fRPV0C,fCentralityV0M) ;
2528
2529   // Flattening
2530   fRPV0A=ApplyFlatteningV0A(fRPV0A,fCentralityV0M) ;
2531   while (fRPV0A<0          ) fRPV0A+=TMath::Pi() ;
2532   while (fRPV0A>TMath::Pi()) fRPV0A-=TMath::Pi() ;
2533
2534   fRPV0C=ApplyFlatteningV0C(fRPV0C,fCentralityV0M) ;
2535   while (fRPV0C<0          ) fRPV0C+=TMath::Pi() ;
2536   while (fRPV0C>TMath::Pi()) fRPV0C-=TMath::Pi() ;
2537
2538   if( fDebug >= 2 )
2539     AliInfo(Form("V0 Reaction Plane after  flattening: A side: %f, C side: %f", fRPV0A, fRPV0C));
2540
2541   FillHistogram("phiRPV0Aflat",fRPV0A,fCentralityV0M) ;
2542   FillHistogram("cos2V0AC",TMath::Cos(2.*(fRPV0A-fRPV0C)),fCentralityV0M) ;
2543   if(fHaveTPCRP){
2544     FillHistogram("phiRPV0ATPC",fRP,fRPV0A,fCentralityV0M) ;
2545     FillHistogram("cos2V0ATPC",TMath::Cos(2.*(fRP-fRPV0A)),fCentralityV0M) ;
2546   }
2547
2548   FillHistogram("phiRPV0Cflat",fRPV0C,fCentralityV0M) ;
2549   if(fHaveTPCRP){
2550     FillHistogram("phiRPV0CTPC",fRP,fRPV0C,fCentralityV0M) ;
2551     FillHistogram("cos2V0CTPC",TMath::Cos(2.*(fRP-fRPV0C)),fCentralityV0M) ;
2552   }
2553 }
2554 //____________________________________________________________________________
2555 void  AliAnalysisTaskPi0Flow::EvalCoreLambdas(AliVCluster * clu, AliVCaloCells * cells,Double_t &m02, Double_t &m20){
2556   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
2557
2558   const Double_t rCut=4.5 ;
2559
2560   const Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
2561 // Calculates the center of gravity in the local PHOS-module coordinates
2562   Float_t wtot = 0;
2563   const Int_t mulDigit=clu->GetNCells() ;
2564   Double_t xc[mulDigit] ;
2565   Double_t zc[mulDigit] ;
2566   Double_t wi[mulDigit] ;
2567   Double_t x = 0 ;
2568   Double_t z = 0 ;
2569   const Double_t logWeight=4.5 ;
2570   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2571     Int_t relid[4] ;
2572     Float_t xi=0. ;
2573     Float_t zi=0. ;
2574     Int_t absId = clu->GetCellAbsId(iDigit) ;
2575     fPHOSGeo->AbsToRelNumbering(absId, relid) ;
2576     fPHOSGeo->RelPosInModule(relid, xi, zi);
2577     xc[iDigit]=xi ;
2578     zc[iDigit]=zi ;
2579     Double_t ei = elist[iDigit]*cells->GetCellAmplitude(absId) ;
2580     wi[iDigit]=0. ;
2581     if (clu->E()>0 && ei>0) {
2582       wi[iDigit] = TMath::Max( 0., logWeight + TMath::Log( ei / clu->E() ) ) ;
2583       Double_t w=wi[iDigit];
2584       x    += xc[iDigit] * w ;
2585       z    += zc[iDigit] * w ;
2586       wtot += w ;
2587     }
2588   }
2589   if (wtot>0) {
2590     x /= wtot ;
2591     z /= wtot ;
2592   }
2593
2594   wtot = 0. ;
2595   Double_t dxx  = 0.;
2596   Double_t dzz  = 0.;
2597   Double_t dxz  = 0.;
2598   Double_t xCut = 0. ;
2599   Double_t zCut = 0. ;
2600   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
2601     Double_t w=wi[iDigit];
2602     if (w>0.) {
2603         Double_t xi= xc[iDigit] ;
2604         Double_t zi= zc[iDigit] ;
2605         if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
2606           xCut += w * xi ;
2607           zCut += w * zi ;
2608           dxx  += w * xi * xi ;
2609           dzz  += w * zi * zi ;
2610           dxz  += w * xi * zi ;
2611           wtot += w ;
2612         }
2613     }
2614
2615   }
2616   if (wtot>0) {
2617     xCut/= wtot ;
2618     zCut/= wtot ;
2619     dxx /= wtot ;
2620     dzz /= wtot ;
2621     dxz /= wtot ;
2622     dxx -= xCut * xCut ;
2623     dzz -= zCut * zCut ;
2624     dxz -= xCut * zCut ;
2625
2626     m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
2627     m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
2628   }
2629   else {
2630     m20=m02=0.;
2631   }
2632
2633 }
2634 //____________________________________________________________________________
2635 Bool_t AliAnalysisTaskPi0Flow::TestCoreLambda(Double_t pt,Double_t l1,Double_t l2){
2636   //Evaluates if lambdas correspond to photon cluster
2637   //Tuned using pp date
2638   //For core radius R=4.5
2639   Double_t   l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
2640   Double_t   l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
2641   Double_t   l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
2642   Double_t   l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
2643   Double_t   c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
2644
2645   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
2646               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
2647               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
2648   return (R2<2.5*2.5) ;
2649 }