PWGGA/PHOSTasks/PHOS_PbPb_MC/AliPHOSHijingEfficiency.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TObjArray.h"
   4 #include "TF1.h"
   5 #include "TFile.h"
   6 #include "TH1F.h"
   7 #include "TH2F.h"
   8 #include "TH2I.h"
   9 #include "TH3F.h"
  10 #include "TParticle.h"
  11 #include "TCanvas.h"
  12 #include "TStyle.h"
  13 #include "TRandom.h"
  14 #include "TROOT.h"
  15 #include "THashList.h"
  16
  17 #include "AliAnalysisManager.h"
  18 #include "AliMCEventHandler.h"
  19 #include "AliMCEvent.h"
  20 #include "AliStack.h"
  21 #include "AliAnalysisTaskSE.h"
  22 #include "AliPHOSHijingEfficiency.h"
  23 #include "AliCaloPhoton.h"
  24 #include "AliPHOSGeometry.h"
  25 #include "AliPHOSEsdCluster.h"
  26 #include "AliPHOSCalibData.h"
  27 #include "AliESDEvent.h"
  28 #include "AliESDCaloCells.h"
  29 #include "AliESDVertex.h"
  30 #include "AliESDtrackCuts.h"
  31 #include "AliLog.h"
  32 #include "AliPID.h"
  33 #include "AliCDBManager.h"
  34 #include "AliCentrality.h"
  35 #include "AliESDtrackCuts.h"
  36 #include "AliEventplane.h"
  37 #include "TProfile.h"
  38 #include "AliOADBContainer.h"
  39
  40 // Calculates pi0 and photon efficiencies using Hijing events
  41 // Also calculates pi0 contamination: non-vertex pi0s
  42 // And photon contamination:
  43 //    - off-vertex photons
  44 //    - misidentified photons
  45 // Authors: Dmitri Peressounko
  46 // Date   : 05.02.2013
  47
  48 ClassImp(AliPHOSHijingEfficiency)
  49
  50 //________________________________________________________________________
  51 Double_t rnlin(Double_t x)
  52 {
  53   //This is standard PHOS non-linearity used in reconstruction
  54   //and re-calibration by 2.2%
  55   return 1.022*(0.0241+1.0504*x+0.000249*x*x) ;
  56
  57 }
  58
  59 //________________________________________________________________________
  60 AliPHOSHijingEfficiency::AliPHOSHijingEfficiency(const char *name)
  61 : AliAnalysisTaskSE(name),
  62   fStack(0x0),
  63   fOutputContainer(0x0),
  64   fPHOSEvent(0x0),
  65   fPHOSCalibData(0x0),
  66   fRP(0.),fRPV0A(0.),fRPV0C(0.),fHaveTPCRP(0),
  67   fRunNumber(0),
  68   fCentrality(0.),
  69   fCenBin(0),
  70   fPHOSGeo(0),
  71   fEventCounter(0)
  72 {
  73   // Constructor
  74   for(Int_t i=0;i<1;i++){
  75     for(Int_t j=0;j<10;j++)
  76       for(Int_t k=0;k<11;k++)
  77         fPHOSEvents[i][j][k]=0 ;
  78   }
  79
  80   // Output slots #0 write into a TH1 container
  81   DefineOutput(1,TList::Class());
  82
  83   // Set bad channel map
  84   char key[55] ;
  85   for(Int_t i=0; i<6; i++){
  86     sprintf(key,"PHOS_BadMap_mod%d",i) ;
  87     fPHOSBadMap[i]=new TH2I(key,"Bad Modules map",64,0.,64.,56,0.,56.) ;
  88   }
  89
  90
  91 }
  92
  93 //________________________________________________________________________
  94 void AliPHOSHijingEfficiency::UserCreateOutputObjects()
  95 {
  96
  97
  98   // Create histograms
  99   // Called once
 100   const Int_t nRuns=200 ;
 101
 102   // ESD histograms
 103   if(fOutputContainer != NULL){
 104     delete fOutputContainer;
 105   }
 106   fOutputContainer = new THashList();
 107   fOutputContainer->SetOwner(kTRUE);
 108   PostData(1, fOutputContainer);
 109
 110   //========QA histograms=======
 111
 112   //Event selection
 113   fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 10,0.,10.,nRuns,0.,float(nRuns))) ;
 114   fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 10,0.,10.)) ;
 115
 116   //vertex distribution
 117   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
 118
 119   //Centrality
 120   fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
 121   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
 122   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
 123   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
 124   fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,20.,100,0.,100.)) ;
 125   fOutputContainer->Add(new TH2F("hCluEvsCluM","ClusterMult vs E",200,0.,20.,100,0.,20.)) ;
 126
 127   fOutputContainer->Add(new TH3F("hCPVr","CPV radius",100,0.,20.,100,0.,20.,10,0.,100.));
 128
 129   //Bad Map
 130   fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 131   fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 132   fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 133
 134   fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 135   fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 136   fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 137
 138   fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 139   fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 140   fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 141
 142   fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 143   fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 144   fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 145
 146   //Single photon and pi0 spectrum
 147   Int_t nPtPhot = 300 ;
 148   Double_t ptPhotMax = 30 ;
 149   Int_t nM       = 500;
 150   Double_t mMin  = 0.0;
 151   Double_t mMax  = 1.0;
 152
 153   //PHOS calibration QA
 154   fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 155   fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 156   fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 157   fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 158   fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 159   fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 160
 161   char key[55] ;
 162   for(Int_t cent=0; cent<6; cent++){
 163
 164     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 165     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 166     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 167
 168     fOutputContainer->Add(new TH1F(Form("hPhotAll_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 169     fOutputContainer->Add(new TH1F(Form("hPhotAllcore_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 170     fOutputContainer->Add(new TH1F(Form("hPhotAllwou_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 171     fOutputContainer->Add(new TH1F(Form("hPhotDisp_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 172     fOutputContainer->Add(new TH1F(Form("hPhotDisp2_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 173     fOutputContainer->Add(new TH1F(Form("hPhotDispcore_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 174     fOutputContainer->Add(new TH1F(Form("hPhotDisp2core_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 175     fOutputContainer->Add(new TH1F(Form("hPhotDispwou_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 176     fOutputContainer->Add(new TH1F(Form("hPhotCPV_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 177     fOutputContainer->Add(new TH1F(Form("hPhotCPVcore_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 178     fOutputContainer->Add(new TH1F(Form("hPhotCPV2_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 179     fOutputContainer->Add(new TH1F(Form("hPhotBoth_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
 180     fOutputContainer->Add(new TH1F(Form("hPhotBothcore_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
 181     fOutputContainer->Add(new TH1F(Form("hPhotBoth2_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
 182     fOutputContainer->Add(new TH1F(Form("hPhotBoth2core_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
 183
 184
 185     fOutputContainer->Add(new TH2F(Form("hPi0All_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 186     fOutputContainer->Add(new TH2F(Form("hPi0Allcore_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 187     fOutputContainer->Add(new TH2F(Form("hPi0Allwou_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 188     fOutputContainer->Add(new TH2F(Form("hPi0Disp_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 189     fOutputContainer->Add(new TH2F(Form("hPi0Disp2_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 190     fOutputContainer->Add(new TH2F(Form("hPi0Dispcore_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 191     fOutputContainer->Add(new TH2F(Form("hPi0Disp2core_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 192     fOutputContainer->Add(new TH2F(Form("hPi0Dispwou_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 193     fOutputContainer->Add(new TH2F(Form("hPi0CPV_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 194     fOutputContainer->Add(new TH2F(Form("hPi0CPVcore_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 195     fOutputContainer->Add(new TH2F(Form("hPi0CPV2_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 196     fOutputContainer->Add(new TH2F(Form("hPi0Both_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 197     fOutputContainer->Add(new TH2F(Form("hPi0Bothcore_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 198     fOutputContainer->Add(new TH2F(Form("hPi0Both2_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 199     fOutputContainer->Add(new TH2F(Form("hPi0Both2core_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 200
 201
 202     snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
 203     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 204     snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
 205     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 206     snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
 207     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 208     snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
 209     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 210     snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
 211     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 212     snprintf(key,55,"hPi0Both2_a07_cen%d",cent) ;
 213     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 214
 215     snprintf(key,55,"hSingleAll_cen%d",cent) ;
 216     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 217     snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
 218     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 219     snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
 220     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 221     snprintf(key,55,"hSingleDisp_cen%d",cent) ;
 222     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 223     snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
 224     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 225     snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
 226     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 227     snprintf(key,55,"hSingleDisp2core_cen%d",cent) ;
 228     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 229     snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
 230     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 231     snprintf(key,55,"hSingleCPV_cen%d",cent) ;
 232     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 233     snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
 234     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 235     snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
 236     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 237     snprintf(key,55,"hSingleBoth_cen%d",cent) ;
 238     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 239     snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
 240     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 241     snprintf(key,55,"hSingleBoth2_cen%d",cent) ;
 242     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 243     snprintf(key,55,"hSingleBoth2core_cen%d",cent) ;
 244     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 245
 246
 247     snprintf(key,55,"hMiPi0All_cen%d",cent) ;
 248     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 249     snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
 250     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 251     snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
 252     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 253     snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
 254     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 255     snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
 256     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 257     snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
 258     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 259     snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
 260     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 261     snprintf(key,55,"hMiPi0Disp2core_cen%d",cent) ;
 262     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 263     snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
 264     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 265     snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
 266     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 267     snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
 268     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 269     snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
 270     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 271     snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
 272     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 273     snprintf(key,55,"hMiPi0Both2_cen%d",cent) ;
 274     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 275     snprintf(key,55,"hMiPi0Both2core_cen%d",cent) ;
 276     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 277
 278     snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
 279     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 280     snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
 281     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 282     snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
 283     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 284     snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
 285     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 286     snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
 287     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 288     snprintf(key,55,"hMiPi0Both2_a07_cen%d",cent) ;
 289     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 290
 291     snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
 292     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 293     snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
 294     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 295     snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
 296     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 297     snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
 298     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 299     snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
 300     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 301     snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
 302     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 303     snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
 304     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 305     snprintf(key,55,"hMiSingleDisp2core_cen%d",cent) ;
 306     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 307     snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
 308     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 309     snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
 310     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 311     snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
 312     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 313     snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
 314     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 315     snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
 316     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 317     snprintf(key,55,"hMiSingleBoth2_cen%d",cent) ;
 318     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 319     snprintf(key,55,"hMiSingleBoth2core_cen%d",cent) ;
 320     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 321   }
 322
 323
 324  //MC
 325   for(Int_t cent=0; cent<6; cent++){
 326     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 327     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 328     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 329     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 330     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 331     fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
 332     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 333     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 334     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 335     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 336     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 337     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 338     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 339     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 340     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 341     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 342     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 343     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 344     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 345     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 346     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 347     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 348     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 349     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 350   }
 351   fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 352   fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 353   fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 354
 355
 356   Int_t nPt      = 200;
 357   Double_t ptMin = 0;
 358   Double_t ptMax = 20;
 359   fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
 360   fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
 361   fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
 362   fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
 363   fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
 364   fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
 365 //  fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
 366 //  fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
 367
 368   for(Int_t cen=0; cen<6; cen++){
 369     fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 370     fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 371     fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 372     fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 373     fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 374     fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 375     fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 376     fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 377     fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 378     fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 379     fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 380     fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 381     fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 382
 383     //pairs from common parents
 384     fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 385     fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 386     fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 387     fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 388     fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 389     fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 390     fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 391
 392     //common parent - pi0
 393     fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 394     fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 395     fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 396     fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 397     fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 398     fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 399     fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 400     fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 401     fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 402
 403   }
 404
 405
 406   //Photon contaminations
 407   fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 408   fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 409   fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 410   fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 411   fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
 412
 413    const Int_t nTypes=24 ;
 414    char partTypes[nTypes][55] ;
 415    snprintf(partTypes[0],55,"hGammaNoPrim") ; //
 416    snprintf(partTypes[1],55,"hGammaPhot") ; //
 417    snprintf(partTypes[2],55,"hGammaEl") ; //
 418    snprintf(partTypes[3],55,"hGammaPi0") ; //
 419    snprintf(partTypes[4],55,"hGammaEta") ; //
 420    snprintf(partTypes[5],55,"hhGammaOmega") ; //
 421    snprintf(partTypes[6],55,"hGammaPipm") ; //
 422    snprintf(partTypes[7],55,"hGammaP") ; //
 423    snprintf(partTypes[8],55,"hGammaPbar") ; //
 424    snprintf(partTypes[9],55,"hGammaN") ; //
 425    snprintf(partTypes[10],55,"hGammaNbar") ; //
 426    snprintf(partTypes[11],55,"hGammaK0S") ; //
 427    snprintf(partTypes[12],55,"hGammaK0L") ; //
 428    snprintf(partTypes[13],55,"hGammaKpm") ; //
 429    snprintf(partTypes[14],55,"hGammaKstar") ; //
 430    snprintf(partTypes[15],55,"hGammaDelta") ; //
 431    snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
 432    snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
 433    snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
 434    snprintf(partTypes[19],55,"hGammaPipmEl") ; //
 435    snprintf(partTypes[20],55,"hGammaPipmOther") ; //
 436    snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
 437    snprintf(partTypes[22],55,"hGammaPipmp") ; //
 438    snprintf(partTypes[23],55,"hGammaPipmn") ; //
 439
 440    const Int_t nPID=12 ;
 441    char cPID[25][12] ;
 442    snprintf(cPID[0],25,"All") ;
 443    snprintf(cPID[1],25,"Allcore") ;
 444    snprintf(cPID[2],25,"CPV") ;
 445    snprintf(cPID[3],25,"CPVcore") ;
 446    snprintf(cPID[4],25,"CPV2") ;
 447    snprintf(cPID[5],25,"CPV2core") ;
 448    snprintf(cPID[6],25,"Disp") ;
 449    snprintf(cPID[7],25,"Dispcore") ;
 450    snprintf(cPID[8],25,"Disp2") ;
 451    snprintf(cPID[9],25,"Disp2core") ;
 452    snprintf(cPID[10],25,"Both") ;
 453    snprintf(cPID[11],25,"Bothcore") ;
 454
 455    for(Int_t itype=0; itype<nTypes; itype++){
 456      for(Int_t iPID=0; iPID<nPID; iPID++){
 457        for(Int_t cen=0; cen<5; cen++){
 458          fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
 459        }
 460      }
 461    }
 462
 463
 464
 465   PostData(1, fOutputContainer);
 466
 467   printf("Init done \n") ;
 468
 469 }
 470
 471 //________________________________________________________________________
 472 void AliPHOSHijingEfficiency::UserExec(Option_t *)
 473 {
 474   // Main loop, called for each event
 475   // Analyze ESD/AOD
 476   fStack=0 ;
 477   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 478     if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
 479       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 480   }
 481
 482
 483   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 484   if (!event) {
 485     Printf("ERROR: Could not retrieve event");
 486     PostData(1, fOutputContainer);
 487     return;
 488   }
 489
 490   fRunNumber=ConvertRunNumber(event->GetRunNumber()) ;
 491   FillHistogram("hSelEvents",0.5,fRunNumber-0.5) ;
 492   FillHistogram("hTotSelEvents",0.5) ;
 493
 494   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 495   char skey[55] ;
 496   if(fEventCounter == 0) {
 497     // Initialize the PHOS geometry
 498     fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
 499
 500     //We have to apply re-calibration for pass1 LCH10h
 501     // Initialize decalibration factors in the form of the OCDB object
 502     const Float_t wDecalib = 0.055;
 503
 504     fPHOSCalibData = new AliPHOSCalibData();
 505     fPHOSCalibData->SetName("PhosCalibData");
 506
 507     for(Int_t module=1; module<=5; module++) {
 508       for(Int_t column=1; column<=56; column++) {
 509         for(Int_t row=1; row<=64; row++) {
 510           Float_t adcChannelEmc = gRandom->Gaus(1.,wDecalib);
 511           fPHOSCalibData->SetADCchannelEmc(module,column,row,adcChannelEmc);
 512         }
 513       }
 514     }
 515     for(Int_t mod=0; mod<5; mod++) {
 516       if(!event->GetPHOSMatrix(mod)) continue;
 517       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 518       Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
 519      }
 520
 521     fEventCounter++ ;
 522   }
 523
 524
 525   // Checks if we have a primary vertex
 526   // Get primary vertices form ESD
 527   const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
 528
 529   Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
 530   Double_t vtx5[3] ={0.,0.,0.};
 531
 532   vtx5[0] = esdVertex5->GetX();
 533   vtx5[1] = esdVertex5->GetY();
 534   vtx5[2] = esdVertex5->GetZ();
 535
 536
 537   FillHistogram("hZvertex",esdVertex5->GetZ(),fRunNumber-0.5);
 538   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 539     PostData(1, fOutputContainer);
 540     return;
 541   }
 542   FillHistogram("hSelEvents",1.5,fRunNumber-0.5) ;
 543   FillHistogram("hTotSelEvents",1.5) ;
 544
 545 /*
 546   if(event->IsPileupFromSPD()){
 547     PostData(1, fOutputContainer);
 548     return;
 549   }
 550 */
 551   FillHistogram("hSelEvents",2.5,fRunNumber-0.5) ;
 552   FillHistogram("hTotSelEvents",2.5) ;
 553
 554
 555   /*
 556   //Vtx class z-bin
 557   Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 558   if(zvtx<0)zvtx=0 ;
 559   if(zvtx>9)zvtx=9 ;
 560   */
 561   //No dependence on zVtx observed, save memory
 562   Int_t zvtx=0 ;
 563
 564   AliCentrality *centrality = event->GetCentrality();
 565   fCentrality=centrality->GetCentralityPercentile("V0M");
 566
 567   if( fCentrality <= 0. || fCentrality>80. ){
 568     PostData(1, fOutputContainer);
 569     return;
 570   }
 571
 572   FillHistogram("hSelEvents",3.5,fRunNumber-0.5) ;
 573   FillHistogram("hTotSelEvents",3.5) ;
 574
 575   if(fCentrality<5.)
 576     fCenBin=0 ;
 577   else if(fCentrality<10.)
 578     fCenBin=1 ;
 579   else if(fCentrality<20.)
 580     fCenBin=2 ;
 581   else if(fCentrality<40.)
 582     fCenBin=3 ;
 583   else if(fCentrality<60.)
 584     fCenBin=4 ;
 585   else if(fCentrality<80.)
 586     fCenBin=5 ;
 587   else{
 588     PostData(1, fOutputContainer);
 589     return;
 590   }
 591
 592   fRP=0. ;
 593
 594   FillHistogram("hSelEvents",4.5,fRunNumber-0.5) ;
 595   FillHistogram("hTotSelEvents",4.5) ;
 596   //All event selections done
 597   FillHistogram("hCentrality",fCentrality,fRunNumber-0.5) ;
 598   //Reaction plane is defined in the range (0;pi)
 599   //We have 10 bins
 600   Double_t averageRP = 0.;//fRPV0A+fRPV0C+fRP ;
 601   Int_t irp=Int_t(10.*(averageRP)/TMath::Pi());
 602   if(irp>9)irp=9 ;
 603
 604   if(!fPHOSEvents[zvtx][fCenBin][irp])
 605     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 606   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 607
 608   ProcessMC() ;
 609
 610   if(fPHOSEvent)
 611     fPHOSEvent->Clear() ;
 612   else
 613     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
 614
 615   //For re-calibration
 616   const Double_t logWeight=4.5 ;
 617
 618   TVector3 vertex(vtx5);
 619
 620   Int_t multClust = event->GetNumberOfCaloClusters();
 621   Int_t inPHOS=0,inMod1=0, inMod2=0, inMod3=0 ;
 622   Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy
 623
 624   AliESDCaloCells * cells = event->GetPHOSCells() ;
 625   FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
 626   FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
 627
 628
 629   TVector3 localPos ;
 630   for (Int_t i=0; i<multClust; i++) {
 631     AliESDCaloCluster *clu = event->GetCaloCluster(i);
 632     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 633
 634     Float_t  position[3];
 635     clu->GetPosition(position);
 636     TVector3 global(position) ;
 637     Int_t relId[4] ;
 638     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 639     Int_t mod  = relId[0] ;
 640     Int_t cellX = relId[2];
 641     Int_t cellZ = relId[3] ;
 642     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 643       continue ;
 644
 645
 646     FillHistogram("hCluEvsClu",clu->E(),clu->GetNCells()) ;
 647     FillHistogram("hCluEvsCluM",clu->E(),clu->GetM02()) ;
 648
 649     //Apply re-Calibreation
 650     AliPHOSEsdCluster cluPHOS1(*clu);
 651     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 652     Reclusterize(&cluPHOS1) ;
 653     cluPHOS1.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
 654     cluPHOS1.SetE(rnlin(cluPHOS1.E()));// Users's nonlinearity
 655     if(cluPHOS1.E()<0.3) continue;
 656     if(clu->GetNCells()<3) continue ;
 657     if(clu->GetM02()<0.2)   continue ;
 658
 659
 660     if(mod==1){
 661       avrgEm1+=cluPHOS1.E() ;
 662       inMod1++ ;
 663     }
 664     if(mod==2){
 665       avrgEm2+=cluPHOS1.E() ;
 666       inMod2++ ;
 667     }
 668     if(mod==3){
 669       avrgEm3+=cluPHOS1.E() ;
 670       inMod3++ ;
 671     }
 672
 673     TLorentzVector pv1 ;
 674     cluPHOS1.GetMomentum(pv1 ,vtx0);
 675
 676     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 677
 678     sprintf(skey,"hCluLowM%d",mod) ;
 679     FillHistogram(skey,cellX,cellZ,1.);
 680     if(cluPHOS1.E()>1.5){
 681       sprintf(skey,"hCluHighM%d",mod) ;
 682       FillHistogram(skey,cellX,cellZ,1.);
 683     }
 684
 685     if(inPHOS>=fPHOSEvent->GetSize()){
 686       fPHOSEvent->Expand(inPHOS+50) ;
 687     }
 688     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 689     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 690     ph->SetModule(mod) ;
 691     pv1*= ecore/pv1.E() ;
 692     ph->SetMomV2(&pv1) ;
 693     ph->SetNCells(clu->GetNCells());
 694     Double_t m02=0.,m20=0.;
 695     EvalLambdas(&cluPHOS1,0,m02, m20);
 696     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 697     EvalLambdas(&cluPHOS1,1,m02, m20);
 698     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 699
 700     Double_t distBC=clu->GetDistanceToBadChannel();
 701     if(distBC>2.)
 702       FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
 703       if(distBC>4.)
 704         FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
 705         if(distBC>6.)
 706           FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
 707
 708     if(ph->IsDispOK()){
 709       sprintf(skey,"hCluDispM%d",mod) ;
 710       FillHistogram(skey,cellX,cellZ,1.);
 711     }
 712     //Track matching
 713     Double_t dx=clu->GetTrackDx() ;
 714     Double_t dz=clu->GetTrackDz() ;
 715
 716     Bool_t cpvBit=kTRUE ; //No track matched by default
 717     Bool_t cpvBit2=kTRUE ; //More Strict criterion
 718     TArrayI * itracks = clu->GetTracksMatched() ;
 719     if(itracks->GetSize()>0){
 720       Int_t iTr = itracks->At(0);
 721       if(iTr>=0 && iTr<event->GetNumberOfTracks()){
 722         AliESDtrack *track = event->GetTrack(iTr) ;
 723         Double_t pt = track->Pt() ;
 724         Short_t charge = track->Charge() ;
 725         Double_t r=TestCPV(dx, dz, pt,charge) ;
 726         cpvBit=(r>2.) ;
 727         cpvBit2=(r>4.) ;
 728         FillHistogram("hCPVr",r,pv1.E(),fCentrality);
 729       }
 730     }
 731     ph->SetCPVBit(cpvBit) ;
 732     ph->SetCPV2Bit(cpvBit2) ;
 733     if(cpvBit){
 734       sprintf(skey,"hCluVetoM%d",mod) ;
 735       FillHistogram(skey,cellX,cellZ,1.);
 736     }
 737     ph->SetEMCx(float(cellX)) ;
 738     ph->SetEMCz(float(cellZ)) ;
 739     ph->SetLambdas(cluPHOS1.GetM20(),cluPHOS1.GetM02()) ;
 740     ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolded
 741     Bool_t sure=  kTRUE;
 742     Int_t primary=FindPrimary(clu,sure) ;  //номер праймари частицы в стеке
 743     ph->SetPrimary(primary) ;
 744     ph->SetWeight(PrimaryWeight(primary)) ;
 745     inPHOS++ ;
 746   }
 747
 748   FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
 749   if(inPHOS==0){
 750     PostData(1, fOutputContainer);
 751     fEventCounter++;
 752     return ;
 753   }
 754
 755
 756   for (Int_t i1=0; i1<inPHOS; i1++) {
 757     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 758     Double_t w=ph1->GetWeight() ;
 759     Double_t pt = ph1->Pt() ;
 760     Double_t ptv =ph1->GetMomV2()->Pt() ;
 761
 762     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
 763     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptv,w) ;
 764     if(ph1->IsntUnfolded()){
 765       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
 766     }
 767     if(ph1->IsCPVOK()){
 768       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
 769       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptv,w) ;
 770     }
 771     if(ph1->IsCPV2OK()){
 772       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
 773     }
 774     if(ph1->IsDispOK()){
 775       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
 776
 777       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptv,w) ;
 778       if(ph1->IsntUnfolded()){
 779         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
 780       }
 781       if(ph1->IsCPVOK()){
 782         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
 783         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptv,w) ;
 784       }
 785     }
 786     if(ph1->IsDisp2OK()){
 787       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
 788       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptv,w) ;
 789       if(ph1->IsCPVOK()){
 790         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
 791         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptv,w) ;
 792       }
 793     }
 794   }
 795
 796
 797   const Double_t alphaCut=0.7 ;
 798   //pi0
 799   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 800     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 801     Double_t w1=ph1->GetWeight() ;
 802     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 803       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 804       Double_t w2=ph2->GetWeight() ;
 805       Double_t w = TMath::Sqrt(w1*w2) ;
 806       TLorentzVector p12  = *ph1  + *ph2;
 807       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 808
 809       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 810       Double_t m = p12.M() ;
 811       Double_t pt = p12.Pt() ;
 812       Double_t mv = p12.M() ;
 813       Double_t ptv = p12.Pt() ;
 814
 815       FillHistogram(Form("hPi0All_cen%d",fCenBin),pt,m,w) ;
 816       FillHistogram(Form("hPi0Allcore_cen%d",fCenBin),ptv,mv,w) ;
 817       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 818         FillHistogram(Form("hPi0Allwou_cen%d",fCenBin),pt,m,w) ;
 819       }
 820
 821
 822       Double_t pt1 = ph1->Pt() ;
 823       Double_t ptv1 = ph1->GetMomV2()->Pt() ;
 824       Double_t pt2 = ph2->Pt() ;
 825       Double_t ptv2 = ph2->GetMomV2()->Pt() ;
 826       FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt1,w1) ;
 827       FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt2,w2) ;
 828       FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv1,w1) ;
 829       FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv2,w2) ;
 830       if(ph1->IsntUnfolded())
 831         FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt1,w1) ;
 832       if(ph2->IsntUnfolded())
 833         FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt2,w2) ;
 834       if(ph1->IsCPVOK()){
 835         FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt1,w1) ;
 836         FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv1,w1) ;
 837       }
 838       if(ph2->IsCPVOK()){
 839         FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt2,w2) ;
 840         FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv2,w2) ;
 841       }
 842       if(ph1->IsCPV2OK()){
 843         FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt1,w1) ;
 844       }
 845       if(ph2->IsCPV2OK()){
 846         FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt2,w2) ;
 847       }
 848       if(ph1->IsDispOK()){
 849         FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt1,w1) ;
 850         if(ph1->IsntUnfolded()){
 851           FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt1,w1) ;
 852         }
 853         FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv1,w1) ;
 854       }
 855       if(ph2->IsDispOK()){
 856         FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt2,w2) ;
 857         if(ph1->IsntUnfolded()){
 858           FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt2,w2) ;
 859         }
 860         FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv2,w2) ;
 861       }
 862       if(ph1->IsDisp2OK()){
 863         FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt1,w1) ;
 864         FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv1,w1) ;
 865       }
 866       if(ph2->IsDisp2OK()){
 867         FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt2,w2) ;
 868         FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv2,w2) ;
 869       }
 870       if(ph1->IsDispOK() && ph1->IsCPVOK()){
 871         FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt1,w1) ;
 872         FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv1,w1) ;
 873       }
 874       if(ph2->IsDispOK() && ph2->IsCPVOK()){
 875         FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt2,w2) ;
 876         FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv2,w2) ;
 877       }
 878       if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
 879         FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt1,w1) ;
 880         FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv1,w1) ;
 881       }
 882       if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
 883         FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt2,w2) ;
 884         FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv2,w2) ;
 885       }
 886
 887
 888       if(a<alphaCut){
 889         FillHistogram(Form("hPi0All_a07_cen%d",fCenBin),pt,m,w) ;
 890       }
 891
 892       if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 CPV
 893         FillHistogram(Form("hPi0CPV_cen%d",fCenBin),pt,m,w) ;
 894         FillHistogram(Form("hPi0CPVcore_cen%d",fCenBin),ptv,mv,w) ;
 895
 896         if(a<alphaCut){
 897           FillHistogram(Form("hPi0CPV_a07_cen%d",fCenBin),pt,m,w) ;
 898         }
 899       }
 900       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 901         FillHistogram(Form("hPi0CPV2_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
 902         if(a<alphaCut){
 903           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCenBin),pt,m,w) ;
 904         }
 905       }
 906       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 907         FillHistogram(Form("hPi0Disp2_cen%d",fCenBin),pt,m,w) ;
 908         FillHistogram(Form("hPi0Disp2core_cen%d",fCenBin),ptv,mv,w) ;
 909         if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
 910           FillHistogram(Form("hPi0Both2_cen%d",fCenBin),pt,m,w) ;
 911           FillHistogram(Form("hPi0Both2core_cen%d",fCenBin),ptv,mv,w) ;
 912           if(a<alphaCut){
 913             FillHistogram(Form("hPi0Both2_a07_cen%d",fCenBin),pt,m,w) ;
 914           }
 915         }
 916       }
 917
 918
 919       if(ph1->IsDispOK() && ph2->IsDispOK()){ //pi0 Disp
 920         FillHistogram(Form("hPi0Disp_cen%d",fCenBin),pt,m,w) ;
 921         FillHistogram(Form("hPi0Dispcore_cen%d",fCenBin),ptv,mv,w) ;
 922         if(a<alphaCut){
 923           FillHistogram(Form("hPi0Disp_a07_cen%d",fCenBin),pt,m,w) ;
 924         }
 925
 926         if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
 927           FillHistogram(Form("hPi0Both_cen%d",fCenBin),pt,m,w) ;
 928           FillHistogram(Form("hPi0Bothcore_cen%d",fCenBin),ptv,mv,w) ;
 929
 930           if(a<alphaCut){
 931             FillHistogram(Form("hPi0Both_a07_cen%d",fCenBin),pt,m,w) ;
 932           }
 933         }
 934       }
 935     } // end of loop i2
 936   } // end of loop i1
 937
 938   //now mixed
 939   for (Int_t i1=0; i1<inPHOS; i1++) {
 940     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 941     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 942       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 943       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 944         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 945         TLorentzVector p12  = *ph1  + *ph2;
 946         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 947
 948         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 949         Double_t m = p12.M() ;
 950         Double_t pt = p12.Pt() ;
 951         Double_t mv = p12.M() ;
 952         Double_t ptv = p12.Pt() ;
 953
 954         FillHistogram(Form("hMiPi0All_cen%d",fCenBin),pt,m) ;
 955         FillHistogram(Form("hMiPi0Allcore_cen%d",fCenBin),ptv,mv) ;
 956         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 957           FillHistogram(Form("hMiPi0Allwou_cen%d",fCenBin),pt,m) ;
 958         }
 959
 960         Double_t pt1 = ph1->Pt() ;
 961         Double_t ptv1 = ph1->GetMomV2()->Pt() ;
 962         Double_t pt2 = ph2->Pt() ;
 963         Double_t ptv2 = ph2->GetMomV2()->Pt() ;
 964         FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt1) ;
 965         FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt2) ;
 966         FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv1) ;
 967         FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv2) ;
 968         if(ph1->IsntUnfolded())
 969           FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt1) ;
 970         if(ph2->IsntUnfolded())
 971           FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt2) ;
 972         if(ph1->IsCPVOK()){
 973           FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt1) ;
 974           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv1) ;
 975         }
 976         if(ph2->IsCPVOK()){
 977           FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt2) ;
 978           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv2) ;
 979         }
 980         if(ph1->IsCPV2OK()){
 981           FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt1) ;
 982         }
 983         if(ph2->IsCPV2OK()){
 984           FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt2) ;
 985         }
 986         if(ph1->IsDispOK()){
 987           FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt1) ;
 988           if(ph1->IsntUnfolded()){
 989             FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt1) ;
 990           }
 991           FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv1) ;
 992         }
 993         if(ph2->IsDispOK()){
 994           FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt2) ;
 995           if(ph1->IsntUnfolded()){
 996             FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt2) ;
 997           }
 998           FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv2) ;
 999         }
1000         if(ph1->IsDisp2OK()){
1001           FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt1) ;
1002           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv1) ;
1003         }
1004         if(ph2->IsDisp2OK()){
1005           FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt2) ;
1006           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv2) ;
1007         }
1008         if(ph1->IsDispOK() && ph1->IsCPVOK()){
1009           FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt1) ;
1010           FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv1) ;
1011         }
1012         if(ph2->IsDispOK() && ph2->IsCPVOK()){
1013           FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt2) ;
1014           FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv2) ;
1015         }
1016         if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1017           FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt1) ;
1018           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv1) ;
1019         }
1020         if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1021           FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt2) ;
1022           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv2) ;
1023         }
1024
1025
1026
1027         if(a<alphaCut){
1028           FillHistogram(Form("hMiPi0All_a07_cen%d",fCenBin),pt,m) ;
1029         }
1030         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1031           FillHistogram(Form("hMiPi0CPV_cen%d",fCenBin),p12.M(), p12.Pt()) ;
1032           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCenBin),ptv,mv) ;
1033
1034           if(a<alphaCut){
1035             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCenBin),pt,m) ;
1036           }
1037         }
1038         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1039           FillHistogram(Form("hMiPi0CPV2_cen%d",fCenBin),pt,m) ;
1040
1041           if(a<alphaCut){
1042             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCenBin),pt,m) ;
1043           }
1044         }
1045         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1046           FillHistogram(Form("hMiPi0Disp2_cen%d",fCenBin),pt,m) ;
1047           FillHistogram(Form("hMiPi0Disp2core_cen%d",fCenBin),ptv,mv) ;
1048           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1049             FillHistogram(Form("hMiPi0Both2_cen%d",fCenBin),pt,m) ;
1050             FillHistogram(Form("hMiPi0Both2core_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1051           }
1052         }
1053         if(ph1->IsDispOK() && ph2->IsDispOK()){
1054
1055           FillHistogram(Form("hMiPi0Disp_cen%d",fCenBin),pt,m) ;
1056           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1057           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1058             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCenBin),p12.M(),p12.Pt()) ;
1059           }
1060
1061           if(a<alphaCut){
1062             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCenBin),pt,m) ;
1063           }
1064           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1065             FillHistogram(Form("hMiPi0Both_cen%d",fCenBin),pt,m) ;
1066             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCenBin),pv12.M(),pv12.Pt()) ;
1067
1068             if(a<alphaCut){
1069               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCenBin),pt,m) ;
1070             }
1071           }
1072         }
1073       } // end of loop i2
1074     }
1075   } // end of loop i1
1076
1077
1078   FillSecondaries() ;
1079
1080
1081   //Now we either add current events to stack or remove
1082   //If no photons in current event - no need to add it to mixed
1083   const Int_t kMixEvents[6]={5,5,5,10,10,30} ;
1084   if(fPHOSEvent->GetEntriesFast()>0){
1085     prevPHOS->AddFirst(fPHOSEvent) ;
1086     fPHOSEvent=0;
1087     if(prevPHOS->GetSize()>kMixEvents[fCenBin]){//Remove redundant events
1088       TClonesArray * tmp = static_cast<TClonesArray*>(prevPHOS->Last()) ;
1089       prevPHOS->RemoveLast() ;
1090       delete tmp ;
1091     }
1092   }
1093   // Post output data.
1094   PostData(1, fOutputContainer);
1095   fEventCounter++;
1096 }
1097 //________________________________________________________________________
1098 void AliPHOSHijingEfficiency::FillSecondaries(){
1099   //Sort secondaires
1100
1101   const Double_t rcut = 1. ;
1102   //Fill spectra of primary particles
1103   //with proper weight
1104   for(Int_t i=0; i<fStack->GetNtrack(); i++){
1105     TParticle * p = fStack->Particle(i) ;
1106     if(p->R()>rcut)
1107       continue ;
1108     if(TMath::Abs(p->Y())>0.5)
1109       continue ;
1110     Double_t w = PrimaryParticleWeight(p) ;
1111     Int_t primPdgCode=p->GetPdgCode() ;
1112       switch(primPdgCode){
1113         case  22: FillHistogram(Form("hPrimPhot_cen%d",fCenBin),p->Pt(),w);
1114                   break ;
1115         case  11:
1116         case -11:
1117                   FillHistogram(Form("hPrimEl_cen%d",fCenBin),p->Pt(),w);
1118                   break ;
1119         case  111:
1120                   FillHistogram(Form("hPrimPi0_cen%d",fCenBin),p->Pt(),w);
1121                   break ;
1122         case  221:
1123                   FillHistogram(Form("hPrimEta_cen%d",fCenBin),p->Pt(),w);
1124                   break ;
1125         case  211:
1126         case -211:
1127                   FillHistogram(Form("hPrimPipm_cen%d",fCenBin),p->Pt(),w);
1128                   break ;
1129         case  2212:  //p
1130                   FillHistogram(Form("hPrimP_cen%d",fCenBin),p->Pt(),w);
1131                   break ;
1132         case -2212:  //pbar
1133                   FillHistogram(Form("hPrimPbar_cen%d",fCenBin),p->Pt(),w);
1134                   break ;
1135         case  2112:  //n
1136                   FillHistogram(Form("hPrimN_cen%d",fCenBin),p->Pt(),w);
1137                   break ;
1138         case -2112:  //nbar
1139                   FillHistogram(Form("hPrimNbar_cen%d",fCenBin),p->Pt(),w);
1140                   break ;
1141         case  310:  //nbar
1142                   FillHistogram(Form("hPrimK0S_cen%d",fCenBin),p->Pt(),w);
1143                   break ;
1144         case  130:  //nbar
1145                   FillHistogram(Form("hPrimK0L_cen%d",fCenBin),p->Pt(),w);
1146                   break ;
1147         case  321:  //K+
1148         case -321:  //K-
1149                   FillHistogram(Form("hPrimKpm_cen%d",fCenBin),p->Pt(),w);
1150                   break ;
1151         default:           //other
1152                   FillHistogram(Form("hPrimOther_cen%d",fCenBin),p->Pt(),w);
1153       }
1154   }
1155
1156   //Origins of secondary pi0s
1157   for(Int_t i=0; i<fStack->GetNtrack(); i++){
1158     TParticle * p = fStack->Particle(i) ;
1159     if(p->GetPdgCode()!=111)
1160       continue ;
1161     FillHistogram("Vertex",p->Pt(),p->R());
1162     if(p->R()<rcut)
1163       continue ;
1164     Double_t phi=p->Phi() ;
1165     while(phi<0.)phi+=TMath::TwoPi() ;
1166     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1167     FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ;
1168     Double_t w = PrimaryParticleWeight(p) ;
1169     FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ;
1170   }
1171
1172   TLorentzVector p1;
1173
1174   Int_t inPHOS=fPHOSEvent->GetEntries() ;
1175   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1176     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
1177     Double_t w1=ph1->GetWeight() ;
1178     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1179       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
1180       TLorentzVector p12  = *ph1  + *ph2;
1181       Double_t w2=ph2->GetWeight() ;
1182       Double_t w = TMath::Sqrt(w1*w2) ;
1183       FillHistogram(Form("hParentAll_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1184       Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
1185       if(prim>-1){
1186         TParticle * particle = (TParticle *)fStack->Particle(prim);
1187         FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ;
1188
1189
1190         Int_t pdgCode=particle->GetPdgCode() ;
1191         if(pdgCode!=111){ //common parent - not pi0
1192           if(pdgCode==22)
1193             FillHistogram(Form("hParentGamma_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1194           else{
1195             if(pdgCode==11 || pdgCode==-11){
1196               FillHistogram(Form("hParentEl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1197             }
1198             else{
1199               if(InPi0mass(p12.M() ,p12.Pt())){
1200                 printf("Common parent: %d \n",pdgCode) ;
1201               }
1202               FillHistogram(Form("hParentOther_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1203             }
1204           }//Not photons
1205         }//Common parent not pi0
1206         else{ //common parent - pi0
1207           FillHistogram(Form("hParentPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1208           FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ;
1209           FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ;
1210           if(particle->R()<rcut && TMath::Abs(particle->Vz())<10.){
1211             FillHistogram(Form("hParentDirPi0_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1212             continue ;
1213           }
1214           //Common particle pi0, created off-vertex
1215           Int_t primPi0=particle->GetFirstMother();
1216           if(primPi0==-1){
1217             FillHistogram(Form("hParentPi0NoPrim_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
1218           }
1219           else{
1220             Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ;
1221             switch(primPdgCode){
1222             case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //eta
1223                       break ;
1224             case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //omega
1225                       break ;
1226             case  211:  //pi+-
1227             case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1228                       break ;
1229             case  321:  //K+-
1230             case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1231                       break ;
1232             case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0S
1233                       break ;
1234             case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // K0L
1235                       break ;
1236             case  2212:  //p
1237             case  2112:  //n
1238                       FillHistogram(Form("hParentPi0pn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1239                       break ;
1240             case -2212:  //pbar
1241             case -2112:  //nbar
1242                       FillHistogram(Form("hParentPi0antipn_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; // pn
1243                       break ;
1244             default:       //other
1245                       FillHistogram(Form("hParentPi0Other_cen%d",fCenBin),p12.M(),p12.Pt(),w) ; //
1246             }//switch
1247           }//pi0 with primary
1248         }//common parent - pi0
1249       }//there is common primary
1250     }//seond photon loop
1251   }//first photon loop
1252
1253
1254   //Now look at photon contaiminations
1255   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1256     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
1257     Int_t iprim = ph1->GetPrimary() ;
1258     if(iprim<0)
1259       FillAllHistograms(Form("hGammaNoPrim_cen%d",fCenBin),ph1) ; //
1260     else{
1261       //Find primary at vertex
1262       TParticle * primPHOS = fStack->Particle(iprim) ;
1263       Int_t iprimV=primPHOS->GetFirstMother();
1264       TParticle * primVtx = primPHOS ;
1265       while((iprimV>-1) && primVtx->R()>rcut){
1266         primVtx = fStack->Particle(iprimV) ;
1267         iprimV=primVtx->GetFirstMother();
1268       }
1269
1270       //photon
1271       Int_t primPdgCode=primVtx->GetPdgCode() ;
1272       switch(primPdgCode){
1273         case  22: FillAllHistograms("hGammaPhot",ph1);
1274                   break ;
1275         case  11:
1276         case -11:
1277                   FillAllHistograms("hGammaEl",ph1);
1278                   break ;
1279         case  111:
1280                   FillAllHistograms("hGammaPi0",ph1);
1281                   break ;
1282         case  221:
1283                   FillAllHistograms("hGammaEta",ph1);
1284                   break ;
1285         case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
1286                   break ;
1287         case  211:
1288         case -211:
1289                   FillAllHistograms("hGammaPipm",ph1);
1290                   //Find particle entered PHOS
1291                   if(primVtx == primPHOS)
1292                     FillAllHistograms("hGammaPipmDirect",ph1);
1293                   else{
1294                     Int_t primPdgPHOS=primPHOS->GetPdgCode() ;
1295                     if(primPdgPHOS==22){
1296                        FillAllHistograms("hGammaPipmGamma",ph1);
1297                        FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R());
1298                        FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R());
1299                        break ;
1300                     }
1301                     if(TMath::Abs(primPdgPHOS)==11){
1302                        FillAllHistograms("hGammaPipmEl",ph1);
1303                        FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
1304                        break ;
1305                     }
1306                     if(TMath::Abs(primPdgPHOS)==2212){
1307                        FillAllHistograms("hGammaPipmp",ph1);
1308                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1309                        break ;
1310                     }
1311                     if(TMath::Abs(primPdgPHOS)==2112){
1312                        FillAllHistograms("hGammaPipmn",ph1);
1313                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1314                        break ;
1315                     }
1316                     FillAllHistograms("hGammaPipmOther",ph1);
1317                     FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());
1318                   }
1319                   break ;
1320         case  2212:  //p
1321                   FillAllHistograms("hGammaP",ph1);
1322                   break ;
1323         case -2212:  //pbar
1324                   FillAllHistograms("hGammaPbar",ph1);
1325                   break ;
1326         case  2112:  //n
1327                   FillAllHistograms("hGammaN",ph1);
1328                   break ;
1329         case -2112:  //nbar
1330                   FillAllHistograms("hGammaNbar",ph1) ; // pn
1331                   break ;
1332         case  310:  //nbar
1333                   FillAllHistograms("hGammaK0S",ph1) ; // pn
1334                   break ;
1335         case  130:  //nbar
1336                   FillAllHistograms("hGammaK0L",ph1) ; // pn
1337                   break ;
1338         case  321:  //K+
1339         case -321:  //K-
1340                   FillAllHistograms("hGammaKpm",ph1) ; // pn
1341                   break ;
1342         case -323:
1343         case  323:
1344         case -313:
1345         case  313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
1346                   break ;
1347
1348         case -2224 : //Deltas
1349         case  2224 : //Deltas
1350         case -2214 : //Deltas
1351         case  2214 : //Deltas
1352         case -2114 : //Deltas
1353         case  2114 : //Deltas
1354         case -1114 : //Deltas
1355         case  1114 : //Deltas
1356                   FillAllHistograms("hGammaDelta",ph1) ; // pn
1357                   break ;
1358         default:           //other
1359             if(primVtx->GetPDG()->Charge())
1360               FillAllHistograms("hGammaOtherCharged",ph1) ; //
1361             else
1362               FillAllHistograms("hGammaOtherNeutral",ph1) ; //
1363       }
1364     }
1365
1366   }//single photons
1367
1368 }
1369 //________________________________________________________________________
1370 void AliPHOSHijingEfficiency::Terminate(Option_t *)
1371 {
1372   // Draw result to the screen
1373   // Called once at the end of the query
1374
1375 }
1376
1377 //________________________________________________________________________
1378 Bool_t AliPHOSHijingEfficiency::IsGoodChannel(const char * det, Int_t mod, Int_t ix, Int_t iz)
1379 {
1380   //Check if this channel belogs to the good ones
1381
1382   if(strcmp(det,"PHOS")==0){
1383     if(mod>5 || mod<1){
1384       AliError(Form("No bad map for PHOS module %d ",mod)) ;
1385       return kTRUE ;
1386     }
1387     if(!fPHOSBadMap[mod]){
1388       AliError(Form("No Bad map for PHOS module %d",mod)) ;
1389       return kTRUE ;
1390     }
1391     if(fPHOSBadMap[mod]->GetBinContent(ix,iz)>0)
1392       return kFALSE ;
1393     else
1394       return kTRUE ;
1395   }
1396   else{
1397     AliError(Form("Can not find bad channels for detector %s ",det)) ;
1398   }
1399   return kTRUE ;
1400 }
1401 //_____________________________________________________________________________
1402 void AliPHOSHijingEfficiency::FillAllHistograms(const char * key,AliCaloPhoton * ph)const{
1403   //Fill All PID histograms
1404
1405   Double_t w=ph->GetWeight() ;
1406   Double_t pt = ph->Pt() ;
1407   Double_t ptC=ph->GetMomV2()->Pt() ;
1408   FillHistogram(Form("%s_All_cen%d",key,fCenBin),pt,w) ;
1409   FillHistogram(Form("%s_Allcore_cen%d",key,fCenBin),ptC,w) ;
1410   if(ph->IsCPVOK()){
1411     FillHistogram(Form("%s_CPV_cen%d",key,fCenBin),pt,w) ;
1412     FillHistogram(Form("%s_CPVcore_cen%d",key,fCenBin),ptC,w) ;
1413   }
1414   if(ph->IsCPV2OK()){
1415     FillHistogram(Form("%s_CPV2_cen%d",key,fCenBin),pt,w) ;
1416     FillHistogram(Form("%s_CPV2core_cen%d",key,fCenBin),ptC,w) ;
1417   }
1418   if(ph->IsDispOK()){
1419     FillHistogram(Form("%s_Disp_cen%d",key,fCenBin),pt,w) ;
1420     FillHistogram(Form("%s_Dispcore_cen%d",key,fCenBin),ptC,w) ;
1421     if(ph->IsDisp2OK()){
1422       FillHistogram(Form("%s_Disp2_cen%d",key,fCenBin),pt,w) ;
1423       FillHistogram(Form("%s_Disp2core_cen%d",key,fCenBin),ptC,w) ;
1424     }
1425     if(ph->IsCPVOK()){
1426       FillHistogram(Form("%s_Both_cen%d",key,fCenBin),pt,w) ;
1427       FillHistogram(Form("%s_Bothcore_cen%d",key,fCenBin),ptC,w) ;
1428     }
1429   }
1430 }
1431 //_____________________________________________________________________________
1432 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x)const{
1433   //FillHistogram
1434   TObject * tmp = fOutputContainer->FindObject(key) ;
1435   if(!tmp){
1436     AliInfo(Form("can not find histogram <%s> ",key)) ;
1437     return ;
1438   }
1439   if(tmp->IsA() == TClass::GetClass("TH1I")){
1440     ((TH1I*)tmp)->Fill(x) ;
1441     return ;
1442   }
1443   if(tmp->IsA() == TClass::GetClass("TH1F")){
1444     ((TH1F*)tmp)->Fill(x) ;
1445     return ;
1446   }
1447   if(tmp->IsA() == TClass::GetClass("TH1D")){
1448     ((TH1D*)tmp)->Fill(x) ;
1449     return ;
1450   }
1451   AliInfo(Form("can not find 1D histogram <%s> ",key)) ;
1452 }
1453 //_____________________________________________________________________________
1454 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y)const{
1455   //FillHistogram
1456   TObject * tmp = fOutputContainer->FindObject(key) ;
1457   if(!tmp){
1458     AliInfo(Form("can not find histogram <%s> ",key)) ;
1459     return ;
1460   }
1461   if(tmp->IsA() == TClass::GetClass("TH1F")){
1462     ((TH1F*)tmp)->Fill(x,y) ;
1463     return ;
1464   }
1465   if(tmp->IsA() == TClass::GetClass("TH2F")){
1466     ((TH2F*)tmp)->Fill(x,y) ;
1467     return ;
1468   }
1469   AliError(Form("Calling FillHistogram with 2 parameters for histo <%s> of type %s",key,tmp->IsA()->GetName())) ;
1470 }
1471
1472 //_____________________________________________________________________________
1473 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z) const{
1474   //Fills 1D histograms with key
1475   TObject * tmp = fOutputContainer->FindObject(key) ;
1476   if(!tmp){
1477     AliInfo(Form("can not find histogram <%s> ",key)) ;
1478     return ;
1479   }
1480   if(tmp->IsA() == TClass::GetClass("TH2F")){
1481     ((TH2F*)tmp)->Fill(x,y,z) ;
1482     return ;
1483   }
1484   if(tmp->IsA() == TClass::GetClass("TH3F")){
1485     ((TH3F*)tmp)->Fill(x,y,z) ;
1486     return ;
1487   }
1488 }
1489 //_____________________________________________________________________________
1490 void AliPHOSHijingEfficiency::FillHistogram(const char * key,Double_t x,Double_t y, Double_t z, Double_t w) const{
1491   //Fills 1D histograms with key
1492   TObject * tmp = fOutputContainer->FindObject(key) ;
1493   if(!tmp){
1494     AliInfo(Form("can not find histogram <%s> ",key)) ;
1495     return ;
1496   }
1497   if(tmp->IsA() == TClass::GetClass("TH3F")){
1498     ((TH3F*)tmp)->Fill(x,y,z,w) ;
1499     return ;
1500   }
1501 }
1502
1503 //___________________________________________________________________________
1504 Int_t AliPHOSHijingEfficiency::ConvertRunNumber(Int_t run){
1505
1506   switch(run){
1507   case  139517 : return 137;
1508   case  139514 : return 136;
1509   case  139513 : return 135;
1510   case  139511 : return 134;
1511   case  139510 : return 133;
1512   case  139507 : return 132;
1513   case  139505 : return 131;
1514   case  139504 : return 130;
1515   case  139503 : return 129;
1516   case  139470 : return 128;
1517   case  139467 : return 127;
1518   case  139466 : return 126;
1519   case  139465 : return 125;
1520   case  139440 : return 124;
1521   case  139439 : return 123;
1522   case  139438 : return 122;
1523   case  139437 : return 121;
1524   case  139360 : return 120;
1525   case  139329 : return 119;
1526   case  139328 : return 118;
1527   case  139314 : return 117;
1528   case  139311 : return 116;
1529   case  139310 : return 115;
1530   case  139309 : return 114;
1531   case  139308 : return 113;
1532   case  139173 : return 112;
1533   case  139172 : return 111;
1534   case  139110 : return 110;
1535   case  139107 : return 109;
1536   case  139105 : return 108;
1537   case  139104 : return 107;
1538   case  139042 : return 106;
1539   case  139038 : return 105;
1540   case  139037 : return 104;
1541   case  139036 : return 103;
1542   case  139029 : return 102;
1543   case  139028 : return 101;
1544   case  138983 : return 100;
1545   case  138982 : return 99;
1546   case  138980 : return 98;
1547   case  138979 : return 97;
1548   case  138978 : return 96;
1549   case  138977 : return 95;
1550   case  138976 : return 94;
1551   case  138973 : return 93;
1552   case  138972 : return 92;
1553   case  138965 : return 91;
1554   case  138924 : return 90;
1555   case  138872 : return 89;
1556   case  138871 : return 88;
1557   case  138870 : return 87;
1558   case  138837 : return 86;
1559   case  138830 : return 85;
1560   case  138828 : return 84;
1561   case  138826 : return 83;
1562   case  138796 : return 82;
1563   case  138795 : return 81;
1564   case  138742 : return 80;
1565   case  138732 : return 79;
1566   case  138730 : return 78;
1567   case  138666 : return 77;
1568   case  138662 : return 76;
1569   case  138653 : return 75;
1570   case  138652 : return 74;
1571   case  138638 : return 73;
1572   case  138624 : return 72;
1573   case  138621 : return 71;
1574   case  138583 : return 70;
1575   case  138582 : return 69;
1576   case  138579 : return 68;
1577   case  138578 : return 67;
1578   case  138534 : return 66;
1579   case  138469 : return 65;
1580   case  138442 : return 64;
1581   case  138439 : return 63;
1582   case  138438 : return 62;
1583   case  138396 : return 61;
1584   case  138364 : return 60;
1585   case  138359 : return 59;
1586   case  138275 : return 58;
1587   case  138225 : return 57;
1588   case  138201 : return 56;
1589   case  138200 : return 55;
1590   case  138197 : return 54;
1591   case  138192 : return 53;
1592   case  138190 : return 52;
1593   case  138154 : return 51;
1594   case  138153 : return 50;
1595   case  138151 : return 49;
1596   case  138150 : return 48;
1597   case  138126 : return 47;
1598   case  138125 : return 46;
1599   case  137848 : return 45;
1600   case  137847 : return 44;
1601   case  137844 : return 43;
1602   case  137843 : return 42;
1603   case  137752 : return 41;
1604   case  137751 : return 40;
1605   case  137748 : return 39;
1606   case  137724 : return 38;
1607   case  137722 : return 37;
1608   case  137718 : return 36;
1609   case  137704 : return 35;
1610   case  137693 : return 34;
1611   case  137692 : return 33;
1612   case  137691 : return 32;
1613   case  137689 : return 31;
1614   case  137686 : return 30;
1615   case  137685 : return 29;
1616   case  137639 : return 28;
1617   case  137638 : return 27;
1618   case  137608 : return 26;
1619   case  137595 : return 25;
1620   case  137549 : return 24;
1621   case  137546 : return 23;
1622   case  137544 : return 22;
1623   case  137541 : return 21;
1624   case  137539 : return 20;
1625   case  137531 : return 19;
1626   case  137530 : return 18;
1627   case  137443 : return 17;
1628   case  137441 : return 16;
1629   case  137440 : return 15;
1630   case  137439 : return 14;
1631   case  137434 : return 13;
1632   case  137432 : return 12;
1633   case  137431 : return 11;
1634   case  137430 : return 10;
1635   case  137366 : return 9;
1636   case  137243 : return 8;
1637   case  137236 : return 7;
1638   case  137235 : return 6;
1639   case  137232 : return 5;
1640   case  137231 : return 4;
1641   case  137165 : return 3;
1642   case  137162 : return 2;
1643   case  137161 : return 1;
1644   default : return 199;
1645   }
1646
1647 }
1648 //_____________________________________________________________________________
1649 Bool_t AliPHOSHijingEfficiency::TestLambda(Double_t pt,Double_t l1,Double_t l2){
1650
1651   //For R=3.5
1652    Double_t  l1Mean  = 1.170014 -0.059465/(1.+0.019343*pt+0.147544*pt*pt) ;
1653    Double_t  l2Mean = 1.626270 + 0.761554*exp(-1.213839*pt)-0.020027*pt ;
1654    Double_t  l1Sigma = 0.133409 + 0.261307*exp(-0.636874*pt)-0.002849*pt ;
1655    Double_t  l2Sigma = 0.289698 + 0.459400*exp(-1.214242*pt)-0.012578*pt ;
1656    Double_t  c=-0.124103 ;
1657 /*
1658   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1659   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1660   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1661   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1662   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1663 */
1664   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1665               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1666               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1667   return (R2<2.5*2.5) ;
1668
1669 }
1670 //_____________________________________________________________________________
1671 Bool_t AliPHOSHijingEfficiency::TestLambda2(Double_t pt,Double_t l1,Double_t l2){
1672
1673 //For R=4.5
1674   Double_t   l1Mean  = 1.150200 + 0.097886/(1.+1.486645*pt+0.000038*pt*pt) ;
1675   Double_t   l2Mean = 1.574706 + 0.997966*exp(-0.895075*pt)-0.010666*pt ;
1676   Double_t   l1Sigma = 0.100255 + 0.337177*exp(-0.517684*pt)+0.001170*pt ;
1677   Double_t   l2Sigma = 0.232580 + 0.573401*exp(-0.735903*pt)-0.002325*pt ;
1678   Double_t   c = -0.110983 -0.017353/(1.-1.836995*pt+0.934517*pt*pt) ;
1679
1680 /*
1681   Double_t l2Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt) ;
1682   Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt ;
1683   Double_t l2Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
1684   Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
1685   Double_t c=-0.35-0.550*TMath::Exp(-0.390730*pt) ;
1686 */
1687   Double_t R2=0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
1688               0.5*(l2-l2Mean)*(l2-l2Mean)/l2Sigma/l2Sigma +
1689               0.5*c*(l1-l1Mean)*(l2-l2Mean)/l1Sigma/l2Sigma ;
1690   return (R2<2.5*2.5) ;
1691
1692 }
1693 //____________________________________________________________________________
1694 Double_t AliPHOSHijingEfficiency::TestCPV(Double_t dx, Double_t dz, Double_t pt, Int_t charge){
1695   //Parameterization of LHC10h period
1696   //_true if neutral_
1697
1698   Double_t meanX=0;
1699   Double_t meanZ=0.;
1700   Double_t sx=TMath::Min(5.4,2.59719e+02*TMath::Exp(-pt/1.02053e-01)+
1701               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);
1702   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) ;
1703   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
1704   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1705
1706   if(mf<0.){ //field --
1707     meanZ = -0.468318 ;
1708     if(charge>0)
1709       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)) ;
1710     else
1711       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)) ;
1712   }
1713   else{ //Field ++
1714     meanZ= -0.468318;
1715     if(charge>0)
1716       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)) ;
1717     else
1718       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)) ;
1719   }
1720
1721   Double_t rz=(dz-meanZ)/sz ;
1722   Double_t rx=(dx-meanX)/sx ;
1723   return TMath::Sqrt(rx*rx+rz*rz) ;
1724 }
1725 //____________________________________________________________________________
1726 Bool_t AliPHOSHijingEfficiency::TestTOF(Double_t t, Double_t e){
1727
1728   Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1729                              +2.24611e-09*2.24611e-09/e
1730                              +5.65054e-09*5.65054e-09/e/e) ;
1731   sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1732   Double_t mean=1.1e-9 ;
1733   if(TMath::Abs(t-mean)<2.*sigma)
1734     return kTRUE ;
1735   else
1736     if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1737       return kTRUE ;
1738
1739   return kFALSE ;
1740 }
1741 //____________________________________________________________________________
1742 Double_t  AliPHOSHijingEfficiency::CoreEnergy(AliESDCaloCluster * clu){
1743   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1744
1745   //Can not use already calculated coordinates?
1746   //They have incidence correction...
1747   const Double_t distanceCut =3.5 ;
1748   const Double_t logWeight=4.5 ;
1749
1750   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1751 // Calculates the center of gravity in the local PHOS-module coordinates
1752   Float_t wtot = 0;
1753   Double_t xc[100]={0} ;
1754   Double_t zc[100]={0} ;
1755   Double_t x = 0 ;
1756   Double_t z = 0 ;
1757   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1758   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1759     Int_t relid[4] ;
1760     Float_t xi ;
1761     Float_t zi ;
1762     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1763     fPHOSGeo->RelPosInModule(relid, xi, zi);
1764     xc[iDigit]=xi ;
1765     zc[iDigit]=zi ;
1766     if (clu->E()>0 && elist[iDigit]>0) {
1767       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1768       x    += xc[iDigit] * w ;
1769       z    += zc[iDigit] * w ;
1770       wtot += w ;
1771     }
1772   }
1773   if (wtot>0) {
1774     x /= wtot ;
1775     z /= wtot ;
1776   }
1777   Double_t coreE=0. ;
1778   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1779     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1780     if(distance < distanceCut)
1781       coreE += elist[iDigit] ;
1782   }
1783   //Apply non-linearity correction
1784   return rnlin(coreE) ;
1785 }
1786 //____________________________________________________________________________
1787 Bool_t  AliPHOSHijingEfficiency::AreNeibors(Int_t id1,Int_t id2){
1788
1789   Int_t relid1[4] ;
1790   fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
1791
1792   Int_t relid2[4] ;
1793   fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
1794
1795   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
1796     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
1797     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
1798
1799     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
1800       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
1801       return 1 ;
1802     }
1803     else {
1804       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
1805         return 0; //  Difference in row numbers is too large to look further
1806     }
1807     return 0 ;
1808
1809   }
1810
1811   return 0 ;
1812 }
1813 //____________________________________________________________________________
1814 void  AliPHOSHijingEfficiency::Reclusterize(AliESDCaloCluster * clu){
1815   //Re-clusterize to make continues cluster
1816
1817   const Int_t oldMulDigit=clu->GetNCells() ;
1818   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1819   UShort_t * dlist = clu->GetCellsAbsId();
1820
1821   Int_t index[oldMulDigit] ;
1822   Bool_t used[oldMulDigit] ;
1823   for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
1824   Int_t inClu=0 ;
1825   Double_t eMax=0. ;
1826   //find maximum
1827   for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
1828     if(eMax<elist[iDigit]){
1829       eMax=elist[iDigit];
1830       index[0]=iDigit ;
1831       inClu=1 ;
1832     }
1833   }
1834   if(inClu==0){ //empty cluster
1835     return ;
1836   }
1837   used[index[0]]=kTRUE ; //mark as used
1838   for(Int_t i=0; i<inClu; i++){
1839     for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
1840        if(used[iDigit]) //already used
1841          continue ;
1842        if(AreNeibors(dlist[index[i]],dlist[iDigit])){
1843          index[inClu]= iDigit ;
1844          inClu++ ;
1845          used[iDigit]=kTRUE ;
1846        }
1847     }
1848   }
1849
1850   if(inClu==oldMulDigit) //no need to modify
1851     return ;
1852
1853   clu->SetNCells(inClu);
1854   //copy
1855   UShort_t tmpD[oldMulDigit] ;
1856   Double_t tmpE[oldMulDigit] ;
1857   for(Int_t i=0; i<oldMulDigit; i++){
1858     tmpD[i]=dlist[i] ;
1859     tmpE[i]=elist[i] ;
1860   }
1861   //change order of digits in list so that
1862   //first inClu cells were true ones
1863   for(Int_t i=0; i<inClu; i++){
1864     dlist[i]=tmpD[index[i]] ;
1865     elist[i]=tmpE[index[i]] ;
1866   }
1867
1868
1869 }
1870
1871 //____________________________________________________________________________
1872 void  AliPHOSHijingEfficiency::EvalLambdas(AliESDCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1873   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1874
1875   Double_t rCut=0. ;
1876   if(iR==0)
1877     rCut=3.5 ;
1878   else
1879     rCut=4.5 ;
1880
1881
1882   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1883 // Calculates the center of gravity in the local PHOS-module coordinates
1884   Float_t wtot = 0;
1885   Double_t xc[100]={0} ;
1886   Double_t zc[100]={0} ;
1887   Double_t x = 0 ;
1888   Double_t z = 0 ;
1889   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1890   const Double_t logWeight=4.5 ;
1891   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1892     Int_t relid[4] ;
1893     Float_t xi ;
1894     Float_t zi ;
1895     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1896     fPHOSGeo->RelPosInModule(relid, xi, zi);
1897     xc[iDigit]=xi ;
1898     zc[iDigit]=zi ;
1899     if (clu->E()>0 && elist[iDigit]>0) {
1900       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1901       x    += xc[iDigit] * w ;
1902       z    += zc[iDigit] * w ;
1903       wtot += w ;
1904     }
1905   }
1906   if (wtot>0) {
1907     x /= wtot ;
1908     z /= wtot ;
1909   }
1910
1911   wtot = 0. ;
1912   Double_t dxx  = 0.;
1913   Double_t dzz  = 0.;
1914   Double_t dxz  = 0.;
1915   Double_t xCut = 0. ;
1916   Double_t zCut = 0. ;
1917   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1918     if (clu->E()>0 && elist[iDigit]>0.) {
1919         Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1920         Double_t xi= xc[iDigit] ;
1921         Double_t zi= zc[iDigit] ;
1922         if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1923           xCut += w * xi ;
1924           zCut += w * zi ;
1925           dxx  += w * xi * xi ;
1926           dzz  += w * zi * zi ;
1927           dxz  += w * xi * zi ;
1928           wtot += w ;
1929         }
1930     }
1931
1932   }
1933   if (wtot>0) {
1934     xCut/= wtot ;
1935     zCut/= wtot ;
1936     dxx /= wtot ;
1937     dzz /= wtot ;
1938     dxz /= wtot ;
1939     dxx -= xCut * xCut ;
1940     dzz -= zCut * zCut ;
1941     dxz -= xCut * zCut ;
1942
1943     m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1944     m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1945   }
1946   else {
1947     m20=m02=0.;
1948   }
1949
1950 }
1951 //___________________________________________________________________________
1952 void AliPHOSHijingEfficiency::ProcessMC(){
1953   //fill histograms for efficiensy etc. calculation
1954   const Double_t rcut = 1. ; //cut for primary particles
1955   //---------First pi0/eta-----------------------------
1956   char partName[10] ;
1957   char hkey[55] ;
1958
1959   if(!fStack) return ;
1960   for(Int_t i=0;i<fStack->GetNtrack();i++){
1961      TParticle* particle =  fStack->Particle(i);
1962     if(particle->GetPdgCode() == 111)
1963       snprintf(partName,10,"pi0") ;
1964     else
1965       if(particle->GetPdgCode() == 221)
1966         snprintf(partName,10,"eta") ;
1967       else
1968         if(particle->GetPdgCode() == 22)
1969            snprintf(partName,10,"gamma") ;
1970         else
1971            continue ;
1972
1973     //Primary particle
1974     Double_t r=particle->R() ;
1975     Double_t pt = particle->Pt() ;
1976     //Distribution over vertex
1977     FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
1978
1979     if(r >rcut)
1980       continue ;
1981
1982     //Total number of pi0 with creation radius <1 cm
1983     Double_t w = PrimaryParticleWeight(particle) ;
1984     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
1985     FillHistogram(hkey,pt,w) ;
1986     if(TMath::Abs(particle->Y())<0.12){
1987       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
1988       FillHistogram(hkey,pt,w) ;
1989     }
1990
1991     snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
1992     FillHistogram(hkey,particle->Y(),w) ;
1993
1994     Double_t phi=particle->Phi() ;
1995     while(phi<0.)phi+=TMath::TwoPi() ;
1996     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1997     snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
1998     FillHistogram(hkey,phi,w) ;
1999   }
2000 }
2001 //___________________________________________________________________________
2002 Int_t AliPHOSHijingEfficiency::FindPrimary(AliESDCaloCluster*clu,  Bool_t&sure){
2003   //Finds primary and estimates if it unique one?
2004   //First check can it be photon/electron
2005   const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
2006   Int_t n=clu->GetNLabels() ;
2007   for(Int_t i=0;  i<n;  i++){
2008     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
2009     Int_t pdg = p->GetPdgCode() ;
2010     if(pdg==22  ||  pdg==11 || pdg == -11){
2011       if(p->Energy()>emFraction*clu->E()){
2012         sure=kTRUE ;
2013         return clu->GetLabelAt(i);
2014       }
2015     }
2016   }
2017
2018   Double_t*  Ekin=  new  Double_t[n] ;
2019   for(Int_t i=0;  i<n;  i++){
2020     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
2021     Ekin[i]=p->P() ;  // estimate of kinetic energy
2022     if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
2023       Ekin[i]+=1.8  ;  //due to annihilation
2024     }
2025   }
2026   Int_t iMax=0;
2027   Double_t eMax=0.,eSubMax=0. ;
2028   for(Int_t i=0;  i<n;  i++){
2029     if(Ekin[i]>eMax){
2030       eSubMax=eMax;
2031       eMax=Ekin[i];
2032       iMax=i;
2033     }
2034   }
2035   if(eSubMax>0.8*eMax)//not obvious primary
2036     sure=kFALSE;
2037   else
2038     sure=kTRUE;
2039   delete  Ekin;
2040   return  clu->GetLabelAt(iMax) ;
2041 }
2042 //___________________________________________________________________________
2043 Double_t AliPHOSHijingEfficiency::PrimaryWeight(Int_t primary){
2044    //Check who is the primary and introduce weight to correct primary spectrum
2045
2046   if(primary<0 || primary>=fStack->GetNtrack())
2047     return 1 ;
2048   //trace primaries up to IP
2049   TParticle* particle =  fStack->Particle(primary);
2050   Double_t r=particle->R() ;
2051   Int_t mother = particle->GetFirstMother() ;
2052   while(mother>-1){
2053     if(r<1. && particle->GetPdgCode()==111)
2054       break ;
2055     particle =  fStack->Particle(mother);
2056     mother = particle->GetFirstMother() ;
2057     r=particle->R() ;
2058   }
2059
2060   return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
2061 }
2062 //________________________________________________________________________
2063 Double_t AliPHOSHijingEfficiency::PrimaryParticleWeight(TParticle * particle){
2064
2065   Int_t pdg = particle->GetPdgCode() ;
2066   Int_t type=0 ;
2067   if(pdg == 111 || TMath::Abs(pdg)==211){
2068     type =1 ;
2069   }
2070   else{
2071     if(TMath::Abs(pdg)<1000){ //Kaon-like
2072       type =2 ;
2073     }
2074     else
2075       type = 3;  //baryons
2076   }
2077
2078   Double_t x = particle->Pt() ;
2079   if(type==1){
2080    if(fCenBin==0) //0-5
2081      return (1.662990+1.140890*x-0.192088*x*x)/(1.-0.806630*x+0.304771*x*x)+0.141690*x ;
2082    if(fCenBin==1) //5-10
2083      return (1.474351+0.791492*x-0.066369*x*x)/(1.-0.839338*x+0.317312*x*x)+0.093289*x ;
2084    if(fCenBin==2) //10-20
2085      return (1.174728+0.959681*x-0.137695*x*x)/(1.-0.788873*x+0.299538*x*x)+0.128759*x ;
2086    if(fCenBin==3) //20-40
2087      return (0.927335+0.475349*x+0.004364*x*x)/(1.-0.817966*x+0.309787*x*x)+0.086899*x ;
2088    if(fCenBin==4) //40-60
2089      return (0.676878+0.190680*x+0.077031*x*x)/(1.-0.790623*x+0.305183*x*x)+0.064510*x ;
2090    if(fCenBin==5) //60-80
2091      return (0.684726-0.606262*x+0.409963*x*x)/(1.-1.080061*x+0.456933*x*x)+0.005151*x ;
2092   }
2093   if(type==2){
2094    if(fCenBin==0) //0-5
2095      return (-0.417131+2.253936*x-0.337731*x*x)/(1.-0.909892*x+0.316820*x*x)+0.157312*x ;
2096    if(fCenBin==1) //5-10
2097      return (-0.352275+1.844466*x-0.248598*x*x)/(1.-0.897048*x+0.316462*x*x)+0.132461*x ;
2098    if(fCenBin==2) //10-20
2099      return (-0.475481+1.975108*x-0.336013*x*x)/(1.-0.801028*x+0.276705*x*x)+0.188164*x ;
2100    if(fCenBin==3) //20-40
2101      return (-0.198954+1.068789*x-0.103540*x*x)/(1.-0.848354*x+0.299209*x*x)+0.112939*x ;
2102    if(fCenBin==4) //40-60
2103      return (-0.111052+0.664041*x-0.019717*x*x)/(1.-0.804916*x+0.300779*x*x)+0.095784*x ;
2104    if(fCenBin==5) //0-5
2105      return (0.202788-0.439832*x+0.564585*x*x)/(1.-1.254029*x+0.679444*x*x)+0.016235*x ;
2106   }
2107   if(type==3){
2108    if(fCenBin==0) //0-5
2109      return (-1.312732+2.743568*x-0.375775*x*x)/(1.-0.717533*x+0.164694*x*x)+0.164445*x ;
2110    if(fCenBin==1) //5-10
2111      return (-1.229425+2.585889*x-0.330164*x*x)/(1.-0.715892*x+0.167386*x*x)+0.133085*x ;
2112    if(fCenBin==2) //10-20
2113      return (-1.135677+2.397489*x-0.320355*x*x)/(1.-0.709312*x+0.164350*x*x)+0.146095*x ;
2114    if(fCenBin==3) //20-40
2115      return (-0.889993+1.928263*x-0.220785*x*x)/(1.-0.715991*x+0.174729*x*x)+0.095098*x ;
2116    if(fCenBin==4) //40-60
2117      return (-0.539237+1.329118*x-0.115439*x*x)/(1.-0.722906*x+0.186832*x*x)+0.059267*x ;
2118    if(fCenBin==5) //60-80
2119      return (-0.518126+1.327628*x-0.130881*x*x)/(1.-0.665649*x+0.184300*x*x)+0.081701*x ;
2120   }
2121   return 1. ;
2122 }
2123 //________________________________________________________________________
2124 Int_t AliPHOSHijingEfficiency::FindCommonParent(Int_t iPart, Int_t jPart){
2125   //check if there is a common parent for particles i and j
2126   // -1: no common parent or wrong iPart/jPart
2127
2128   if(iPart==-1 || iPart>=fStack->GetNtrack() ||
2129      jPart==-1 || jPart>=fStack->GetNtrack()) return -1;
2130
2131   Int_t iprim1=iPart;
2132   while(iprim1>-1){
2133      Int_t iprim2=jPart;
2134      while(iprim2>-1){
2135        if(iprim1==iprim2)
2136          return iprim1 ;
2137        iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother();
2138      }
2139      iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother();
2140   }
2141   return -1;
2142 }
2143 //________________________________________________________________________
2144 Bool_t AliPHOSHijingEfficiency::HaveParent(Int_t iPart, Int_t pdgParent){
2145   //check if there is a common parent for particles i and j
2146   // -1: no common parent or wrong iPart/jPart
2147
2148   if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1;
2149
2150   Int_t iprim1=iPart;
2151   while(iprim1>-1){
2152     TParticle * tmp = fStack->Particle(iprim1) ;
2153     if(tmp->GetPdgCode()==pdgParent)
2154       return kTRUE ;
2155     iprim1=tmp->GetFirstMother();
2156   }
2157   return kFALSE;
2158 }
2159 //________________________________________________________________________
2160 Bool_t AliPHOSHijingEfficiency::InPi0mass(Double_t m, Double_t /*pt*/){
2161
2162  return TMath::Abs(m-0.135)<0.007*2.5 ;
2163 }