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