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   for(Int_t i=0; i<6; i++){
  85     fPHOSBadMap[i]=new TH2I(Form("PHOS_BadMap_mod%d",i),"Bad Modules map",64,0.,64.,56,0.,56.) ;
  86   }
  87
  88
  89 }
  90
  91 //________________________________________________________________________
  92 void AliPHOSHijingEfficiency::UserCreateOutputObjects()
  93 {
  94
  95
  96   // Create histograms
  97   // Called once
  98   const Int_t nRuns=200 ;
  99
 100   // ESD histograms
 101   if(fOutputContainer != NULL){
 102     delete fOutputContainer;
 103   }
 104   fOutputContainer = new THashList();
 105   fOutputContainer->SetOwner(kTRUE);
 106   PostData(1, fOutputContainer);
 107
 108   //========QA histograms=======
 109
 110   //Event selection
 111   fOutputContainer->Add(new TH2F("hSelEvents","Event selection", 10,0.,10.,nRuns,0.,float(nRuns))) ;
 112   fOutputContainer->Add(new TH1F("hTotSelEvents","Event selection", 10,0.,10.)) ;
 113
 114   //vertex distribution
 115   fOutputContainer->Add(new TH2F("hZvertex","Z vertex position", 50,-25.,25.,nRuns,0.,float(nRuns))) ;
 116
 117   //Centrality
 118   fOutputContainer->Add(new TH2F("hCentrality","Event centrality", 100,0.,100.,nRuns,0.,float(nRuns))) ;
 119   fOutputContainer->Add(new TH2F("hCenPHOS","Centrality vs PHOSclusters", 100,0.,100.,200,0.,200.)) ;
 120   fOutputContainer->Add(new TH2F("hCenPHOSCells","Centrality vs PHOS cells", 100,0.,100.,100,0.,1000.)) ;
 121   fOutputContainer->Add(new TH2F("hCenTrack","Centrality vs tracks", 100,0.,100.,100,0.,15000.)) ;
 122   fOutputContainer->Add(new TH2F("hCluEvsClu","ClusterMult vs E",200,0.,20.,100,0.,100.)) ;
 123   fOutputContainer->Add(new TH2F("hCluEvsCluM","ClusterMult vs E",200,0.,20.,100,0.,20.)) ;
 124
 125   fOutputContainer->Add(new TH3F("hCPVr","CPV radius",100,0.,20.,100,0.,20.,10,0.,100.));
 126
 127   //Bad Map
 128   fOutputContainer->Add(new TH2F("hCluLowM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 129   fOutputContainer->Add(new TH2F("hCluLowM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 130   fOutputContainer->Add(new TH2F("hCluLowM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 131
 132   fOutputContainer->Add(new TH2F("hCluHighM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 133   fOutputContainer->Add(new TH2F("hCluHighM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 134   fOutputContainer->Add(new TH2F("hCluHighM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 135
 136   fOutputContainer->Add(new TH2F("hCluVetoM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 137   fOutputContainer->Add(new TH2F("hCluVetoM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 138   fOutputContainer->Add(new TH2F("hCluVetoM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 139
 140   fOutputContainer->Add(new TH2F("hCluDispM1","Cell (X,Z), M1" ,64,0.5,64.5, 56,0.5,56.5));
 141   fOutputContainer->Add(new TH2F("hCluDispM2","Cell (X,Z), M2" ,64,0.5,64.5, 56,0.5,56.5));
 142   fOutputContainer->Add(new TH2F("hCluDispM3","Cell (X,Z), M3" ,64,0.5,64.5, 56,0.5,56.5));
 143
 144   //Single photon and pi0 spectrum
 145   Int_t nPtPhot = 300 ;
 146   Double_t ptPhotMax = 30 ;
 147   Int_t nM       = 500;
 148   Double_t mMin  = 0.0;
 149   Double_t mMax  = 1.0;
 150
 151   //PHOS calibration QA
 152   fOutputContainer->Add(new TH2F("hPi0M11","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 153   fOutputContainer->Add(new TH2F("hPi0M12","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 154   fOutputContainer->Add(new TH2F("hPi0M13","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 155   fOutputContainer->Add(new TH2F("hPi0M22","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 156   fOutputContainer->Add(new TH2F("hPi0M23","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 157   fOutputContainer->Add(new TH2F("hPi0M33","Pairs in modules",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 158
 159   char key[55] ;
 160   for(Int_t cent=0; cent<6; cent++){
 161
 162     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad2_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 163     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad4_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 164     fOutputContainer->Add(new TH1F(Form("hPhotAll_DistBad6_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 165
 166     fOutputContainer->Add(new TH1F(Form("hPhotAll_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 167     fOutputContainer->Add(new TH1F(Form("hPhotAllcore_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 168     fOutputContainer->Add(new TH1F(Form("hPhotAllwou_cen%d",cent),"All clusters",nPtPhot,0.,ptPhotMax));
 169     fOutputContainer->Add(new TH1F(Form("hPhotDisp_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 170     fOutputContainer->Add(new TH1F(Form("hPhotDisp2_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 171     fOutputContainer->Add(new TH1F(Form("hPhotDispcore_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 172     fOutputContainer->Add(new TH1F(Form("hPhotDisp2core_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 173     fOutputContainer->Add(new TH1F(Form("hPhotDispwou_cen%d",cent),"Disp clusters",nPtPhot,0.,ptPhotMax));
 174     fOutputContainer->Add(new TH1F(Form("hPhotCPV_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 175     fOutputContainer->Add(new TH1F(Form("hPhotCPVcore_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 176     fOutputContainer->Add(new TH1F(Form("hPhotCPV2_cen%d",cent),"CPV clusters",nPtPhot,0.,ptPhotMax));
 177     fOutputContainer->Add(new TH1F(Form("hPhotBoth_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
 178     fOutputContainer->Add(new TH1F(Form("hPhotBothcore_cen%d",cent),"Both clusters",nPtPhot,0.,ptPhotMax));
 179     fOutputContainer->Add(new TH1F(Form("hPhotBoth2_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
 180     fOutputContainer->Add(new TH1F(Form("hPhotBoth2core_cen%d",cent),"Both2 clusters",nPtPhot,0.,ptPhotMax));
 181
 182
 183     fOutputContainer->Add(new TH2F(Form("hPi0All_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 184     fOutputContainer->Add(new TH2F(Form("hPi0Allcore_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 185     fOutputContainer->Add(new TH2F(Form("hPi0Allwou_cen%d",cent),"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 186     fOutputContainer->Add(new TH2F(Form("hPi0Disp_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 187     fOutputContainer->Add(new TH2F(Form("hPi0Disp2_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 188     fOutputContainer->Add(new TH2F(Form("hPi0Dispcore_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 189     fOutputContainer->Add(new TH2F(Form("hPi0Disp2core_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 190     fOutputContainer->Add(new TH2F(Form("hPi0Dispwou_cen%d",cent),"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 191     fOutputContainer->Add(new TH2F(Form("hPi0CPV_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 192     fOutputContainer->Add(new TH2F(Form("hPi0CPVcore_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 193     fOutputContainer->Add(new TH2F(Form("hPi0CPV2_cen%d",cent),"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 194     fOutputContainer->Add(new TH2F(Form("hPi0Both_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 195     fOutputContainer->Add(new TH2F(Form("hPi0Bothcore_cen%d",cent),"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 196     fOutputContainer->Add(new TH2F(Form("hPi0Both2_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 197     fOutputContainer->Add(new TH2F(Form("hPi0Both2core_cen%d",cent),"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 198
 199
 200     snprintf(key,55,"hPi0All_a07_cen%d",cent) ;
 201     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 202     snprintf(key,55,"hPi0Disp_a07_cen%d",cent) ;
 203     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 204     snprintf(key,55,"hPi0CPV_a07_cen%d",cent) ;
 205     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 206     snprintf(key,55,"hPi0CPV2_a07_cen%d",cent) ;
 207     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 208     snprintf(key,55,"hPi0Both_a07_cen%d",cent) ;
 209     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 210     snprintf(key,55,"hPi0Both2_a07_cen%d",cent) ;
 211     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 212
 213     snprintf(key,55,"hSingleAll_cen%d",cent) ;
 214     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 215     snprintf(key,55,"hSingleAllcore_cen%d",cent) ;
 216     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 217     snprintf(key,55,"hSingleAllwou_cen%d",cent) ;
 218     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 219     snprintf(key,55,"hSingleDisp_cen%d",cent) ;
 220     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 221     snprintf(key,55,"hSingleDisp2_cen%d",cent) ;
 222     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 223     snprintf(key,55,"hSingleDispcore_cen%d",cent) ;
 224     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 225     snprintf(key,55,"hSingleDisp2core_cen%d",cent) ;
 226     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 227     snprintf(key,55,"hSingleDispwou_cen%d",cent) ;
 228     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 229     snprintf(key,55,"hSingleCPV_cen%d",cent) ;
 230     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 231     snprintf(key,55,"hSingleCPVcore_cen%d",cent) ;
 232     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 233     snprintf(key,55,"hSingleCPV2_cen%d",cent) ;
 234     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 235     snprintf(key,55,"hSingleBoth_cen%d",cent) ;
 236     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 237     snprintf(key,55,"hSingleBothcore_cen%d",cent) ;
 238     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 239     snprintf(key,55,"hSingleBoth2_cen%d",cent) ;
 240     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 241     snprintf(key,55,"hSingleBoth2core_cen%d",cent) ;
 242     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 243
 244
 245     snprintf(key,55,"hMiPi0All_cen%d",cent) ;
 246     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 247     snprintf(key,55,"hMiPi0Allcore_cen%d",cent) ;
 248     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 249     snprintf(key,55,"hMiPi0Allwou_cen%d",cent) ;
 250     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 251     snprintf(key,55,"hMiPi0Disp_cen%d",cent) ;
 252     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 253     snprintf(key,55,"hMiPi0Disp2_cen%d",cent) ;
 254     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 255     snprintf(key,55,"hMiPi0Dispwou_cen%d",cent) ;
 256     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 257     snprintf(key,55,"hMiPi0Dispcore_cen%d",cent) ;
 258     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 259     snprintf(key,55,"hMiPi0Disp2core_cen%d",cent) ;
 260     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 261     snprintf(key,55,"hMiPi0CPV_cen%d",cent) ;
 262     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 263     snprintf(key,55,"hMiPi0CPVcore_cen%d",cent) ;
 264     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 265     snprintf(key,55,"hMiPi0CPV2_cen%d",cent) ;
 266     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 267     snprintf(key,55,"hMiPi0Both_cen%d",cent) ;
 268     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 269     snprintf(key,55,"hMiPi0Bothcore_cen%d",cent) ;
 270     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 271     snprintf(key,55,"hMiPi0Both2_cen%d",cent) ;
 272     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 273     snprintf(key,55,"hMiPi0Both2core_cen%d",cent) ;
 274     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 275
 276     snprintf(key,55,"hMiPi0All_a07_cen%d",cent) ;
 277     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 278     snprintf(key,55,"hMiPi0Disp_a07_cen%d",cent) ;
 279     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 280     snprintf(key,55,"hMiPi0CPV_a07_cen%d",cent) ;
 281     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 282     snprintf(key,55,"hMiPi0CPV2_a07_cen%d",cent) ;
 283     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 284     snprintf(key,55,"hMiPi0Both_a07_cen%d",cent) ;
 285     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 286     snprintf(key,55,"hMiPi0Both2_a07_cen%d",cent) ;
 287     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 288
 289     snprintf(key,55,"hMiSingleAll_cen%d",cent) ;
 290     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 291     snprintf(key,55,"hMiSingleAllwou_cen%d",cent) ;
 292     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 293     snprintf(key,55,"hMiSingleAllcore_cen%d",cent) ;
 294     fOutputContainer->Add(new TH2F(key,"All clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 295     snprintf(key,55,"hMiSingleDisp_cen%d",cent) ;
 296     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 297     snprintf(key,55,"hMiSingleDisp2_cen%d",cent) ;
 298     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 299     snprintf(key,55,"hMiSingleDispwou_cen%d",cent) ;
 300     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 301     snprintf(key,55,"hMiSingleDispcore_cen%d",cent) ;
 302     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 303     snprintf(key,55,"hMiSingleDisp2core_cen%d",cent) ;
 304     fOutputContainer->Add(new TH2F(key,"Disp clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 305     snprintf(key,55,"hMiSingleCPV_cen%d",cent) ;
 306     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 307     snprintf(key,55,"hMiSingleCPVcore_cen%d",cent) ;
 308     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 309     snprintf(key,55,"hMiSingleCPV2_cen%d",cent) ;
 310     fOutputContainer->Add(new TH2F(key,"CPV clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 311     snprintf(key,55,"hMiSingleBoth_cen%d",cent) ;
 312     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 313     snprintf(key,55,"hMiSingleBothcore_cen%d",cent) ;
 314     fOutputContainer->Add(new TH2F(key,"Both clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 315     snprintf(key,55,"hMiSingleBoth2_cen%d",cent) ;
 316     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 317     snprintf(key,55,"hMiSingleBoth2core_cen%d",cent) ;
 318     fOutputContainer->Add(new TH2F(key,"Both2 clusters",nM,mMin,mMax,nPtPhot,0.,ptPhotMax));
 319   }
 320
 321
 322  //MC
 323   for(Int_t cent=0; cent<6; cent++){
 324     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
 325     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 326     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 327     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 328     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 329     fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
 330     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 331     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 332     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 333     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 334     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 335     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 336     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 337     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 338     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 339     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 340     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 341     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 342     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 343     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 344     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 345     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 346     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 347     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 348   }
 349   fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 350   fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 351   fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 352
 353
 354   Int_t nPt      = 200;
 355   Double_t ptMin = 0;
 356   Double_t ptMax = 20;
 357   fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
 358   fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
 359   fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
 360   fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
 361   fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
 362   fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
 363 //  fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
 364 //  fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
 365
 366   for(Int_t cen=0; cen<6; cen++){
 367     fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 368     fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 369     fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 370     fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 371     fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 372     fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 373     fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 374     fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 375     fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 376     fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 377     fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 378     fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 379     fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 380
 381     //pairs from common parents
 382     fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 383     fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 384     fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 385     fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 386     fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 387     fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 388     fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 389
 390     //common parent - pi0
 391     fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 392     fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 393     fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 394     fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 395     fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 396     fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 397     fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 398     fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 399     fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 400
 401   }
 402
 403
 404   //Photon contaminations
 405   fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 406   fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 407   fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 408   fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 409   fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
 410
 411    const Int_t nTypes=24 ;
 412    char partTypes[nTypes][55] ;
 413    snprintf(partTypes[0],55,"hGammaNoPrim") ; //
 414    snprintf(partTypes[1],55,"hGammaPhot") ; //
 415    snprintf(partTypes[2],55,"hGammaEl") ; //
 416    snprintf(partTypes[3],55,"hGammaPi0") ; //
 417    snprintf(partTypes[4],55,"hGammaEta") ; //
 418    snprintf(partTypes[5],55,"hhGammaOmega") ; //
 419    snprintf(partTypes[6],55,"hGammaPipm") ; //
 420    snprintf(partTypes[7],55,"hGammaP") ; //
 421    snprintf(partTypes[8],55,"hGammaPbar") ; //
 422    snprintf(partTypes[9],55,"hGammaN") ; //
 423    snprintf(partTypes[10],55,"hGammaNbar") ; //
 424    snprintf(partTypes[11],55,"hGammaK0S") ; //
 425    snprintf(partTypes[12],55,"hGammaK0L") ; //
 426    snprintf(partTypes[13],55,"hGammaKpm") ; //
 427    snprintf(partTypes[14],55,"hGammaKstar") ; //
 428    snprintf(partTypes[15],55,"hGammaDelta") ; //
 429    snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
 430    snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
 431    snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
 432    snprintf(partTypes[19],55,"hGammaPipmEl") ; //
 433    snprintf(partTypes[20],55,"hGammaPipmOther") ; //
 434    snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
 435    snprintf(partTypes[22],55,"hGammaPipmp") ; //
 436    snprintf(partTypes[23],55,"hGammaPipmn") ; //
 437
 438    const Int_t nPID=12 ;
 439    char cPID[12][25] ;
 440    snprintf(cPID[0],25,"All") ;
 441    snprintf(cPID[1],25,"Allcore") ;
 442    snprintf(cPID[2],25,"CPV") ;
 443    snprintf(cPID[3],25,"CPVcore") ;
 444    snprintf(cPID[4],25,"CPV2") ;
 445    snprintf(cPID[5],25,"CPV2core") ;
 446    snprintf(cPID[6],25,"Disp") ;
 447    snprintf(cPID[7],25,"Dispcore") ;
 448    snprintf(cPID[8],25,"Disp2") ;
 449    snprintf(cPID[9],25,"Disp2core") ;
 450    snprintf(cPID[10],25,"Both") ;
 451    snprintf(cPID[11],25,"Bothcore") ;
 452
 453    for(Int_t itype=0; itype<nTypes; itype++){
 454      for(Int_t iPID=0; iPID<nPID; iPID++){
 455        for(Int_t cen=0; cen<5; cen++){
 456          fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
 457        }
 458      }
 459    }
 460
 461
 462
 463   PostData(1, fOutputContainer);
 464
 465   printf("Init done \n") ;
 466
 467 }
 468
 469 //________________________________________________________________________
 470 void AliPHOSHijingEfficiency::UserExec(Option_t *)
 471 {
 472   // Main loop, called for each event
 473   // Analyze ESD/AOD
 474   fStack=0 ;
 475   if(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()){
 476     if(static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
 477       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
 478   }
 479
 480
 481   AliESDEvent *event = dynamic_cast<AliESDEvent*>(InputEvent());
 482   if (!event) {
 483     Printf("ERROR: Could not retrieve event");
 484     PostData(1, fOutputContainer);
 485     return;
 486   }
 487
 488   fRunNumber=ConvertRunNumber(event->GetRunNumber()) ;
 489   FillHistogram("hSelEvents",0.5,fRunNumber-0.5) ;
 490   FillHistogram("hTotSelEvents",0.5) ;
 491
 492   // Get PHOS rotation matrices from ESD and set them to the PHOS geometry
 493   if(fEventCounter == 0) {
 494     // Initialize the PHOS geometry
 495     fPHOSGeo = AliPHOSGeometry::GetInstance("IHEP") ;
 496
 497     //We have to apply re-calibration for pass1 LCH10h
 498     // Initialize decalibration factors in the form of the OCDB object
 499     const Float_t wDecalib = 0.055;
 500
 501     fPHOSCalibData = new AliPHOSCalibData();
 502     fPHOSCalibData->SetName("PhosCalibData");
 503
 504     for(Int_t module=1; module<=5; module++) {
 505       for(Int_t column=1; column<=56; column++) {
 506         for(Int_t row=1; row<=64; row++) {
 507           Float_t adcChannelEmc = gRandom->Gaus(1.,wDecalib);
 508           fPHOSCalibData->SetADCchannelEmc(module,column,row,adcChannelEmc);
 509         }
 510       }
 511     }
 512     for(Int_t mod=0; mod<5; mod++) {
 513       if(!event->GetPHOSMatrix(mod)) continue;
 514       fPHOSGeo->SetMisalMatrix(event->GetPHOSMatrix(mod),mod) ;
 515       Printf("PHOS geo matrix %p for module # %d is set\n", event->GetPHOSMatrix(mod), mod);
 516      }
 517
 518     fEventCounter++ ;
 519   }
 520
 521
 522   // Checks if we have a primary vertex
 523   // Get primary vertices form ESD
 524   const AliESDVertex *esdVertex5 = event->GetPrimaryVertex();
 525
 526   Double_t vtx0[3] = {0,0,0}; // don't rely on ESD vertex, assume (0,0,0)
 527   Double_t vtx5[3] ={0.,0.,0.};
 528
 529   vtx5[0] = esdVertex5->GetX();
 530   vtx5[1] = esdVertex5->GetY();
 531   vtx5[2] = esdVertex5->GetZ();
 532
 533
 534   FillHistogram("hZvertex",esdVertex5->GetZ(),fRunNumber-0.5);
 535   if (TMath::Abs(esdVertex5->GetZ()) > 10. ){
 536     PostData(1, fOutputContainer);
 537     return;
 538   }
 539   FillHistogram("hSelEvents",1.5,fRunNumber-0.5) ;
 540   FillHistogram("hTotSelEvents",1.5) ;
 541
 542 /*
 543   if(event->IsPileupFromSPD()){
 544     PostData(1, fOutputContainer);
 545     return;
 546   }
 547 */
 548   FillHistogram("hSelEvents",2.5,fRunNumber-0.5) ;
 549   FillHistogram("hTotSelEvents",2.5) ;
 550
 551
 552   /*
 553   //Vtx class z-bin
 554   Int_t zvtx = (Int_t)((vtx5[2]+10.)/2.) ;
 555   if(zvtx<0)zvtx=0 ;
 556   if(zvtx>9)zvtx=9 ;
 557   */
 558   //No dependence on zVtx observed, save memory
 559   Int_t zvtx=0 ;
 560
 561   AliCentrality *centrality = event->GetCentrality();
 562   fCentrality=centrality->GetCentralityPercentile("V0M");
 563
 564   if( fCentrality <= 0. || fCentrality>80. ){
 565     PostData(1, fOutputContainer);
 566     return;
 567   }
 568
 569   FillHistogram("hSelEvents",3.5,fRunNumber-0.5) ;
 570   FillHistogram("hTotSelEvents",3.5) ;
 571
 572   if(fCentrality<5.)
 573     fCenBin=0 ;
 574   else if(fCentrality<10.)
 575     fCenBin=1 ;
 576   else if(fCentrality<20.)
 577     fCenBin=2 ;
 578   else if(fCentrality<40.)
 579     fCenBin=3 ;
 580   else if(fCentrality<60.)
 581     fCenBin=4 ;
 582   else if(fCentrality<80.)
 583     fCenBin=5 ;
 584   else{
 585     PostData(1, fOutputContainer);
 586     return;
 587   }
 588
 589   fRP=0. ;
 590
 591   FillHistogram("hSelEvents",4.5,fRunNumber-0.5) ;
 592   FillHistogram("hTotSelEvents",4.5) ;
 593   //All event selections done
 594   FillHistogram("hCentrality",fCentrality,fRunNumber-0.5) ;
 595   //Reaction plane is defined in the range (0;pi)
 596   //We have 10 bins
 597   Double_t averageRP = 0.;//fRPV0A+fRPV0C+fRP ;
 598   Int_t irp=Int_t(10.*(averageRP)/TMath::Pi());
 599   if(irp>9)irp=9 ;
 600
 601   if(!fPHOSEvents[zvtx][fCenBin][irp])
 602     fPHOSEvents[zvtx][fCenBin][irp]=new TList() ;
 603   TList * prevPHOS = fPHOSEvents[zvtx][fCenBin][irp] ;
 604
 605   ProcessMC() ;
 606
 607   if(fPHOSEvent)
 608     fPHOSEvent->Clear() ;
 609   else
 610     fPHOSEvent = new TClonesArray("AliCaloPhoton",200) ;
 611
 612   //For re-calibration
 613   const Double_t logWeight=4.5 ;
 614
 615   TVector3 vertex(vtx5);
 616
 617   Int_t multClust = event->GetNumberOfCaloClusters();
 618   Int_t inPHOS=0,inMod1=0, inMod2=0, inMod3=0 ;
 619   Double_t avrgEm1=0,avrgEm2=0,avrgEm3=0; //average cluster energy
 620
 621   AliESDCaloCells * cells = event->GetPHOSCells() ;
 622   FillHistogram("hCenPHOSCells",fCentrality,cells->GetNumberOfCells()) ;
 623   FillHistogram("hCenTrack",fCentrality,event->GetNumberOfTracks()) ;
 624
 625
 626   TVector3 localPos ;
 627   for (Int_t i=0; i<multClust; i++) {
 628     AliESDCaloCluster *clu = event->GetCaloCluster(i);
 629     if ( !clu->IsPHOS() || clu->E()<0.3) continue;
 630
 631     Float_t  position[3];
 632     clu->GetPosition(position);
 633     TVector3 global(position) ;
 634     Int_t relId[4] ;
 635     fPHOSGeo->GlobalPos2RelId(global,relId) ;
 636     Int_t mod  = relId[0] ;
 637     Int_t cellX = relId[2];
 638     Int_t cellZ = relId[3] ;
 639     if ( !IsGoodChannel("PHOS",mod,cellX,cellZ) )
 640       continue ;
 641
 642
 643     FillHistogram("hCluEvsClu",clu->E(),clu->GetNCells()) ;
 644     FillHistogram("hCluEvsCluM",clu->E(),clu->GetM02()) ;
 645
 646     //Apply re-Calibreation
 647     AliPHOSEsdCluster cluPHOS1(*clu);
 648     cluPHOS1.Recalibrate(fPHOSCalibData,cells); // modify the cell energies
 649     Reclusterize(&cluPHOS1) ;
 650     cluPHOS1.EvalAll(logWeight,vertex);         // recalculate the cluster parameters
 651     cluPHOS1.SetE(rnlin(cluPHOS1.E()));// Users's nonlinearity
 652     if(cluPHOS1.E()<0.3) continue;
 653     if(clu->GetNCells()<3) continue ;
 654     if(clu->GetM02()<0.2)   continue ;
 655
 656
 657     if(mod==1){
 658       avrgEm1+=cluPHOS1.E() ;
 659       inMod1++ ;
 660     }
 661     if(mod==2){
 662       avrgEm2+=cluPHOS1.E() ;
 663       inMod2++ ;
 664     }
 665     if(mod==3){
 666       avrgEm3+=cluPHOS1.E() ;
 667       inMod3++ ;
 668     }
 669
 670     TLorentzVector pv1 ;
 671     cluPHOS1.GetMomentum(pv1 ,vtx0);
 672
 673     Double_t ecore=CoreEnergy(&cluPHOS1) ;
 674
 675     FillHistogram(Form("hCluLowM%d",mod),cellX,cellZ,1.);
 676     if(cluPHOS1.E()>1.5){
 677       FillHistogram(Form("hCluHighM%d",mod),cellX,cellZ,1.);
 678     }
 679
 680     if(inPHOS>=fPHOSEvent->GetSize()){
 681       fPHOSEvent->Expand(inPHOS+50) ;
 682     }
 683     new((*fPHOSEvent)[inPHOS]) AliCaloPhoton(pv1.X(),pv1.Py(),pv1.Z(),pv1.E()) ;
 684     AliCaloPhoton * ph = (AliCaloPhoton*)fPHOSEvent->At(inPHOS) ;
 685     ph->SetModule(mod) ;
 686     pv1*= ecore/pv1.E() ;
 687     ph->SetMomV2(&pv1) ;
 688     ph->SetNCells(clu->GetNCells());
 689     Double_t m02=0.,m20=0.;
 690     EvalLambdas(&cluPHOS1,0,m02, m20);
 691     ph->SetDispBit(TestLambda(clu->E(),m20,m02)) ;
 692     EvalLambdas(&cluPHOS1,1,m02, m20);
 693     ph->SetDisp2Bit(TestLambda2(clu->E(),m20,m02)) ;
 694
 695     Double_t distBC=clu->GetDistanceToBadChannel();
 696     if(distBC>2.){
 697       FillHistogram(Form("hPhotAll_DistBad2_cen%d",fCenBin),ph->Pt()) ;
 698       if(distBC>4.){
 699         FillHistogram(Form("hPhotAll_DistBad4_cen%d",fCenBin),ph->Pt()) ;
 700         if(distBC>6.)
 701           FillHistogram(Form("hPhotAll_DistBad6_cen%d",fCenBin),ph->Pt()) ;
 702       }
 703     }
 704     if(ph->IsDispOK()){
 705       FillHistogram(Form("hCluDispM%d",mod),cellX,cellZ,1.);
 706     }
 707     //Track matching
 708     Double_t dx=clu->GetTrackDx() ;
 709     Double_t dz=clu->GetTrackDz() ;
 710
 711     Bool_t cpvBit=kTRUE ; //No track matched by default
 712     Bool_t cpvBit2=kTRUE ; //More Strict criterion
 713     TArrayI * itracks = clu->GetTracksMatched() ;
 714     if(itracks->GetSize()>0){
 715       Int_t iTr = itracks->At(0);
 716       if(iTr>=0 && iTr<event->GetNumberOfTracks()){
 717         AliESDtrack *track = event->GetTrack(iTr) ;
 718         Double_t pt = track->Pt() ;
 719         Short_t charge = track->Charge() ;
 720         Double_t r=TestCPV(dx, dz, pt,charge) ;
 721         cpvBit=(r>2.) ;
 722         cpvBit2=(r>4.) ;
 723         FillHistogram("hCPVr",r,pv1.E(),fCentrality);
 724       }
 725     }
 726     ph->SetCPVBit(cpvBit) ;
 727     ph->SetCPV2Bit(cpvBit2) ;
 728     if(cpvBit){
 729       FillHistogram(Form("hCluVetoM%d",mod),cellX,cellZ,1.);
 730     }
 731     ph->SetEMCx(float(cellX)) ;
 732     ph->SetEMCz(float(cellZ)) ;
 733     ph->SetLambdas(cluPHOS1.GetM20(),cluPHOS1.GetM02()) ;
 734     ph->SetUnfolded(clu->GetNExMax()<2); // Remember, if it is unfolded
 735     Bool_t sure=  kTRUE;
 736     Int_t primary=FindPrimary(clu,sure) ;  //номер праймари частицы в стеке
 737     ph->SetPrimary(primary) ;
 738     ph->SetWeight(PrimaryWeight(primary)) ;
 739     inPHOS++ ;
 740   }
 741
 742   FillHistogram("hCenPHOS",fCentrality,inPHOS) ;
 743   if(inPHOS==0){
 744     PostData(1, fOutputContainer);
 745     fEventCounter++;
 746     return ;
 747   }
 748
 749
 750   for (Int_t i1=0; i1<inPHOS; i1++) {
 751     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 752     Double_t w=ph1->GetWeight() ;
 753     Double_t pt = ph1->Pt() ;
 754     Double_t ptv =ph1->GetMomV2()->Pt() ;
 755
 756     FillHistogram(Form("hPhotAll_cen%d",fCenBin),pt,w) ;
 757     FillHistogram(Form("hPhotAllcore_cen%d",fCenBin),ptv,w) ;
 758     if(ph1->IsntUnfolded()){
 759       FillHistogram(Form("hPhotAllwou_cen%d",fCenBin),pt,w) ;
 760     }
 761     if(ph1->IsCPVOK()){
 762       FillHistogram(Form("hPhotCPV_cen%d",fCenBin),pt,w) ;
 763       FillHistogram(Form("hPhotCPVcore_cen%d",fCenBin),ptv,w) ;
 764     }
 765     if(ph1->IsCPV2OK()){
 766       FillHistogram(Form("hPhotCPV2_cen%d",fCenBin),pt,w) ;
 767     }
 768     if(ph1->IsDispOK()){
 769       FillHistogram(Form("hPhotDisp_cen%d",fCenBin),pt,w) ;
 770
 771       FillHistogram(Form("hPhotDispcore_cen%d",fCenBin),ptv,w) ;
 772       if(ph1->IsntUnfolded()){
 773         FillHistogram(Form("hPhotDispwou_cen%d",fCenBin),pt,w) ;
 774       }
 775       if(ph1->IsCPVOK()){
 776         FillHistogram(Form("hPhotBoth_cen%d",fCenBin),pt,w) ;
 777         FillHistogram(Form("hPhotBothcore_cen%d",fCenBin),ptv,w) ;
 778       }
 779     }
 780     if(ph1->IsDisp2OK()){
 781       FillHistogram(Form("hPhotDisp2_cen%d",fCenBin),pt,w) ;
 782       FillHistogram(Form("hPhotDisp2core_cen%d",fCenBin),ptv,w) ;
 783       if(ph1->IsCPVOK()){
 784         FillHistogram(Form("hPhotBoth2_cen%d",fCenBin),pt,w) ;
 785         FillHistogram(Form("hPhotBoth2core_cen%d",fCenBin),ptv,w) ;
 786       }
 787     }
 788   }
 789
 790
 791   const Double_t alphaCut=0.7 ;
 792   //pi0
 793   for (Int_t i1=0; i1<inPHOS-1; i1++) {
 794     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 795     Double_t w1=ph1->GetWeight() ;
 796     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
 797       AliCaloPhoton * ph2=(AliCaloPhoton*)fPHOSEvent->At(i2) ;
 798       Double_t w2=ph2->GetWeight() ;
 799       Double_t w = w1*w2 ;
 800       if(FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary())!=-1)
 801         w=w1 ;
 802
 803       TLorentzVector p12  = *ph1  + *ph2;
 804       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 805
 806       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 807       Double_t m = p12.M() ;
 808       Double_t pt = p12.Pt() ;
 809       Double_t mv = p12.M() ;
 810       Double_t ptv = p12.Pt() ;
 811
 812       FillHistogram(Form("hPi0All_cen%d",fCenBin),m,pt,w) ;
 813       FillHistogram(Form("hPi0Allcore_cen%d",fCenBin),mv,ptv,w) ;
 814       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 815         FillHistogram(Form("hPi0Allwou_cen%d",fCenBin),m,pt,w) ;
 816       }
 817
 818
 819       Double_t pt1 = ph1->Pt() ;
 820       Double_t ptv1 = ph1->GetMomV2()->Pt() ;
 821       Double_t pt2 = ph2->Pt() ;
 822       Double_t ptv2 = ph2->GetMomV2()->Pt() ;
 823       FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt1,w) ;
 824       FillHistogram(Form("hSingleAll_cen%d",fCenBin),m,pt2,w) ;
 825       FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ;
 826       FillHistogram(Form("hSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ;
 827       if(ph1->IsntUnfolded())
 828         FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt1,w) ;
 829       if(ph2->IsntUnfolded())
 830         FillHistogram(Form("hSingleAllwou_cen%d",fCenBin),m,pt2,w) ;
 831       if(ph1->IsCPVOK()){
 832         FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt1,w) ;
 833         FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ;
 834       }
 835       if(ph2->IsCPVOK()){
 836         FillHistogram(Form("hSingleCPV_cen%d",fCenBin),m,pt2,w) ;
 837         FillHistogram(Form("hSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ;
 838       }
 839       if(ph1->IsCPV2OK()){
 840         FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt1,w) ;
 841       }
 842       if(ph2->IsCPV2OK()){
 843         FillHistogram(Form("hSingleCPV2_cen%d",fCenBin),m,pt2,w) ;
 844       }
 845       if(ph1->IsDispOK()){
 846         FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt1,w) ;
 847         if(ph1->IsntUnfolded()){
 848           FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt1,w) ;
 849         }
 850         FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ;
 851       }
 852       if(ph2->IsDispOK()){
 853         FillHistogram(Form("hSingleDisp_cen%d",fCenBin),m,pt2,w) ;
 854         if(ph1->IsntUnfolded()){
 855           FillHistogram(Form("hSingleDispwou_cen%d",fCenBin),m,pt2,w) ;
 856         }
 857         FillHistogram(Form("hSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ;
 858       }
 859       if(ph1->IsDisp2OK()){
 860         FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt1,w) ;
 861         FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ;
 862       }
 863       if(ph2->IsDisp2OK()){
 864         FillHistogram(Form("hSingleDisp2_cen%d",fCenBin),m,pt2,w) ;
 865         FillHistogram(Form("hSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ;
 866       }
 867       if(ph1->IsDispOK() && ph1->IsCPVOK()){
 868         FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt1,w) ;
 869         FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ;
 870       }
 871       if(ph2->IsDispOK() && ph2->IsCPVOK()){
 872         FillHistogram(Form("hSingleBoth_cen%d",fCenBin),m,pt2,w) ;
 873         FillHistogram(Form("hSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ;
 874       }
 875       if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
 876         FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt1,w) ;
 877         FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ;
 878       }
 879       if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
 880         FillHistogram(Form("hSingleBoth2_cen%d",fCenBin),m,pt2,w) ;
 881         FillHistogram(Form("hSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ;
 882       }
 883
 884
 885       if(a<alphaCut){
 886         FillHistogram(Form("hPi0All_a07_cen%d",fCenBin),m,pt,w) ;
 887       }
 888
 889       if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 CPV
 890         FillHistogram(Form("hPi0CPV_cen%d",fCenBin),m,pt,w) ;
 891         FillHistogram(Form("hPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ;
 892
 893         if(a<alphaCut){
 894           FillHistogram(Form("hPi0CPV_a07_cen%d",fCenBin),m,pt,w) ;
 895         }
 896       }
 897       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 898         FillHistogram(Form("hPi0CPV2_cen%d",fCenBin),p12.M(),p12.Pt(),w) ;
 899         if(a<alphaCut){
 900           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCenBin),m,pt,w) ;
 901         }
 902       }
 903       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 904         FillHistogram(Form("hPi0Disp2_cen%d",fCenBin),m,pt,w) ;
 905         FillHistogram(Form("hPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ;
 906         if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
 907           FillHistogram(Form("hPi0Both2_cen%d",fCenBin),m,pt,w) ;
 908           FillHistogram(Form("hPi0Both2core_cen%d",fCenBin),mv,ptv,w) ;
 909           if(a<alphaCut){
 910             FillHistogram(Form("hPi0Both2_a07_cen%d",fCenBin),m,pt,w) ;
 911           }
 912         }
 913       }
 914
 915
 916       if(ph1->IsDispOK() && ph2->IsDispOK()){ //pi0 Disp
 917         FillHistogram(Form("hPi0Disp_cen%d",fCenBin),m,pt,w) ;
 918         FillHistogram(Form("hPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ;
 919         if(a<alphaCut){
 920           FillHistogram(Form("hPi0Disp_a07_cen%d",fCenBin),m,pt,w) ;
 921         }
 922
 923         if(ph1->IsCPVOK() && ph2->IsCPVOK()){ //pi0 Both
 924           FillHistogram(Form("hPi0Both_cen%d",fCenBin),m,pt,w) ;
 925           FillHistogram(Form("hPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ;
 926
 927           if(a<alphaCut){
 928             FillHistogram(Form("hPi0Both_a07_cen%d",fCenBin),m,pt,w) ;
 929           }
 930         }
 931       }
 932     } // end of loop i2
 933   } // end of loop i1
 934
 935   //now mixed
 936   for (Int_t i1=0; i1<inPHOS; i1++) {
 937     AliCaloPhoton * ph1=(AliCaloPhoton*)fPHOSEvent->At(i1) ;
 938     Double_t w1=ph1->GetWeight() ;
 939     for(Int_t ev=0; ev<prevPHOS->GetSize();ev++){
 940       TClonesArray * mixPHOS = static_cast<TClonesArray*>(prevPHOS->At(ev)) ;
 941       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 942         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 943         Double_t w2=ph2->GetWeight() ;
 944         TLorentzVector p12  = *ph1  + *ph2;
 945         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 946
 947         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 948         Double_t m = p12.M() ;
 949         Double_t pt = p12.Pt() ;
 950         Double_t mv = p12.M() ;
 951         Double_t ptv = p12.Pt() ;
 952         Double_t w=w1*w2 ;
 953
 954         FillHistogram(Form("hMiPi0All_cen%d",fCenBin),m,pt,w) ;
 955         FillHistogram(Form("hMiPi0Allcore_cen%d",fCenBin),mv,ptv,w) ;
 956         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 957           FillHistogram(Form("hMiPi0Allwou_cen%d",fCenBin),m,pt,w) ;
 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,w) ;
 965         FillHistogram(Form("hMiSingleAll_cen%d",fCenBin),m,pt2,w) ;
 966         FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv1,w) ;
 967         FillHistogram(Form("hMiSingleAllcore_cen%d",fCenBin),mv,ptv2,w) ;
 968         if(ph1->IsntUnfolded())
 969           FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt1,w) ;
 970         if(ph2->IsntUnfolded())
 971           FillHistogram(Form("hMiSingleAllwou_cen%d",fCenBin),m,pt2,w) ;
 972         if(ph1->IsCPVOK()){
 973           FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt1,w) ;
 974           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv1,w) ;
 975         }
 976         if(ph2->IsCPVOK()){
 977           FillHistogram(Form("hMiSingleCPV_cen%d",fCenBin),m,pt2,w) ;
 978           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCenBin),mv,ptv2,w) ;
 979         }
 980         if(ph1->IsCPV2OK()){
 981           FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt1,w) ;
 982         }
 983         if(ph2->IsCPV2OK()){
 984           FillHistogram(Form("hMiSingleCPV2_cen%d",fCenBin),m,pt2,w) ;
 985         }
 986         if(ph1->IsDispOK()){
 987           FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt1,w) ;
 988           if(ph1->IsntUnfolded()){
 989             FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt1,w) ;
 990           }
 991           FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv1,w) ;
 992         }
 993         if(ph2->IsDispOK()){
 994           FillHistogram(Form("hMiSingleDisp_cen%d",fCenBin),m,pt2,w) ;
 995           if(ph2->IsntUnfolded()){
 996             FillHistogram(Form("hMiSingleDispwou_cen%d",fCenBin),m,pt2,w) ;
 997           }
 998           FillHistogram(Form("hMiSingleDispcore_cen%d",fCenBin),mv,ptv2,w) ;
 999         }
1000         if(ph1->IsDisp2OK()){
1001           FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt1,w) ;
1002           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv1,w) ;
1003         }
1004         if(ph2->IsDisp2OK()){
1005           FillHistogram(Form("hMiSingleDisp2_cen%d",fCenBin),m,pt2,w) ;
1006           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCenBin),mv,ptv2,w) ;
1007         }
1008         if(ph1->IsDispOK() && ph1->IsCPVOK()){
1009           FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt1,w) ;
1010           FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv1,w) ;
1011         }
1012         if(ph2->IsDispOK() && ph2->IsCPVOK()){
1013           FillHistogram(Form("hMiSingleBoth_cen%d",fCenBin),m,pt2,w) ;
1014           FillHistogram(Form("hMiSingleBothcore_cen%d",fCenBin),mv,ptv2,w) ;
1015         }
1016         if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
1017           FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt1,w) ;
1018           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv1,w) ;
1019         }
1020         if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
1021           FillHistogram(Form("hMiSingleBoth2_cen%d",fCenBin),m,pt2,w) ;
1022           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCenBin),mv,ptv2,w) ;
1023         }
1024
1025
1026
1027         if(a<alphaCut){
1028           FillHistogram(Form("hMiPi0All_a07_cen%d",fCenBin),m,pt,w) ;
1029         }
1030         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1031           FillHistogram(Form("hMiPi0CPV_cen%d",fCenBin),m,pt,w) ;
1032           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCenBin),mv,ptv,w) ;
1033
1034           if(a<alphaCut){
1035             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCenBin),m,pt,w) ;
1036           }
1037         }
1038         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
1039           FillHistogram(Form("hMiPi0CPV2_cen%d",fCenBin),m,pt,w) ;
1040
1041           if(a<alphaCut){
1042             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCenBin),m,pt,w) ;
1043           }
1044         }
1045         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
1046           FillHistogram(Form("hMiPi0Disp2_cen%d",fCenBin),m,pt,w) ;
1047           FillHistogram(Form("hMiPi0Disp2core_cen%d",fCenBin),mv,ptv,w) ;
1048           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1049             FillHistogram(Form("hMiPi0Both2_cen%d",fCenBin),m,pt,w) ;
1050             FillHistogram(Form("hMiPi0Both2core_cen%d",fCenBin),mv,ptv,w) ;
1051           }
1052         }
1053         if(ph1->IsDispOK() && ph2->IsDispOK()){
1054
1055           FillHistogram(Form("hMiPi0Disp_cen%d",fCenBin),m,pt,w) ;
1056           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCenBin),mv,ptv,w) ;
1057           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
1058             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCenBin),m,pt,w) ;
1059           }
1060
1061           if(a<alphaCut){
1062             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCenBin),m,pt,w) ;
1063           }
1064           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
1065             FillHistogram(Form("hMiPi0Both_cen%d",fCenBin),m,pt,w) ;
1066             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCenBin),mv,ptv,w) ;
1067
1068             if(a<alphaCut){
1069               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCenBin),m,pt,w) ;
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   if(!event) //not possible, but required by Coverity
1705     return 999.;
1706   Double_t mf = event->GetMagneticField(); //Positive for ++ and negative for --
1707
1708   if(mf<0.){ //field --
1709     meanZ = -0.468318 ;
1710     if(charge>0)
1711       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)) ;
1712     else
1713       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)) ;
1714   }
1715   else{ //Field ++
1716     meanZ= -0.468318;
1717     if(charge>0)
1718       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)) ;
1719     else
1720       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)) ;
1721   }
1722
1723   Double_t rz=(dz-meanZ)/sz ;
1724   Double_t rx=(dx-meanX)/sx ;
1725   return TMath::Sqrt(rx*rx+rz*rz) ;
1726 }
1727 //____________________________________________________________________________
1728 Bool_t AliPHOSHijingEfficiency::TestTOF(Double_t t, Double_t e){
1729
1730   Double_t sigma = TMath::Sqrt(2.23183e-09*2.23183e-09
1731                              +2.24611e-09*2.24611e-09/e
1732                              +5.65054e-09*5.65054e-09/e/e) ;
1733   sigma=TMath::Min(20.e-9,sigma) ; //for the soft (<400 MeV) part
1734   Double_t mean=1.1e-9 ;
1735   if(TMath::Abs(t-mean)<2.*sigma)
1736     return kTRUE ;
1737   else
1738     if(TMath::Abs(t-mean+100.e-9)<2.*sigma)
1739       return kTRUE ;
1740
1741   return kFALSE ;
1742 }
1743 //____________________________________________________________________________
1744 Double_t  AliPHOSHijingEfficiency::CoreEnergy(AliESDCaloCluster * clu){
1745   //calculate energy of the cluster in the circle with radius distanceCut around the maximum
1746
1747   //Can not use already calculated coordinates?
1748   //They have incidence correction...
1749   const Double_t distanceCut =3.5 ;
1750   const Double_t logWeight=4.5 ;
1751
1752   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1753 // Calculates the center of gravity in the local PHOS-module coordinates
1754   Float_t wtot = 0;
1755   Double_t xc[100]={0} ;
1756   Double_t zc[100]={0} ;
1757   Double_t x = 0 ;
1758   Double_t z = 0 ;
1759   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1760   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1761     Int_t relid[4] ;
1762     Float_t xi ;
1763     Float_t zi ;
1764     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1765     fPHOSGeo->RelPosInModule(relid, xi, zi);
1766     xc[iDigit]=xi ;
1767     zc[iDigit]=zi ;
1768     if (clu->E()>0 && elist[iDigit]>0) {
1769       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1770       x    += xc[iDigit] * w ;
1771       z    += zc[iDigit] * w ;
1772       wtot += w ;
1773     }
1774   }
1775   if (wtot>0) {
1776     x /= wtot ;
1777     z /= wtot ;
1778   }
1779   Double_t coreE=0. ;
1780   for(Int_t iDigit=0; iDigit < mulDigit; iDigit++) {
1781     Double_t distance = TMath::Sqrt((xc[iDigit]-x)*(xc[iDigit]-x)+(zc[iDigit]-z)*(zc[iDigit]-z)) ;
1782     if(distance < distanceCut)
1783       coreE += elist[iDigit] ;
1784   }
1785   //Apply non-linearity correction
1786   return rnlin(coreE) ;
1787 }
1788 //____________________________________________________________________________
1789 Bool_t  AliPHOSHijingEfficiency::AreNeibors(Int_t id1,Int_t id2){
1790
1791   Int_t relid1[4] ;
1792   fPHOSGeo->AbsToRelNumbering(id1, relid1) ;
1793
1794   Int_t relid2[4] ;
1795   fPHOSGeo->AbsToRelNumbering(id2, relid2) ;
1796
1797   if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
1798     Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
1799     Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
1800
1801     if (( coldiff <= 1 )  && ( rowdiff <= 1 )){   //At least common vertex
1802       //    if (( relid1[2]==relid2[2] && coldiff <= 1 )  || ( relid1[3]==relid2[3] &&  rowdiff <= 1 )){ //common side
1803       return 1 ;
1804     }
1805     else {
1806       if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
1807         return 0; //  Difference in row numbers is too large to look further
1808     }
1809     return 0 ;
1810
1811   }
1812
1813   return 0 ;
1814 }
1815 //____________________________________________________________________________
1816 void  AliPHOSHijingEfficiency::Reclusterize(AliESDCaloCluster * clu){
1817   //Re-clusterize to make continues cluster
1818
1819   const Int_t oldMulDigit=clu->GetNCells() ;
1820   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1821   UShort_t * dlist = clu->GetCellsAbsId();
1822
1823   Int_t index[oldMulDigit] ;
1824   Bool_t used[oldMulDigit] ;
1825   for(Int_t i=0; i<oldMulDigit; i++) used[i]=0 ;
1826   Int_t inClu=0 ;
1827   Double_t eMax=0. ;
1828   //find maximum
1829   for(Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
1830     if(eMax<elist[iDigit]){
1831       eMax=elist[iDigit];
1832       index[0]=iDigit ;
1833       inClu=1 ;
1834     }
1835   }
1836   if(inClu==0){ //empty cluster
1837     return ;
1838   }
1839   used[index[0]]=kTRUE ; //mark as used
1840   for(Int_t i=0; i<inClu; i++){
1841     for(Int_t iDigit=0 ;iDigit<oldMulDigit; iDigit++){
1842        if(used[iDigit]) //already used
1843          continue ;
1844        if(AreNeibors(dlist[index[i]],dlist[iDigit])){
1845          index[inClu]= iDigit ;
1846          inClu++ ;
1847          used[iDigit]=kTRUE ;
1848        }
1849     }
1850   }
1851
1852   if(inClu==oldMulDigit) //no need to modify
1853     return ;
1854
1855   clu->SetNCells(inClu);
1856   //copy
1857   UShort_t tmpD[oldMulDigit] ;
1858   Double_t tmpE[oldMulDigit] ;
1859   for(Int_t i=0; i<oldMulDigit; i++){
1860     tmpD[i]=dlist[i] ;
1861     tmpE[i]=elist[i] ;
1862   }
1863   //change order of digits in list so that
1864   //first inClu cells were true ones
1865   for(Int_t i=0; i<inClu; i++){
1866     dlist[i]=tmpD[index[i]] ;
1867     elist[i]=tmpE[index[i]] ;
1868   }
1869
1870
1871 }
1872
1873 //____________________________________________________________________________
1874 void  AliPHOSHijingEfficiency::EvalLambdas(AliESDCaloCluster * clu, Int_t iR,Double_t &m02, Double_t &m20){
1875   //calculate dispecrsion of the cluster in the circle with radius distanceCut around the maximum
1876
1877   Double_t rCut=0. ;
1878   if(iR==0)
1879     rCut=3.5 ;
1880   else
1881     rCut=4.5 ;
1882
1883
1884   Double32_t * elist = clu->GetCellsAmplitudeFraction() ;
1885 // Calculates the center of gravity in the local PHOS-module coordinates
1886   Float_t wtot = 0;
1887   Double_t xc[100]={0} ;
1888   Double_t zc[100]={0} ;
1889   Double_t x = 0 ;
1890   Double_t z = 0 ;
1891   Int_t mulDigit=TMath::Min(100,clu->GetNCells()) ;
1892   const Double_t logWeight=4.5 ;
1893   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1894     Int_t relid[4] ;
1895     Float_t xi ;
1896     Float_t zi ;
1897     fPHOSGeo->AbsToRelNumbering(clu->GetCellAbsId(iDigit), relid) ;
1898     fPHOSGeo->RelPosInModule(relid, xi, zi);
1899     xc[iDigit]=xi ;
1900     zc[iDigit]=zi ;
1901     if (clu->E()>0 && elist[iDigit]>0) {
1902       Float_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1903       x    += xc[iDigit] * w ;
1904       z    += zc[iDigit] * w ;
1905       wtot += w ;
1906     }
1907   }
1908   if (wtot>0) {
1909     x /= wtot ;
1910     z /= wtot ;
1911   }
1912
1913   wtot = 0. ;
1914   Double_t dxx  = 0.;
1915   Double_t dzz  = 0.;
1916   Double_t dxz  = 0.;
1917   Double_t xCut = 0. ;
1918   Double_t zCut = 0. ;
1919   for(Int_t iDigit=0; iDigit<mulDigit; iDigit++) {
1920     if (clu->E()>0 && elist[iDigit]>0.) {
1921         Double_t w = TMath::Max( 0., logWeight + TMath::Log( elist[iDigit] / clu->E() ) ) ;
1922         Double_t xi= xc[iDigit] ;
1923         Double_t zi= zc[iDigit] ;
1924         if((xi-x)*(xi-x)+(zi-z)*(zi-z) < rCut*rCut){
1925           xCut += w * xi ;
1926           zCut += w * zi ;
1927           dxx  += w * xi * xi ;
1928           dzz  += w * zi * zi ;
1929           dxz  += w * xi * zi ;
1930           wtot += w ;
1931         }
1932     }
1933
1934   }
1935   if (wtot>0) {
1936     xCut/= wtot ;
1937     zCut/= wtot ;
1938     dxx /= wtot ;
1939     dzz /= wtot ;
1940     dxz /= wtot ;
1941     dxx -= xCut * xCut ;
1942     dzz -= zCut * zCut ;
1943     dxz -= xCut * zCut ;
1944
1945     m02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1946     m20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
1947   }
1948   else {
1949     m20=m02=0.;
1950   }
1951
1952 }
1953 //___________________________________________________________________________
1954 void AliPHOSHijingEfficiency::ProcessMC(){
1955   //fill histograms for efficiensy etc. calculation
1956   const Double_t rcut = 1. ; //cut for primary particles
1957   //---------First pi0/eta-----------------------------
1958   char partName[10] ;
1959   char hkey[55] ;
1960
1961   if(!fStack) return ;
1962   for(Int_t i=0;i<fStack->GetNtrack();i++){
1963      TParticle* particle =  fStack->Particle(i);
1964     if(particle->GetPdgCode() == 111)
1965       snprintf(partName,10,"pi0") ;
1966     else
1967       if(particle->GetPdgCode() == 221)
1968         snprintf(partName,10,"eta") ;
1969       else
1970         if(particle->GetPdgCode() == 22)
1971            snprintf(partName,10,"gamma") ;
1972         else
1973            continue ;
1974
1975     //Primary particle
1976     Double_t r=particle->R() ;
1977     Double_t pt = particle->Pt() ;
1978     //Distribution over vertex
1979     FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
1980
1981     if(r >rcut)
1982       continue ;
1983
1984     //Total number of pi0 with creation radius <1 cm
1985     Double_t w = PrimaryParticleWeight(particle) ;
1986     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCenBin) ;
1987     FillHistogram(hkey,pt,w) ;
1988     if(TMath::Abs(particle->Y())<0.12){
1989       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCenBin) ;
1990       FillHistogram(hkey,pt,w) ;
1991     }
1992
1993     snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCenBin) ;
1994     FillHistogram(hkey,particle->Y(),w) ;
1995
1996     Double_t phi=particle->Phi() ;
1997     while(phi<0.)phi+=TMath::TwoPi() ;
1998     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1999     snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCenBin) ;
2000     FillHistogram(hkey,phi,w) ;
2001   }
2002 }
2003 //___________________________________________________________________________
2004 Int_t AliPHOSHijingEfficiency::FindPrimary(AliESDCaloCluster*clu,  Bool_t&sure){
2005   //Finds primary and estimates if it unique one?
2006   //First check can it be photon/electron
2007   const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
2008   Int_t n=clu->GetNLabels() ;
2009   for(Int_t i=0;  i<n;  i++){
2010     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
2011     Int_t pdg = p->GetPdgCode() ;
2012     if(pdg==22  ||  pdg==11 || pdg == -11){
2013       if(p->Energy()>emFraction*clu->E()){
2014         sure=kTRUE ;
2015         return clu->GetLabelAt(i);
2016       }
2017     }
2018   }
2019
2020   Double_t*  Ekin=  new  Double_t[n] ;
2021   for(Int_t i=0;  i<n;  i++){
2022     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
2023     Ekin[i]=p->P() ;  // estimate of kinetic energy
2024     if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
2025       Ekin[i]+=1.8  ;  //due to annihilation
2026     }
2027   }
2028   Int_t iMax=0;
2029   Double_t eMax=0.,eSubMax=0. ;
2030   for(Int_t i=0;  i<n;  i++){
2031     if(Ekin[i]>eMax){
2032       eSubMax=eMax;
2033       eMax=Ekin[i];
2034       iMax=i;
2035     }
2036   }
2037   if(eSubMax>0.8*eMax)//not obvious primary
2038     sure=kFALSE;
2039   else
2040     sure=kTRUE;
2041   delete[]  Ekin;
2042   return  clu->GetLabelAt(iMax) ;
2043 }
2044 //___________________________________________________________________________
2045 Double_t AliPHOSHijingEfficiency::PrimaryWeight(Int_t primary){
2046    //Check who is the primary and introduce weight to correct primary spectrum
2047
2048   if(primary<0 || primary>=fStack->GetNtrack())
2049     return 1 ;
2050   //trace primaries up to IP
2051   TParticle* particle =  fStack->Particle(primary);
2052   Double_t r=particle->R() ;
2053   Int_t mother = particle->GetFirstMother() ;
2054   while(mother>-1){
2055     if(r<1. && particle->GetPdgCode()==111)
2056       break ;
2057     particle =  fStack->Particle(mother);
2058     mother = particle->GetFirstMother() ;
2059     r=particle->R() ;
2060   }
2061
2062   return TMath::Max(0.,PrimaryParticleWeight(particle)) ;
2063 }
2064 //________________________________________________________________________
2065 Double_t AliPHOSHijingEfficiency::PrimaryParticleWeight(TParticle * particle){
2066
2067   Int_t pdg = particle->GetPdgCode() ;
2068   Int_t type=0 ;
2069   if(pdg == 111 || TMath::Abs(pdg)==211){
2070     type =1 ;
2071   }
2072   else{
2073     if(TMath::Abs(pdg)<1000){ //Kaon-like
2074       type =2 ;
2075     }
2076     else
2077       type = 3;  //baryons
2078   }
2079
2080   Double_t x = particle->Pt() ;
2081   if(type==1){
2082    if(fCenBin==0) //0-5
2083      return (1.662990+1.140890*x-0.192088*x*x)/(1.-0.806630*x+0.304771*x*x)+0.141690*x ;
2084    if(fCenBin==1) //5-10
2085      return (1.474351+0.791492*x-0.066369*x*x)/(1.-0.839338*x+0.317312*x*x)+0.093289*x ;
2086    if(fCenBin==2) //10-20
2087      return (1.174728+0.959681*x-0.137695*x*x)/(1.-0.788873*x+0.299538*x*x)+0.128759*x ;
2088    if(fCenBin==3) //20-40
2089      return (0.927335+0.475349*x+0.004364*x*x)/(1.-0.817966*x+0.309787*x*x)+0.086899*x ;
2090    if(fCenBin==4) //40-60
2091      return (0.676878+0.190680*x+0.077031*x*x)/(1.-0.790623*x+0.305183*x*x)+0.064510*x ;
2092    if(fCenBin==5) //60-80
2093      return (0.684726-0.606262*x+0.409963*x*x)/(1.-1.080061*x+0.456933*x*x)+0.005151*x ;
2094   }
2095   if(type==2){
2096    if(fCenBin==0) //0-5
2097      return (-0.417131+2.253936*x-0.337731*x*x)/(1.-0.909892*x+0.316820*x*x)+0.157312*x ;
2098    if(fCenBin==1) //5-10
2099      return (-0.352275+1.844466*x-0.248598*x*x)/(1.-0.897048*x+0.316462*x*x)+0.132461*x ;
2100    if(fCenBin==2) //10-20
2101      return (-0.475481+1.975108*x-0.336013*x*x)/(1.-0.801028*x+0.276705*x*x)+0.188164*x ;
2102    if(fCenBin==3) //20-40
2103      return (-0.198954+1.068789*x-0.103540*x*x)/(1.-0.848354*x+0.299209*x*x)+0.112939*x ;
2104    if(fCenBin==4) //40-60
2105      return (-0.111052+0.664041*x-0.019717*x*x)/(1.-0.804916*x+0.300779*x*x)+0.095784*x ;
2106    if(fCenBin==5) //0-5
2107      return (0.202788-0.439832*x+0.564585*x*x)/(1.-1.254029*x+0.679444*x*x)+0.016235*x ;
2108   }
2109   if(type==3){
2110    if(fCenBin==0) //0-5
2111      return (-1.312732+2.743568*x-0.375775*x*x)/(1.-0.717533*x+0.164694*x*x)+0.164445*x ;
2112    if(fCenBin==1) //5-10
2113      return (-1.229425+2.585889*x-0.330164*x*x)/(1.-0.715892*x+0.167386*x*x)+0.133085*x ;
2114    if(fCenBin==2) //10-20
2115      return (-1.135677+2.397489*x-0.320355*x*x)/(1.-0.709312*x+0.164350*x*x)+0.146095*x ;
2116    if(fCenBin==3) //20-40
2117      return (-0.889993+1.928263*x-0.220785*x*x)/(1.-0.715991*x+0.174729*x*x)+0.095098*x ;
2118    if(fCenBin==4) //40-60
2119      return (-0.539237+1.329118*x-0.115439*x*x)/(1.-0.722906*x+0.186832*x*x)+0.059267*x ;
2120    if(fCenBin==5) //60-80
2121      return (-0.518126+1.327628*x-0.130881*x*x)/(1.-0.665649*x+0.184300*x*x)+0.081701*x ;
2122   }
2123   return 1. ;
2124 }
2125 //________________________________________________________________________
2126 Int_t AliPHOSHijingEfficiency::FindCommonParent(Int_t iPart, Int_t jPart){
2127   //check if there is a common parent for particles i and j
2128   // -1: no common parent or wrong iPart/jPart
2129
2130   if(iPart==-1 || iPart>=fStack->GetNtrack() ||
2131      jPart==-1 || jPart>=fStack->GetNtrack()) return -1;
2132
2133   Int_t iprim1=iPart;
2134   while(iprim1>-1){
2135      Int_t iprim2=jPart;
2136      while(iprim2>-1){
2137        if(iprim1==iprim2)
2138          return iprim1 ;
2139        iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother();
2140      }
2141      iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother();
2142   }
2143   return -1;
2144 }
2145 //________________________________________________________________________
2146 Bool_t AliPHOSHijingEfficiency::HaveParent(Int_t iPart, Int_t pdgParent){
2147   //check if there is a common parent for particles i and j
2148   // -1: no common parent or wrong iPart/jPart
2149
2150   if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1;
2151
2152   Int_t iprim1=iPart;
2153   while(iprim1>-1){
2154     TParticle * tmp = fStack->Particle(iprim1) ;
2155     if(tmp->GetPdgCode()==pdgParent)
2156       return kTRUE ;
2157     iprim1=tmp->GetFirstMother();
2158   }
2159   return kFALSE;
2160 }
2161 //________________________________________________________________________
2162 Bool_t AliPHOSHijingEfficiency::InPi0mass(Double_t m, Double_t /*pt*/){
2163
2164  return TMath::Abs(m-0.135)<0.007*2.5 ;
2165 }