PWGGA/EMCALTasks/AliAnalysisTaskEMCALMesonGGSDMpPb.cxx

   1 #include "AliAnalysisTaskEMCALMesonGGSDMpPb.h"
   2
   3 // ROOT includes
   4 #include <vector>
   5 #include <Riostream.h>
   6 #include <TChain.h>
   7 #include <TTree.h>
   8 #include <TF1.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH3F.h>
  12 #include <TH1D.h>
  13 #include <TH2D.h>
  14 #include <TH3D.h>
  15 #include <TCanvas.h>
  16 #include <TList.h>
  17 #include <TFile.h>
  18 #include <TLorentzVector.h>
  19 #include <TNtuple.h>
  20 #include <TRandom3.h>
  21 #include <TGeoManager.h>
  22 #include <TGeoMatrix.h>
  23 #include <TGeoBBox.h>
  24 #include <TArrayI.h>
  25 #include <TArrayF.h>
  26 #include <TObjArray.h>
  27
  28 // STEER? includes
  29 #include "AliAnalysisTaskSE.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliVCluster.h"
  32 #include "AliVCaloCells.h"
  33 #include "AliLog.h"
  34 #include "AliPID.h"
  35 #include "AliStack.h"
  36 #include "AliESDtrack.h"
  37 #include "AliESDtrackCuts.h"
  38 #include "AliESDEvent.h"
  39 #include "AliAODEvent.h"
  40 #include "AliMCEvent.h"
  41 #include "AliInputEventHandler.h"
  42 #include "AliESDInputHandler.h"
  43 #include "AliAODInputHandler.h"
  44 #include "AliAODTrack.h"
  45 #include "AliExternalTrackParam.h"
  46 #include "AliESDfriendTrack.h"
  47 #include "AliTrackerBase.h"
  48
  49 // EMCAL includes
  50 #include "AliEMCALRecoUtils.h"
  51 #include "AliEMCALGeometry.h"
  52 #include "AliTrackerBase.h"
  53 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
  54 #include "AliEMCALPIDUtils.h"
  55 #include "AliExternalTrackParam.h"
  56
  57 #include "AliCentrality.h"
  58
  59 using std::cout;
  60 using std::endl;
  61
  62
  63 ClassImp(AliAnalysisTaskEMCALMesonGGSDMpPb)
  64
  65 //________________________________________________________________________
  66 AliAnalysisTaskEMCALMesonGGSDMpPb::AliAnalysisTaskEMCALMesonGGSDMpPb() :
  67   AliAnalysisTaskSE(),
  68   fOutput(0),
  69   fMcMode(0),
  70   fRecalibrator(0),
  71   fdRmin_ClustTrack(0),
  72   fPhimin(0),
  73   fPhimax(0),
  74   fEtamin(0),
  75   fEtamax(0),
  76   fTrackCuts(0),
  77   fEsdEv(0),
  78   fAodEv(0),
  79   h1_zvtx(0),
  80   h1_trigger(0),
  81   h1_centrality(0),
  82   h2_PhiEtaCluster(0),
  83   h2_PhiEtaClusterCut(0),
  84   h2_PhiEtaMaxCell(0),
  85   h2_PhiEtaMaxCellCut(0),
  86   h2_gE_RecTruth(0),
  87   h2_eop_E(0),
  88   h2_eop_pT(0),
  89   h2_E_time(0),
  90   h2_Pi0TruthPhiEta(0),
  91   h2_PriPi0TruthPhiEta(0),
  92   h2_Pi0TruthPhiEtaEmcal(0),
  93   h2_PriPi0TruthPhiEtaEmcal(0),
  94   h2_Pi0TruthPhiEta_Phi2piEta065(0),
  95   h2_Pi0TruthPhiEta_Phi2piEta1(0),
  96   h2_TruthPhotonsPhiEta(0),
  97   h2_PhotonsPhiEtaIsEmcal(0),
  98   TriggerList(0),
  99   fHelperClass(0)
 100 {
 101   // Dummy constructor ALWAYS needed for I/O.
 102   for(int i=0; i<cent_bins; i++){
 103     h1_nClusters[i] = 0;
 104     h1_M[i] = 0;
 105     h1_M_mix[i] = 0;
 106     h1_E[i] = 0;
 107     h1_dR_ClustTrk[i] = 0;
 108     h1_Pi0TruthPt[i] = 0;
 109     h1_PriPi0TruthPt[i] = 0;
 110     h1_Pi0TruthPtEmcal[i] = 0;
 111     h1_PriPi0TruthPtEmcal[i] = 0;
 112     h1_Pi0TruthPtPhi2piEta065[i] = 0;
 113     h1_Pi0TruthPtPhi2piEta1[i] = 0;
 114     h1_TruthPhotonsEmcal[i] = 0;
 115     h1_PhotonsEmcal[i] = 0;
 116     h1_PhotonsNCellsCut[i] = 0;
 117     h1_PhotonsTrackMatchCut[i] = 0;
 118     h1_PhotonsAllCut[i] = 0;
 119     h1_dR_RealMC[i] = 0;
 120     h1_Chi2[i] = 0;
 121     h1_nTrkMatch[i] = 0;
 122     h1_nCells[i] = 0;
 123     h1_ClusterDisp[i] = 0;
 124     h2_Ellipse[i] = 0;
 125     h2_EtaPt[i] = 0;
 126     h3_MptAsymm[i] = 0;
 127     h3_MptAsymm_mix[i] = 0;
 128     h2_dphi_deta[i] = 0;
 129     h2_dphi_deta_mix[i] = 0;
 130     h2_DispRes[i] = 0;
 131     h2_cells_M02[i] = 0;
 132   }
 133 }
 134
 135 //________________________________________________________________________
 136 AliAnalysisTaskEMCALMesonGGSDMpPb::AliAnalysisTaskEMCALMesonGGSDMpPb(const char *name) :
 137   AliAnalysisTaskSE(name),
 138   fOutput(0),
 139   fMcMode(0),
 140   fRecalibrator(0),
 141   fdRmin_ClustTrack(0),
 142   fPhimin(0),
 143   fPhimax(0),
 144   fEtamin(0),
 145   fEtamax(0),
 146   fTrackCuts(0),
 147   fEsdEv(0),
 148   fAodEv(0),
 149   h1_zvtx(0),
 150   h1_trigger(0),
 151   h1_centrality(0),
 152   h2_PhiEtaCluster(0),
 153   h2_PhiEtaClusterCut(0),
 154   h2_PhiEtaMaxCell(0),
 155   h2_PhiEtaMaxCellCut(0),
 156   h2_gE_RecTruth(0),
 157   h2_eop_E(0),
 158   h2_eop_pT(0),
 159   h2_E_time(0),
 160   h2_Pi0TruthPhiEta(0),
 161   h2_PriPi0TruthPhiEta(0),
 162   h2_Pi0TruthPhiEtaEmcal(0),
 163   h2_PriPi0TruthPhiEtaEmcal(0),
 164   h2_Pi0TruthPhiEta_Phi2piEta065(0),
 165   h2_Pi0TruthPhiEta_Phi2piEta1(0),
 166   h2_TruthPhotonsPhiEta(0),
 167   h2_PhotonsPhiEtaIsEmcal(0),
 168   TriggerList(0),
 169   fHelperClass(0)
 170 {
 171   // Constructor
 172   // Define input and output slots here (never in the dummy constructor)
 173   // Input slot #0 works with a TChain - it is connected to the default input container
 174   // Output slot #1 writes into a TH1 container
 175   DefineOutput(1, TList::Class());                                            // for output list
 176   for(int i=0; i<cent_bins; i++){
 177     h1_nClusters[i] = 0;
 178     h1_M[i] = 0;
 179     h1_M_mix[i] = 0;
 180     h1_E[i] = 0;
 181     h1_dR_ClustTrk[i] = 0;
 182     h1_Pi0TruthPt[i] = 0;
 183     h1_PriPi0TruthPt[i] = 0;
 184     h1_Pi0TruthPtEmcal[i] = 0;
 185     h1_PriPi0TruthPtEmcal[i] = 0;
 186     h1_Pi0TruthPtPhi2piEta065[i] = 0;
 187     h1_Pi0TruthPtPhi2piEta1[i] = 0;
 188     h1_TruthPhotonsEmcal[i] = 0;
 189     h1_PhotonsEmcal[i] = 0;
 190     h1_PhotonsNCellsCut[i] = 0;
 191     h1_PhotonsTrackMatchCut[i] = 0;
 192     h1_PhotonsAllCut[i] = 0;
 193     h1_dR_RealMC[i] = 0;
 194     h1_Chi2[i] = 0;
 195     h1_nTrkMatch[i] = 0;
 196     h1_nCells[i] = 0;
 197     h1_ClusterDisp[i] = 0;
 198     h2_Ellipse[i] = 0;
 199     h2_EtaPt[i] = 0;
 200     h3_MptAsymm[i] = 0;
 201     h3_MptAsymm_mix[i] = 0;
 202     h2_dphi_deta[i] = 0;
 203     h2_dphi_deta_mix[i] = 0;
 204     h2_DispRes[i] = 0;
 205     h2_cells_M02[i] = 0;
 206   }
 207 }
 208
 209 //________________________________________________________________________
 210 AliAnalysisTaskEMCALMesonGGSDMpPb::~AliAnalysisTaskEMCALMesonGGSDMpPb()
 211 {
 212   // Destructor. Clean-up the output list, but not the histograms that are put inside
 213   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 214   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 215     delete fOutput;
 216   }
 217   delete fTrackCuts;
 218 }
 219
 220 //________________________________________________________________________
 221 void AliAnalysisTaskEMCALMesonGGSDMpPb::UserCreateOutputObjects()
 222 {
 223   // Create histograms
 224   // Called once (on the worker node)
 225
 226   fOutput = new TList();
 227   fOutput->SetOwner();  // IMPORTANT!
 228
 229   fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 230
 231   cout << "__________AliAnalysisTaskEMCALMesonGGSDMpPb: Input settings__________" << endl;
 232   cout << " fMcMode:             " << fMcMode   << endl;
 233   cout << " fRecalibrator:       " << fRecalibrator << endl;
 234   cout << " dRmin_ClustTrack:    " << fdRmin_ClustTrack << endl;
 235   cout << " phi range:           " << fPhimin << ", " << fPhimax << endl;
 236   cout << " eta range:           " << fEtamin << ", " << fEtamax << endl;
 237   cout << " number of zvtx bins: " << zvtx_bins << endl;
 238   cout << " number of mult bins: " << mult_bins << endl;
 239   cout << " poolDepth:           " << poolDepth << endl;
 240   cout << endl;
 241
 242   char saythis1[500];
 243   char saythis2[500];
 244
 245   double TotalNBins = 0.0;
 246
 247   // Create histograms
 248   Int_t nClustersbins = 501;
 249   Float_t nClusterslow = -0.5, nClustersup = 500.5;
 250   for(int i=0; i<cent_bins; i++){
 251     sprintf(saythis1,"h1_nClusters_%d",i);
 252     sprintf(saythis2,"# of clusters");
 253     h1_nClusters[i] = new TH1F(saythis1, saythis2, nClustersbins, nClusterslow, nClustersup);
 254     h1_nClusters[i]->GetXaxis()->SetTitle("number of clusters/evt");
 255     h1_nClusters[i]->GetYaxis()->SetTitle("counts");
 256     h1_nClusters[i]->SetMarkerStyle(kFullCircle);
 257     TotalNBins+=nClustersbins;
 258   }
 259
 260   Int_t nZvertexbins = 501;
 261   Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
 262   h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
 263   h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
 264   h1_zvtx->GetYaxis()->SetTitle("counts");
 265   h1_zvtx->SetMarkerStyle(kFullCircle);
 266   TotalNBins+=nZvertexbins;
 267
 268   h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
 269   TotalNBins+=1001;
 270
 271   h1_centrality = new TH1F("h1_centrality", "centrality", 1001,-0.1,100.1);
 272   TotalNBins+=1001;
 273
 274   Int_t Mbins = 3000;
 275   Float_t Mlow = 0.0, Mup = 3.0;
 276   for(int i=0; i<cent_bins; i++){
 277     sprintf(saythis1,"h1_M_%d",i);
 278     sprintf(saythis2,"Invariant Mass");
 279     h1_M[i] = new TH1F(saythis1, saythis2, Mbins, Mlow, Mup);
 280     h1_M[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 281     h1_M[i]->GetYaxis()->SetTitle("counts");
 282     h1_M[i]->SetMarkerStyle(kFullCircle);
 283     TotalNBins+=Mbins;
 284   }
 285
 286   for(int i=0; i<cent_bins; i++){
 287     sprintf(saythis1,"h1_M_mix_%d",i);
 288     sprintf(saythis2,"Invariant Mass (mixed events)");
 289     h1_M_mix[i] = new TH1F(saythis1, saythis2, Mbins, Mlow, Mup);
 290     h1_M_mix[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 291     h1_M_mix[i]->GetYaxis()->SetTitle("counts");
 292     h1_M_mix[i]->SetMarkerStyle(kFullCircle);
 293     TotalNBins+=Mbins;
 294   }
 295
 296   Int_t ptbins = 2000;
 297   Float_t ptlow = 0.0, ptup = 20.0;
 298   Int_t Ebins = 500;
 299   Float_t Elow = 0.0, Eup = 20.0;
 300   for(int i=0; i<cent_bins; i++){
 301     sprintf(saythis1,"h1_E_%d",i);
 302     sprintf(saythis2,"Cluster Energy in EMCal");
 303     h1_E[i] = new TH1F(saythis1, saythis2, Ebins, Elow, Eup);
 304     h1_E[i]->GetXaxis()->SetTitle("E [GeV]");
 305     h1_E[i]->GetYaxis()->SetTitle("counts");
 306     h1_E[i]->SetMarkerStyle(kFullCircle);
 307     TotalNBins+=Ebins;
 308   }
 309
 310   h2_PhiEtaCluster = new TH2F("h2_PhiEtaCluster", "cluster phi vs eta", 400,1.362,3.178, 300,-0.728,0.728);
 311   h2_PhiEtaCluster->GetXaxis()->SetTitle("#phi [rad]");
 312   h2_PhiEtaCluster->GetYaxis()->SetTitle("#eta");
 313   h2_PhiEtaCluster->GetZaxis()->SetTitle("hits");
 314   h2_PhiEtaCluster->SetMarkerStyle(kFullCircle);
 315   TotalNBins+=400*300;
 316
 317   h2_PhiEtaClusterCut = new TH2F("h2_PhiEtaClusterCut", "cluster phi vs eta (after cuts)", 400,1.362,3.178, 300,-0.728,0.728);
 318   h2_PhiEtaClusterCut->GetXaxis()->SetTitle("#phi [rad]");
 319   h2_PhiEtaClusterCut->GetYaxis()->SetTitle("#eta");
 320   h2_PhiEtaClusterCut->GetZaxis()->SetTitle("hits");
 321   h2_PhiEtaClusterCut->SetMarkerStyle(kFullCircle);
 322   TotalNBins+=400*300;
 323
 324 // eta binning
 325   Double_t EtaBins[97] = {-0.66687,-0.653,-0.63913,-0.62526,-0.61139,-0.59752,-0.58365,-0.56978,-0.55591,-0.54204,-0.52817,-0.5143,-0.50043,-0.48656,-0.47269,-0.45882,-0.44495,-0.43108,-0.41721,-0.40334,-0.38947,-0.3756,-0.36173,-0.34786,-0.33399,-0.32012,-0.30625,-0.29238,-0.27851,-0.26464,-0.25077,-0.2369,-0.22303,-0.20916,-0.19529,-0.18142,-0.16755,-0.15368,-0.13981,-0.12594,-0.11207,-0.0982,-0.08433,-0.07046,-0.05659,-0.04272,-0.02885,-0.01498,-0.00111,0.01276,0.02663,0.0405,0.05437,0.06824,0.08211,0.09598,0.10985,0.12372,0.13759,0.15146,0.16533,0.1792,0.19307,0.20694,0.22081,0.23468,0.24855,0.26242,0.27629,0.29016,0.30403,0.3179,0.33177,0.34564,0.35951,0.37338,0.38725,0.40112,0.41499,0.42886,0.44273,0.4566,0.47047,0.48434,0.49821,0.51208,0.52595,0.53982,0.55369,0.56756,0.58143,0.5953,0.60917,0.62304,0.63691,0.65078,0.66465};
 326
 327   // phi binning
 328   Double_t PhiBins[125] = {1.408,1.4215,1.435,1.4485,1.462,1.4755,1.489,1.5025,1.516,1.5295,1.543,1.5565,1.57,1.5835,1.597,1.6105,1.624,1.6375,1.651,1.6645,1.678,1.6915,1.705,1.7185,1.732, 1.758,1.7715,1.785,1.7985,1.812,1.8255,1.839,1.8525,1.866,1.8795,1.893,1.9065,1.92,1.9335,1.947,1.9605,1.974,1.9875,2.001,2.0145,2.028,2.0415,2.055,2.0685,2.082,2.108,2.1215,2.135,2.1485,2.162,2.1755,2.189,2.2025,2.216,2.2295,2.243,2.2565,2.27,2.2835,2.297,2.3105,2.324,2.3375,2.351,2.3645,2.378,2.3915,2.405,2.4185,2.432,2.456,2.4695,2.483,2.4965,2.51,2.5235,2.537,2.5505,2.564,2.5775,2.591,2.6045,2.618,2.6315,2.645,2.6585,2.672,2.6855,2.699,2.7125,2.726,2.7395,2.753,2.7665,2.78,2.804,2.8175,2.831,2.8445,2.858,2.8715,2.885,2.8985,2.912,2.9255,2.939,2.9525,2.966,2.9795,2.993,3.0065,3.02,3.0335,3.047,3.0605,3.074,3.0875,3.101,3.1145,3.128};
 329
 330   h2_PhiEtaMaxCell = new TH2F("h2_PhiEtaMaxCell", "maxcell phi vs eta", 124,PhiBins, 96,EtaBins);
 331   h2_PhiEtaMaxCell->GetXaxis()->SetTitle("#phi [rad]");
 332   h2_PhiEtaMaxCell->GetYaxis()->SetTitle("#eta");
 333   h2_PhiEtaMaxCell->GetZaxis()->SetTitle("hits");
 334   h2_PhiEtaMaxCell->SetMarkerStyle(kFullCircle);
 335   TotalNBins+=96*124;
 336
 337   h2_PhiEtaMaxCellCut = new TH2F("h2_PhiEtaMaxCellCut", "maxcell phi vs eta (after cuts)", 124,PhiBins, 96,EtaBins);
 338   h2_PhiEtaMaxCellCut->GetXaxis()->SetTitle("#phi [rad]");
 339   h2_PhiEtaMaxCellCut->GetYaxis()->SetTitle("#eta");
 340   h2_PhiEtaMaxCellCut->GetZaxis()->SetTitle("hits");
 341   h2_PhiEtaMaxCellCut->SetMarkerStyle(kFullCircle);
 342   TotalNBins+=96*124;
 343
 344   for(int i=0; i<cent_bins; i++){
 345     sprintf(saythis1,"h1_dR_ClustTrk_%d",i);
 346     sprintf(saythis2,"Cluster-Track matching");
 347     h1_dR_ClustTrk[i] = new TH1F(saythis1, saythis2, 5000, -0.01, 5);
 348     h1_dR_ClustTrk[i]->GetXaxis()->SetTitle("dR [sqrt(d#phi^{2}+d#eta^{2})]");
 349     h1_dR_ClustTrk[i]->GetYaxis()->SetTitle("N");
 350     h1_dR_ClustTrk[i]->SetMarkerStyle(kFullCircle);
 351     TotalNBins+=5000;
 352   }
 353
 354   h2_gE_RecTruth = new TH2F("h2_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 500,0,2);
 355   h2_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
 356   h2_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
 357   h2_gE_RecTruth->GetZaxis()->SetTitle("counts");
 358   h2_gE_RecTruth->SetMarkerStyle(kFullCircle);
 359   TotalNBins+=Ebins*500;
 360
 361   h2_eop_E = new TH2F("h2_eop_E","E/p vs E (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
 362   h2_eop_E->GetXaxis()->SetTitle("cluster Energy [GeV]");
 363   h2_eop_E->GetYaxis()->SetTitle("E/p");
 364   TotalNBins+=Ebins*1200;
 365
 366   h2_eop_pT = new TH2F("h2_eop_pT","E/p vs p_{T} (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
 367   h2_eop_pT->GetXaxis()->SetTitle("cluster Energy [GeV]");
 368   h2_eop_pT->GetYaxis()->SetTitle("E/p");
 369   TotalNBins+=Ebins*1200;
 370
 371   h2_E_time = new TH2F("h2_E_time","cluster energy vs time", Ebins, Elow, Eup, 1000,-1e-6,1e-6);
 372   h2_E_time->GetXaxis()->SetTitle("cluster Energy [GeV]");
 373   h2_E_time->GetYaxis()->SetTitle("time [s]");
 374   TotalNBins+=Ebins*1000;
 375
 376   for(int i=0; i<cent_bins; i++){
 377     sprintf(saythis1,"h1_Pi0TruthPt_%d",i);
 378     sprintf(saythis2,"P_{T} distribution for Truth Pi0's");
 379     h1_Pi0TruthPt[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 380     h1_Pi0TruthPt[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 381     h1_Pi0TruthPt[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 382     h1_Pi0TruthPt[i]->SetMarkerStyle(kFullCircle);
 383     TotalNBins+=ptbins;
 384   }
 385
 386   for(int i=0; i<cent_bins; i++){
 387     sprintf(saythis1,"h1_PriPi0TruthPt_%d",i);
 388     sprintf(saythis2,"P_{T} distribution for Truth Primary Pi0's");
 389     h1_PriPi0TruthPt[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 390     h1_PriPi0TruthPt[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 391     h1_PriPi0TruthPt[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 392     h1_PriPi0TruthPt[i]->SetMarkerStyle(kFullCircle);
 393     TotalNBins+=ptbins;
 394   }
 395
 396   for(int i=0; i<cent_bins; i++){
 397     sprintf(saythis1,"h1_Pi0TruthPtEmcal_%d",i);
 398     sprintf(saythis2,"P_{T} distribution for Truth Pi0's (hit EMCal)");
 399     h1_Pi0TruthPtEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 400     h1_Pi0TruthPtEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 401     h1_Pi0TruthPtEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 402     h1_Pi0TruthPtEmcal[i]->SetMarkerStyle(kFullCircle);
 403     TotalNBins+=ptbins;
 404   }
 405
 406   for(int i=0; i<cent_bins; i++){
 407     sprintf(saythis1,"h1_PriPi0TruthPtEmcal_%d",i);
 408     sprintf(saythis2,"P_{T} distribution for Truth Primary Pi0's (hit EMCal)");
 409     h1_PriPi0TruthPtEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 410     h1_PriPi0TruthPtEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 411     h1_PriPi0TruthPtEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 412     h1_PriPi0TruthPtEmcal[i]->SetMarkerStyle(kFullCircle);
 413     TotalNBins+=ptbins;
 414   }
 415
 416   for(int i=0; i<cent_bins; i++){
 417     sprintf(saythis1,"h1_Pi0TruthPtPhi2piEta065_%d",i);
 418     sprintf(saythis2,"P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]");
 419     h1_Pi0TruthPtPhi2piEta065[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 420     h1_Pi0TruthPtPhi2piEta065[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 421     h1_Pi0TruthPtPhi2piEta065[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 422     h1_Pi0TruthPtPhi2piEta065[i]->SetMarkerStyle(kFullCircle);
 423     TotalNBins+=ptbins;
 424   }
 425
 426   for(int i=0; i<cent_bins; i++){
 427     sprintf(saythis1,"h1_Pi0TruthPtPhi2piEta1_%d",i);
 428     sprintf(saythis2,"P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]");
 429     h1_Pi0TruthPtPhi2piEta1[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 430     h1_Pi0TruthPtPhi2piEta1[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 431     h1_Pi0TruthPtPhi2piEta1[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 432     h1_Pi0TruthPtPhi2piEta1[i]->SetMarkerStyle(kFullCircle);
 433     TotalNBins+=ptbins;
 434   }
 435
 436   h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 437   h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 438   h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 439   TotalNBins+=380*200;
 440
 441   h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 442   h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 443   h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 444   TotalNBins+=380*200;
 445
 446   h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-1.5,1.5);
 447   h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 448   h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 449   TotalNBins+=380*150;
 450
 451   h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
 452   h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 453   h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 454   TotalNBins+=380*150;
 455
 456   h2_Pi0TruthPhiEta_Phi2piEta065 = new TH2F("h2_Pi0TruthPhiEta_Phi2piEta065",
 457                                             "Pi0Truth Phi vs Eta [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", 380,-0.02,6.30, 150,-5,5);
 458   h2_Pi0TruthPhiEta_Phi2piEta065->GetXaxis()->SetTitle("#phi [rad]");
 459   h2_Pi0TruthPhiEta_Phi2piEta065->GetYaxis()->SetTitle("#eta ");
 460   TotalNBins+=380*150;
 461
 462   h2_Pi0TruthPhiEta_Phi2piEta1 = new TH2F("h2_Pi0TruthPhiEta_Phi2piEta1",
 463                                             "Pi0Truth Phi vs Eta [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", 380,-0.02,6.30, 150,-1.5,1.5);
 464   h2_Pi0TruthPhiEta_Phi2piEta1->GetXaxis()->SetTitle("#phi [rad]");
 465   h2_Pi0TruthPhiEta_Phi2piEta1->GetYaxis()->SetTitle("#eta ");
 466   TotalNBins+=380*150;
 467
 468   for(int i=0; i<cent_bins; i++){
 469     sprintf(saythis1,"h1_TruthPhotonsEmcal_%d",i);
 470     sprintf(saythis2,"P_{T} distribution for photons (in EMCal)");
 471     h1_TruthPhotonsEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 472     h1_TruthPhotonsEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 473     h1_TruthPhotonsEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 474     h1_TruthPhotonsEmcal[i]->SetMarkerStyle(kFullCircle);
 475     TotalNBins+=ptbins;
 476   }
 477
 478   h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
 479                                    "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
 480   h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 481   h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
 482   h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
 483   TotalNBins+=380*150;
 484
 485   for(int i=0; i<cent_bins; i++){
 486     sprintf(saythis1,"h1_PhotonsEmcal_%d",i);
 487     sprintf(saythis2,"P_{T} distribution for photons (in EMCal)");
 488     h1_PhotonsEmcal[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 489     h1_PhotonsEmcal[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 490     h1_PhotonsEmcal[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 491     h1_PhotonsEmcal[i]->SetMarkerStyle(kFullCircle);
 492     TotalNBins+=ptbins;
 493   }
 494
 495   for(int i=0; i<cent_bins; i++){
 496     sprintf(saythis1,"h1_PhotonsNCellsCut_%d",i);
 497     sprintf(saythis2,"P_{T} distribution for #gamma's that survive NCells cut");
 498     h1_PhotonsNCellsCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 499     h1_PhotonsNCellsCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 500     h1_PhotonsNCellsCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 501     h1_PhotonsNCellsCut[i]->SetMarkerStyle(kFullCircle);
 502     TotalNBins+=ptbins;
 503   }
 504
 505   for(int i=0; i<cent_bins; i++){
 506     sprintf(saythis1,"h1_PhotonsTrackMatchCut_%d",i);
 507     sprintf(saythis2,"P_{T} distribution for #gamma's that survive TrackMatch cut");
 508     h1_PhotonsTrackMatchCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 509     h1_PhotonsTrackMatchCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 510     h1_PhotonsTrackMatchCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 511     h1_PhotonsTrackMatchCut[i]->SetMarkerStyle(kFullCircle);
 512     TotalNBins+=ptbins;
 513   }
 514
 515   for(int i=0; i<cent_bins; i++){
 516     sprintf(saythis1,"h1_PhotonsAllCut_%d",i);
 517     sprintf(saythis2,"P_{T} distribution for #gamma's that survive All cut");
 518     h1_PhotonsAllCut[i] = new TH1F(saythis1, saythis2, ptbins, ptlow, ptup);
 519     h1_PhotonsAllCut[i]->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 520     h1_PhotonsAllCut[i]->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 521     h1_PhotonsAllCut[i]->SetMarkerStyle(kFullCircle);
 522     TotalNBins+=ptbins;
 523   }
 524
 525   h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
 526                                      "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 100,-1.0,1.0);
 527   h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
 528   h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
 529   TotalNBins+=380*100;
 530
 531   for(int i=0; i<cent_bins; i++){
 532     sprintf(saythis1,"h1_dR_RealMC_%d",i);
 533     sprintf(saythis2,"P_{T} distribution for #gamma's that survive All cut");
 534     h1_dR_RealMC[i] = new TH1F(saythis1, saythis2, 2000, -0.01, 10);
 535     h1_dR_RealMC[i]->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
 536     h1_dR_RealMC[i]->GetYaxis()->SetTitle("N");
 537     h1_dR_RealMC[i]->SetMarkerStyle(kFullCircle);
 538     TotalNBins+=2000;
 539   }
 540
 541   Int_t chi2bins = 100;
 542   Float_t chi2low = -2, chi2up = 2;
 543   for(int i=0; i<cent_bins; i++){
 544     sprintf(saythis1,"h1_Chi2_%d",i);
 545     sprintf(saythis2,"#chi^{2} distribution for reconstructed");
 546     h1_Chi2[i] = new TH1F(saythis1,saythis2,chi2bins, chi2low, chi2up);
 547     h1_Chi2[i]->GetXaxis()->SetTitle("#chi^{2}");
 548     h1_Chi2[i]->GetYaxis()->SetTitle("counts");
 549     TotalNBins+=chi2bins;
 550   }
 551
 552   for(int i=0; i<cent_bins; i++){
 553     sprintf(saythis1,"h1_nTrkMatch_%d",i);
 554     sprintf(saythis2,"number of matched tracks");
 555     h1_nTrkMatch[i] = new TH1F(saythis1,saythis2,14, -1.5, 5.5);
 556     h1_nTrkMatch[i]->GetXaxis()->SetTitle("nTracksMatched");
 557     h1_nTrkMatch[i]->GetYaxis()->SetTitle("counts");
 558     TotalNBins+=14;
 559   }
 560
 561   for(int i=0; i<cent_bins; i++){
 562     sprintf(saythis1,"h1_ClusterDisp_%d",i);
 563     sprintf(saythis2,"Dispersion of CaloCluster");
 564     h1_ClusterDisp[i] = new TH1F(saythis1,saythis2,1000, -1, 3);
 565     h1_ClusterDisp[i]->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
 566     h1_ClusterDisp[i]->GetYaxis()->SetTitle("counts");
 567     TotalNBins+=1000;
 568   }
 569
 570   for(int i=0; i<cent_bins; i++){
 571     sprintf(saythis1,"h2_Ellipse_%d",i);
 572     sprintf(saythis2,"Ellipse axis M20 vs M02");
 573     h2_Ellipse[i] = new TH2F(saythis1,saythis2,500, -0.01, 1, 500, -0.01, 1);
 574     h2_Ellipse[i]->GetXaxis()->SetTitle("cluster->GetM20()");
 575     h2_Ellipse[i]->GetYaxis()->SetTitle("cluster->GetM02()");
 576     h2_Ellipse[i]->GetZaxis()->SetTitle("counts");
 577     TotalNBins+=500*500;
 578   }
 579
 580   Int_t etabins = 150;
 581   Float_t etalow = -1.5, etaup = 1.5;
 582   for(int i=0; i<cent_bins; i++){
 583     sprintf(saythis1,"h2_EtaPt_%d",i);
 584     sprintf(saythis2,"Cluster Energy vs ");
 585     h2_EtaPt[i] = new TH2F(saythis1,saythis2,etabins, etalow, etaup, ptbins, ptlow, ptup);
 586     h2_EtaPt[i]->GetXaxis()->SetTitle("E [GeV]");
 587     h2_EtaPt[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 588     TotalNBins+=etabins*ptbins;
 589   }
 590
 591   for(int i=0; i<cent_bins; i++){
 592     sprintf(saythis1,"h3_MptAsymm_%d",i);
 593     sprintf(saythis2,"mass vs p_{T} vs Asymm cut");
 594     h3_MptAsymm[i] = new TH3F(saythis1,saythis2,Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
 595     h3_MptAsymm[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 596     h3_MptAsymm[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 597     h3_MptAsymm[i]->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
 598     TotalNBins+=Mbins*ptbins*3.0;
 599   }
 600
 601   for(int i=0; i<cent_bins; i++){
 602     sprintf(saythis1,"h3_MptAsymm_mix_%d",i);
 603     sprintf(saythis2,"mass vs p_{T} vs Asymm cut (mixed events)");
 604     h3_MptAsymm_mix[i] = new TH3F(saythis1,saythis2,Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
 605     h3_MptAsymm_mix[i]->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 606     h3_MptAsymm_mix[i]->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 607     h3_MptAsymm_mix[i]->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
 608     TotalNBins+=Mbins*ptbins*3.0;
 609   }
 610
 611   for(int i=0; i<cent_bins; i++){
 612     sprintf(saythis1,"h2_dphi_deta_%d",i);
 613     sprintf(saythis2,"#Delta#phi vs #Delta#eta");
 614     h2_dphi_deta[i] = new TH2F(saythis1,saythis2, 349,-1.5,5, 400,-2.0,2.0);
 615     h2_dphi_deta[i]->GetXaxis()->SetTitle("#Delta#phi");
 616     h2_dphi_deta[i]->GetYaxis()->SetTitle("#Delta#eta");
 617     TotalNBins+=349*400;
 618   }
 619
 620   for(int i=0; i<cent_bins; i++){
 621     sprintf(saythis1,"h2_dphi_deta_mix_%d",i);
 622     sprintf(saythis2,"#Delta#phi vs #Delta#eta (mixed events)");
 623     h2_dphi_deta_mix[i] = new TH2F(saythis1,saythis2, 349,-1.5,5, 400,-2.0,2.0);
 624     h2_dphi_deta_mix[i]->GetXaxis()->SetTitle("#Delta#phi");
 625     h2_dphi_deta_mix[i]->GetYaxis()->SetTitle("#Delta#eta");
 626     TotalNBins+=349*400;
 627   }
 628
 629   for(int i=0; i<cent_bins; i++){
 630     sprintf(saythis1,"h2_DispRes_%d",i);
 631     sprintf(saythis2,"zvtx info");
 632     h2_DispRes[i] = new TH2F(saythis1, saythis2, 500,-0.01,1, 500,-0.1,2);
 633     h2_DispRes[i]->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
 634     h2_DispRes[i]->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
 635     h2_DispRes[i]->GetZaxis()->SetTitle("counts");
 636     TotalNBins+=500*500;
 637   }
 638
 639   for(int i=0; i<cent_bins; i++){
 640     sprintf(saythis1,"h2_cells_M02_%d",i);
 641     sprintf(saythis2,"nCells vs M02");
 642     h2_cells_M02[i] = new TH2F(saythis1, saythis2, 204,-1.5,100.5, 500,-1,1.5);
 643     h2_cells_M02[i]->GetXaxis()->SetTitle("nCells");
 644     h2_cells_M02[i]->GetYaxis()->SetTitle("M02");
 645     h2_cells_M02[i]->GetZaxis()->SetTitle("counts");
 646     TotalNBins+=204*500;
 647   }
 648
 649   cout << endl << "Total number of bins in booked histograms:  " << TotalNBins << endl << endl;
 650
 651   // Initialize helper class (for vertex selection & pile up correction)
 652   fHelperClass = new AliAnalysisUtils();
 653
 654   //TFile *f = OpenFile(1);
 655   //TDirectory::TContext context(f);
 656
 657   fOutput->Add(h1_zvtx);
 658   fOutput->Add(h1_trigger);
 659   fOutput->Add(h1_centrality);
 660   fOutput->Add(h2_PhiEtaCluster);
 661   fOutput->Add(h2_PhiEtaClusterCut);
 662   fOutput->Add(h2_PhiEtaMaxCell);
 663   fOutput->Add(h2_PhiEtaMaxCellCut);
 664   fOutput->Add(h2_gE_RecTruth);
 665   fOutput->Add(h2_eop_E);
 666   fOutput->Add(h2_eop_pT);
 667   fOutput->Add(h2_E_time);
 668
 669   for(int i=0; i<cent_bins; i++){
 670     fOutput->Add(h1_nClusters[i]);
 671     fOutput->Add(h1_M[i]);
 672     fOutput->Add(h1_M_mix[i]);
 673     fOutput->Add(h1_E[i]);
 674     fOutput->Add(h1_dR_ClustTrk[i]);
 675     fOutput->Add(h1_Pi0TruthPt[i]);
 676     fOutput->Add(h1_PriPi0TruthPt[i]);
 677     fOutput->Add(h1_Pi0TruthPtEmcal[i]);
 678     fOutput->Add(h1_PriPi0TruthPtEmcal[i]);
 679     fOutput->Add(h1_Pi0TruthPtPhi2piEta065[i]);
 680     fOutput->Add(h1_Pi0TruthPtPhi2piEta1[i]);
 681     fOutput->Add(h1_TruthPhotonsEmcal[i]);
 682     fOutput->Add(h1_PhotonsEmcal[i]);
 683     fOutput->Add(h1_PhotonsNCellsCut[i]);
 684     fOutput->Add(h1_PhotonsTrackMatchCut[i]);
 685     fOutput->Add(h1_PhotonsAllCut[i]);
 686     fOutput->Add(h1_dR_RealMC[i]);
 687     fOutput->Add(h1_Chi2[i]);
 688     fOutput->Add(h1_nTrkMatch[i]);
 689     fOutput->Add(h1_ClusterDisp[i]);
 690     fOutput->Add(h2_Ellipse[i]);
 691     fOutput->Add(h2_EtaPt[i]);
 692     fOutput->Add(h3_MptAsymm[i]);
 693     fOutput->Add(h3_MptAsymm_mix[i]);
 694     fOutput->Add(h2_dphi_deta[i]);
 695     fOutput->Add(h2_dphi_deta_mix[i]);
 696     fOutput->Add(h2_DispRes[i]);
 697     fOutput->Add(h2_cells_M02[i]);
 698   }
 699   fOutput->Add(h2_Pi0TruthPhiEta);
 700   fOutput->Add(h2_PriPi0TruthPhiEta);
 701   fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
 702   fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
 703   fOutput->Add(h2_Pi0TruthPhiEta_Phi2piEta065);
 704   fOutput->Add(h2_Pi0TruthPhiEta_Phi2piEta1);
 705   fOutput->Add(h2_TruthPhotonsPhiEta);
 706   fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
 707
 708   // Post data for ALL output slots >0 here,
 709   // To get at least an empty histogram
 710   // 1 is the outputnumber of a certain weg of task 1
 711   PostData(1, fOutput);
 712 }
 713
 714 //________________________________________________________________________
 715 void AliAnalysisTaskEMCALMesonGGSDMpPb::UserExec(Option_t *)
 716 {
 717   // Main loop Called for each event
 718
 719   AliMCEvent *mcEvent = MCEvent();
 720   Bool_t isMC = bool(mcEvent);//is this the right way to do this?
 721
 722   TRandom3 randy; randy.SetSeed(0);
 723   unsigned int iskip = -1;
 724   TLorentzVector ParentMix;
 725
 726   double recalScale = 1.0;
 727
 728   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 729
 730   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
 731   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
 732   if (!aodH && !esdH)  Printf("ERROR: Could not get ESD or AODInputHandler");
 733
 734   if(esdH)      fEsdEv = esdH->GetEvent();
 735   else if(aodH) fAodEv = aodH->GetEvent();
 736   else{
 737     AliFatal("Neither ESD nor AOD event found");
 738     return;
 739   }
 740
 741   // get pointer to reconstructed event
 742   AliVEvent *event = InputEvent();
 743   if (!event){
 744     AliError("Pointer == 0, this can not happen!");  return;}
 745   //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
 746   //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 747   //if (!fEsdEv){
 748   //AliError("Cannot get the ESD event");  return;}
 749
 750   fHelperClass->SetCutOnZVertexSPD(kFALSE);//does the zvtx have to match the spd vertex?
 751   fHelperClass->SetMaxVtxZ(1.0e6);//i set this myself later..
 752   // simply makes sure that there is at least 1 contributer to the zvtx determination.
 753   // this should only remove the *extra* events at zvtx==0.
 754   if(!fHelperClass->IsVertexSelected2013pA(event))
 755     return;
 756
 757   Int_t iTrigger = 0;
 758   if (fEsdEv)       iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
 759   else if (fAodEv)  iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
 760
 761   char saythis[500];
 762   Int_t iTriggerBin = 0;
 763   for(unsigned long j=0; j<TriggerList.size(); j++){
 764     if(iTrigger==TriggerList[j])
 765       iTriggerBin=j+1;
 766   }
 767   if(iTriggerBin==0){
 768     TriggerList.push_back(iTrigger);
 769     iTriggerBin=TriggerList.size();
 770   }
 771
 772   h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
 773   sprintf(saythis,"%d",iTrigger);
 774   h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
 775
 776   Double_t centralityVZERO=0.0;
 777   Int_t centBin = 0;
 778   AliCentrality *aliCent=NULL;
 779
 780   Int_t nclusters=0;
 781   if(fEsdEv){
 782     //Int_t evtN      = fEsdEv->GetEventNumberInFile();
 783     //Int_t ntracks   = fEsdEv->GetNumberOfTracks();
 784     nclusters = fEsdEv->GetNumberOfCaloClusters();
 785     aliCent   = fEsdEv->GetCentrality();
 786   }
 787   else if(fAodEv){
 788     //Int_t evtN      = fAodEv->GetEventNumberInFile();
 789     //Int_t ntracks   = fAodEv->GetNumberOfTracks();
 790     nclusters = fAodEv->GetNumberOfCaloClusters();
 791     aliCent   = fAodEv->GetCentrality();
 792   }
 793   //centBin = aliCent->GetCentralityClass10("V0M");
 794   //centralityVZERO = aliCent->GetCentralityPercentile("V0M");
 795   //centBin = aliCent->GetCentralityClass10("V0C");
 796   //centralityVZERO = aliCent->GetCentralityPercentile("V0C");
 797   centBin = aliCent->GetCentralityClass10("V0A");
 798   centralityVZERO = aliCent->GetCentralityPercentile("V0A");
 799
 800   if     (centralityVZERO<20.0)
 801     centBin = 0;
 802   else if(centralityVZERO<40.0)
 803     centBin = 1;
 804   else if(centralityVZERO<60.0)
 805     centBin = 2;
 806   else
 807     centBin = 3;
 808
 809   //cout << "Centrality:  " << centBin << "    " << centralityVZERO << endl;
 810
 811   if (fEsdEv){
 812     if(!(fEsdEv->GetPrimaryVertex()->GetStatus()))   return;
 813   }
 814   //else if (fAodEv){
 815   //if(!(fAodEv->GetPrimaryVertex()->GetStatus()))   return;
 816   //}
 817
 818   Double_t vertDisp=0.0;
 819   Double_t vertZres=0.0;
 820   Bool_t vertIsfromZ=0;
 821   if (fEsdEv){
 822     vertDisp    = fEsdEv->GetPrimaryVertex()->GetDispersion();
 823     vertZres    = fEsdEv->GetPrimaryVertex()->GetZRes();
 824     vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
 825   }
 826   else if (fAodEv){
 827     vertDisp    = 0;
 828     vertZres    = 0;
 829     vertIsfromZ = 0;
 830   }
 831
 832   h2_DispRes[centBin]->Fill(vertDisp, vertZres);
 833   // if vertex is from spd vertexZ, require more stringent cut
 834   if (vertIsfromZ) {
 835     if (vertDisp>0.02 ||  vertZres>0.25 )
 836       return; // bad vertex from VertexerZ
 837   }
 838
 839   // EMCal cluster loop for reconstructed event
 840   //numberofclusters set above!
 841   TLorentzVector Photon1, Photon2, Parent;
 842   Double_t vertex[3];
 843   Double_t E1=0.0;
 844   Double_t vertZ=0.0;
 845   if (fEsdEv)       vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
 846   else if (fAodEv)  vertZ = fAodEv->GetPrimaryVertex()->GetZ();
 847
 848   h1_zvtx->Fill(vertZ);
 849   //zvertex cut:
 850   if(fabs(vertZ)>10.0)
 851     return;
 852
 853   h1_nClusters[centBin]->Fill(nclusters);
 854   h1_centrality->Fill(centralityVZERO);
 855
 856   int izvtx = GetZvtxBin(vertZ);
 857   int imult = GetMultBin(nclusters);
 858
 859   //cout << iskip << " " << izvtx << " " << imult << endl;
 860   //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
 861
 862
 863
 864   //######################### ~~~~~~~~~~~ ##################################
 865   //######################### STARTING MC ##################################
 866   //######################### ~~~~~~~~~~~ ##################################
 867
 868   if(isMC){
 869     int isPrimary = 0;
 870
 871     if (!mcEvent){
 872       cout << "no MC event" << endl;
 873       return;
 874     }
 875
 876     const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
 877     if (!evtVtx)
 878       return;
 879
 880     mcEvent->PreReadAll();
 881
 882     Int_t nTracksMC  = mcEvent->GetNumberOfTracks();
 883     Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
 884     //cout << "We have  " << nPTracksMC << "  primaries of  " << nTracksMC << "  total tracks." << endl;
 885
 886     for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
 887       AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
 888       if (!mcP)
 889         continue;
 890
 891
 892       // it's a pion !!
 893       if(mcP->PdgCode() != 111)
 894         continue;
 895
 896       /*
 897       // primary particle
 898       Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
 899                                 (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
 900       if(dR <= 0.01)  isPrimary = 1;
 901       else            isPrimary = 0;
 902       */
 903
 904       if(iTrack<nPTracksMC)  isPrimary = 1;
 905       else                   isPrimary = 0;
 906
 907       h1_Pi0TruthPt    [centBin]->Fill(mcP->Pt());
 908       h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
 909
 910       if(isPrimary==1){
 911         h1_PriPi0TruthPt    [centBin]->Fill(mcP->Pt());
 912         h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
 913       }
 914
 915       if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
 916         continue;
 917
 918       h1_Pi0TruthPtPhi2piEta1    [centBin]->Fill(mcP->Pt());
 919       h2_Pi0TruthPhiEta_Phi2piEta1->Fill(mcP->Phi(),mcP->Eta());
 920
 921       if(mcP->Eta()>fEtamin && mcP->Eta()<fEtamax){
 922         h1_Pi0TruthPtPhi2piEta065    [centBin]->Fill(mcP->Pt());
 923         h2_Pi0TruthPhiEta_Phi2piEta065->Fill(mcP->Phi(),mcP->Eta());
 924       }
 925
 926
 927       Int_t d1 = mcP->GetFirstDaughter();
 928       Int_t d2 = mcP->GetLastDaughter();
 929
 930       if (d1<0)  continue;
 931       if (d2<0)  d2=d1;
 932       if (d2-d1 != 1)  continue;
 933
 934       bool bacc = true;
 935       bool binp = true;
 936       for (Int_t i=d1;i<=d2;++i){
 937         const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(i));
 938         Double_t eta_d = dmc->Eta();
 939         Double_t phi_d = dmc->Phi();
 940         if(!(dmc->PdgCode()==22)){
 941           binp = false;
 942         }
 943         if(!(dmc->PdgCode()==22 && eta_d>fEtamin && eta_d<fEtamax && phi_d>fPhimin && phi_d<fPhimax)){
 944           bacc = false;
 945         }
 946       }
 947
 948       if(binp && bacc){// 2 Photons hit the EMCAL!
 949
 950         for (Int_t j=d1;j<=d2;++j){//both truth photons.
 951
 952           const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(j));
 953           Double_t eta_d = dmc->Eta();
 954           Double_t phi_d = dmc->Phi();
 955
 956           if( dmc->PdgCode()==22 &&
 957               dmc->Eta()>fEtamin && dmc->Eta()<fEtamax &&
 958               dmc->Phi()>fPhimin && dmc->Phi()<fPhimax ){
 959             h1_TruthPhotonsEmcal[centBin]->Fill(dmc->Pt());
 960             h2_TruthPhotonsPhiEta->Fill(dmc->Phi(),dmc->Eta());
 961           }
 962
 963           for(int i=0; i<nclusters; i++) {
 964
 965             Bool_t matches_pion_photon = 0;
 966
 967             AliESDCaloCluster* esdCluster=NULL;
 968             AliAODCaloCluster* aodCluster=NULL;
 969             if (fEsdEv)       esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
 970             else if (fAodEv)  aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
 971
 972             Double_t clustMC_phi, clustMC_eta;
 973
 974             if(fEsdEv){
 975
 976               if(esdCluster->IsEMCAL()){
 977
 978                 Float_t pos[3] = {0,0,0};
 979                 esdCluster->GetPosition(pos);
 980                 TVector3 vpos(pos);
 981                 //h1_Phi->Fill(vpos.Phi());
 982                 clustMC_phi = vpos.Phi();
 983                 clustMC_eta = vpos.Eta();
 984
 985                 Double_t dR = TMath::Sqrt((eta_d-clustMC_eta)*(eta_d-clustMC_eta) +
 986                                           (phi_d-clustMC_phi)*(phi_d-clustMC_phi));
 987                 h1_dR_RealMC[centBin]->Fill(dR);
 988                 if(dR<=0.04) matches_pion_photon = 1;
 989
 990                 vpos.Delete();
 991               }
 992               if(matches_pion_photon){
 993                 if(esdCluster->IsEMCAL()){
 994                   h1_PhotonsEmcal[centBin]->Fill(esdCluster->E());
 995                   h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
 996                 }
 997                 if(esdCluster->IsEMCAL() && esdCluster->GetNCells()>=2)
 998                   h1_PhotonsNCellsCut[centBin]->Fill(esdCluster->E());
 999                 if(esdCluster->IsEMCAL() && esdCluster->GetNTracksMatched()==0)
1000                   h1_PhotonsTrackMatchCut[centBin]->Fill(esdCluster->E());
1001                 if(esdCluster->IsEMCAL() && esdCluster->GetNCells()>=2 && esdCluster->GetNTracksMatched()==0)
1002                   h1_PhotonsAllCut[centBin]->Fill(esdCluster->E());
1003               }//if(matches_pion_photon)
1004
1005             }//if(fEsdEv)
1006             else if(fAodEv){
1007
1008               if(aodCluster->IsEMCAL()){
1009
1010                 Float_t pos[3] = {0,0,0};
1011                 aodCluster->GetPosition(pos);
1012                 TVector3 vpos(pos);
1013                 //h1_Phi->Fill(vpos.Phi());
1014                 clustMC_phi = vpos.Phi();
1015                 clustMC_eta = vpos.Eta();
1016
1017                 Double_t dR = TMath::Sqrt((eta_d-clustMC_eta)*(eta_d-clustMC_eta) +
1018                                           (phi_d-clustMC_phi)*(phi_d-clustMC_phi));
1019                 h1_dR_RealMC[centBin]->Fill(dR);
1020                 if(dR<=0.04) matches_pion_photon = 1;
1021
1022                 vpos.Delete();
1023               }
1024               if(matches_pion_photon){
1025                 if(aodCluster->IsEMCAL()){
1026                   h1_PhotonsEmcal[centBin]->Fill(aodCluster->E());
1027                   h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1028                 }
1029                 if(aodCluster->IsEMCAL() && aodCluster->GetNCells()>=2)
1030                   h1_PhotonsNCellsCut[centBin]->Fill(aodCluster->E());
1031                 if(aodCluster->IsEMCAL() && aodCluster->GetNTracksMatched()==0)
1032                   h1_PhotonsTrackMatchCut[centBin]->Fill(aodCluster->E());
1033                 if(aodCluster->IsEMCAL() && aodCluster->GetNCells()>=2 && aodCluster->GetNTracksMatched()==0)
1034                   h1_PhotonsAllCut[centBin]->Fill(aodCluster->E());
1035
1036               }//if(matches_pion_photon)
1037
1038             }//if(fAodEv)
1039
1040           }//loop over nclusters.
1041
1042         }//both truth photons.
1043
1044       }// 2 Photons hit the EMCAL!
1045
1046
1047       if(binp && bacc){// 2 Photons hit the EMCAL!
1048         h1_Pi0TruthPtEmcal    [centBin]->Fill(mcP->Pt());
1049         h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1050
1051         if(isPrimary==1){
1052           h1_PriPi0TruthPtEmcal    [centBin]->Fill(mcP->Pt());
1053           h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1054         }
1055
1056       }//2 photons hit the EMCAL!
1057
1058     }//for(nTracksMC)
1059
1060   }//if(isMC)
1061
1062   //######################### ~~~~~~~~~~~~ ##################################
1063   //######################### DONE WITH MC ##################################
1064   //######################### ~~~~~~~~~~~~ ##################################
1065
1066
1067   for(int i=0; i<nclusters; i++) {
1068
1069     AliESDCaloCluster* esdCluster=NULL;
1070     AliAODCaloCluster* aodCluster=NULL;
1071     if (fEsdEv)       esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1072     else if (fAodEv)  aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1073     if(!esdCluster && !aodCluster) {
1074       AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i));
1075       continue;
1076     }
1077
1078     if(fEsdEv){
1079
1080       recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1081
1082       //uncomment this to do the track matching (1 of 3 lines, esd part)!!
1083       //Bool_t MatchesToTrack = 0;
1084       if(esdCluster->IsEMCAL()){
1085
1086         Float_t pos[3] = {0,0,0};
1087         Short_t maxCellID = -1;
1088         Float_t celleta, cellphi;
1089         esdCluster->GetPosition(pos);
1090         TVector3 clusterPosition(pos);
1091         h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1092         GetMaxCellEnergy(esdCluster, maxCellID);
1093         AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1094         fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1095         h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1096
1097         // _______________Track loop for reconstructed event_____________
1098         for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) {
1099           AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track
1100           if(!esdTrack) {
1101             AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk));
1102             continue;
1103           }
1104
1105           Double_t posTrk[3] = {0,0,0};
1106           esdTrack->GetXYZ(posTrk);
1107           TVector3 vposTrk(posTrk);
1108
1109           Double_t fMass          = 0.139;
1110           Double_t fStepSurface   = 20.;
1111           Float_t etaproj, phiproj, pttrackproj;
1112
1113           AliExternalTrackParam *trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1114           if(!trackParam) continue;
1115           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1116
1117           double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1118                                     (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1119
1120           h1_dR_ClustTrk[centBin]->Fill(dR_clusttrk);
1121
1122           //uncomment this to do the track matching (2 of 3 lines)!!
1123           //if(dR_clusttrk<fdRmin_ClustTrack)
1124           //MatchesToTrack = 1;
1125
1126         }//_____________________________nTracks__________________________
1127
1128         h2_cells_M02  [centBin]->Fill(esdCluster->GetNCells(),esdCluster->GetM02());
1129         h2_Ellipse    [centBin]->Fill(esdCluster->GetM20(),esdCluster->GetM02());
1130         h1_Chi2       [centBin]->Fill(esdCluster->Chi2());//always -1.
1131         h1_nTrkMatch  [centBin]->Fill(esdCluster->GetNTracksMatched());
1132         h1_ClusterDisp[centBin]->Fill(esdCluster->GetDispersion());
1133         h2_E_time              ->Fill(esdCluster->E(),esdCluster->GetTOF());
1134
1135         TArrayI *TrackLabels = esdCluster->GetTracksMatched();
1136         if(TrackLabels){
1137           if(TrackLabels->GetSize()>0){
1138             Int_t trackindex = TrackLabels->At(0);
1139             AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track
1140
1141             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1142             h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P());
1143             h2_eop_pT->Fill(matchingT->Pt(),            esdCluster->E()*recalScale/matchingT->P());
1144           }
1145         }
1146
1147         //uncomment this to do the track matching (3 of 3 lines)!!
1148         //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){
1149         if(isGoodEsdCluster(esdCluster)){
1150           recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1151           E1 = esdCluster->E()*recalScale;// TOTAL HACK - JJ
1152           fEsdEv->GetVertex()->GetXYZ(vertex);
1153           esdCluster->GetMomentum(Photon1,vertex);
1154           Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1155           Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1156           Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1157           Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1158           h1_E[centBin]->Fill(E1);
1159           h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1160           h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1161         }
1162         clusterPosition.Delete();
1163       }//if(esdCluster->isEMCAL())
1164     }//if(fEsdEv)
1165     else if(fAodEv){
1166
1167       recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1168
1169       //uncomment this to do the track matching (1 of 3 lines, aod part)!!
1170       //Bool_t MatchesToTrack = 0;
1171       if(aodCluster->IsEMCAL()){
1172
1173         Float_t pos[3] = {0,0,0};
1174         Short_t maxCellID = -1;
1175         Float_t celleta, cellphi;
1176         aodCluster->GetPosition(pos);
1177         TVector3 clusterPosition(pos);
1178         h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1179         GetMaxCellEnergy(aodCluster, maxCellID);
1180         AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1181         fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1182         h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1183
1184         // _______________Track loop for reconstructed event_____________
1185         for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) {
1186           AliAODTrack* aodTrack = fAodEv->GetTrack(itrk); // pointer to reconstructed to track
1187           if(!aodTrack) {
1188             AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk));
1189             continue;
1190           }
1191
1192           Double_t posTrk[3] = {0,0,0};
1193           aodTrack->GetXYZ(posTrk);
1194           TVector3 vposTrk(posTrk);
1195
1196           Double_t fMass          = 0.139;
1197           Double_t fStepSurface   = 20.;
1198           Float_t etaproj, phiproj, pttrackproj;
1199
1200           AliExternalTrackParam *trackParam =  const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam());
1201           if(!trackParam) continue;
1202           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1203
1204           double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1205                                     (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1206
1207           h1_dR_ClustTrk[centBin]->Fill(dR_clusttrk);
1208
1209           //uncomment this to do the track matching (2 of 3 lines, aod part)!!
1210           //if(dR_clusttrk<fdRmin_ClustTrack)
1211           //MatchesToTrack = 1;
1212
1213
1214         }//_____________________________nTracks__________________________
1215
1216         h2_cells_M02  [centBin]->Fill(aodCluster->GetNCells(),aodCluster->GetM02());
1217         h2_Ellipse    [centBin]->Fill(aodCluster->GetM20(),aodCluster->GetM02());
1218         h1_Chi2       [centBin]->Fill(aodCluster->Chi2());//always -1.
1219         h1_nTrkMatch  [centBin]->Fill(aodCluster->GetNTracksMatched());
1220         h1_ClusterDisp[centBin]->Fill(aodCluster->GetDispersion());
1221         h2_E_time              ->Fill(aodCluster->E(),aodCluster->GetTOF());
1222
1223         // #################################################
1224         // track matching eop histograms are handled here...
1225         // #################################################
1226
1227         //uncomment this to do the track matching (3 of 3 lines, aod part)!!
1228         //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){
1229         if(isGoodAodCluster(aodCluster)){
1230           recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1231           E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ
1232           fAodEv->GetVertex(0)->GetXYZ(vertex);
1233           aodCluster->GetMomentum(Photon1,vertex);
1234           Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1235           Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1236           Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1237           Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1238           h1_E[centBin]->Fill(E1);
1239           h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1240           h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1241         }
1242         clusterPosition.Delete();
1243       }//if(aodCluster->IsEMCAL())
1244     }//if(fAodEv)
1245
1246   }//loop over nclusters.
1247
1248   //Make same event pions...
1249   for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1250     for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){
1251       Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j];
1252       Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1253       Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1254       Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/
1255                               (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt());
1256       Int_t asymCut = 0;
1257       if     (pairasym<0.1)  asymCut = 1;
1258       else if(pairasym<0.7)  asymCut = 2;
1259       else                   asymCut = 3;
1260
1261       h1_M        [centBin]->Fill(Parent.M());
1262       h3_MptAsymm [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1263       h2_dphi_deta[centBin]->Fill(deltaphi,deltaeta);
1264     }
1265   }
1266
1267   //Make mixed event...
1268   for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1269     for(unsigned int ipool=1; ipool<poolDepth; ipool++){
1270       for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){
1271         iskip = randy.Integer(Photons[0][izvtx][imult].size());
1272         if(j==iskip) continue;
1273         Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j];
1274         Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1275         Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1276         Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/
1277                                 (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt());
1278         Int_t asymCut = 0;
1279         if     (pairasym<0.1)  asymCut = 1;
1280         else if(pairasym<0.7)  asymCut = 2;
1281         else                   asymCut = 3;
1282
1283         h1_M_mix        [centBin]->Fill(Parent.M());
1284         h3_MptAsymm_mix [centBin]->Fill(Parent.M(),Parent.Pt(),asymCut);
1285         h2_dphi_deta_mix[centBin]->Fill(deltaphi,deltaeta);
1286       }
1287     }
1288   }
1289
1290   for(int ipool=poolDepth-1; ipool>0; ipool--){
1291     Photons[ipool][izvtx][imult].clear();
1292     for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++)
1293       Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]);
1294   }
1295   Photons[0][izvtx][imult].clear();
1296
1297
1298
1299   // NEW HISTO should be filled before this point, as PostData puts the
1300   // information for this iteration of the UserExec in the container
1301   PostData(1, fOutput);
1302   }
1303
1304 //________________________________________________________________________
1305 void AliAnalysisTaskEMCALMesonGGSDMpPb::Terminate(Option_t *) //specify what you want to have done
1306 {
1307   // Called once at the end of the query.
1308
1309 }
1310
1311 //________________________________________________________________________
1312 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetZvtxBin(Double_t vertZ)
1313 {
1314
1315   int izvtx = -1;
1316
1317   if     (vertZ<-35)
1318     izvtx=0;
1319   else if(vertZ<-30)
1320     izvtx=1;
1321   else if(vertZ<-25)
1322     izvtx=2;
1323   else if(vertZ<-20)
1324     izvtx=3;
1325   else if(vertZ<-15)
1326     izvtx=4;
1327   else if(vertZ<-10)
1328     izvtx=5;
1329   else if(vertZ< -5)
1330     izvtx=6;
1331   else if(vertZ<  0)
1332     izvtx=7;
1333   else if(vertZ<  5)
1334     izvtx=8;
1335   else if(vertZ< 10)
1336     izvtx=9;
1337   else if(vertZ< 15)
1338     izvtx=10;
1339   else if(vertZ< 20)
1340     izvtx=11;
1341   else if(vertZ< 25)
1342     izvtx=12;
1343   else if(vertZ< 30)
1344     izvtx=13;
1345   else if(vertZ< 35)
1346     izvtx=14;
1347   else
1348     izvtx=15;
1349
1350   return izvtx;
1351 }
1352
1353 //________________________________________________________________________
1354 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMultBin(Int_t mult){
1355
1356   int imult = -1;
1357
1358   if     (mult<2)
1359     imult=0;
1360   else if(mult<25)
1361     imult=mult-2;
1362   else
1363     imult=24;
1364
1365   return imult;
1366 }
1367
1368 //________________________________________________________________________
1369 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodEsdCluster(AliESDCaloCluster* esdclust){
1370
1371   int pass = 1;
1372   int nMinCells  = 2;
1373   double MinE    = 0.4;
1374   //double MinErat = 0;
1375   //double MinEcc  = 0;
1376
1377   if (!esdclust)
1378     pass = 0;
1379   if (!esdclust->IsEMCAL())
1380     pass = 0;
1381   if (esdclust->E()<MinE)
1382     pass = 0;
1383   if (esdclust->GetNCells()<nMinCells)
1384     pass = 0;
1385   //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
1386   //pass = 0;
1387   //if (esdclust->Chi2()<MinEcc) // eccentricity cut
1388   //pass = 0;//this is always -1.
1389
1390   //if(esdclust->GetM02()<0.1)
1391   //  pass = 0;
1392   //if(esdclust->GetM02()>0.5)
1393   //  pass = 0;
1394
1395   Float_t pos[3] = {0,0,0};
1396   esdclust->GetPosition(pos);
1397   TVector3 clusterPosition(pos);
1398   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1399      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
1400     pass = 0;
1401   clusterPosition.Delete();
1402
1403   //DOING THIS BY HAND NOW...
1404   //if(!esdclust->GetNTracksMatched()==0)
1405   //pass = 0;
1406
1407   return pass;
1408 }
1409
1410 //________________________________________________________________________
1411 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::isGoodAodCluster(AliAODCaloCluster* aodclust){
1412
1413   int pass = 1;
1414   int nMinCells  = 2;
1415   double MinE    = 0.4;
1416   //double MinErat = 0;
1417   //double MinEcc  = 0;
1418
1419   if (!aodclust)
1420     pass = 0;
1421   if (!aodclust->IsEMCAL())
1422     pass = 0;
1423   if (aodclust->E()<MinE)
1424     pass = 0;
1425   if (aodclust->GetNCells()<nMinCells)
1426     pass = 0;
1427   //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1428   //pass = 0;
1429   //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1430   //pass = 0;//this is always -1.
1431
1432   //if(aodclust->GetM02()<0.1)
1433   //pass = 0;
1434   //if(aodclust->GetM02()>0.5)
1435   //pass = 0;
1436
1437   Float_t pos[3] = {0,0,0};
1438   aodclust->GetPosition(pos);
1439   TVector3 clusterPosition(pos);
1440   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1441      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
1442     pass = 0;
1443   clusterPosition.Delete();
1444
1445   //DOING THIS BY HAND NOW...
1446   //if(!aodclust->GetNTracksMatched()==0)
1447   //pass = 0;
1448
1449   return pass;
1450 }
1451
1452 //________________________________________________________________________
1453 Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1454
1455   double dphi = p1.Phi() - p2.Phi();
1456
1457   if(dphi<0.5*TMath::Pi())
1458     dphi = dphi + 2.0*TMath::Pi();
1459
1460   if(dphi>1.5*TMath::Pi())
1461     dphi = dphi - 2.0*TMath::Pi();
1462
1463   return dphi;
1464 }
1465
1466 //________________________________________________________________________
1467 Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1468
1469   double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1470
1471   return deta;
1472 }
1473
1474
1475 //________________________________________________________________________
1476 Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1477
1478   double recalibfactor = 0.0;
1479
1480   if(iCalib==0){// no recalibration!
1481     recalibfactor = 1.0;
1482   }
1483   else if(iCalib==1){// just a scale factor:
1484     recalibfactor = 0.984;
1485   }
1486   else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1487     Double_t p[3] = {0.96968, -2.68720, -0.831607};
1488     recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1489   }
1490   else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1491     Double_t p[7] = {1.00000e+00, 3.04925e-02, 4.69043e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.00046e+00};
1492     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1493   }
1494   else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1495     Double_t p[7] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01};
1496     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1497   }
1498   else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1499     Double_t p[10] = {9.78672e-01, 2.39745e-01, 6.41199e-01, 9.13538e-02, 1.46058e+02, 1.99469e+01, 9.72716e-01, 0.96968, -2.68720, -0.831607};
1500     recalibfactor = ( (p[6]/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) ) / ( p[7] + exp(p[8] + p[9]*energy/2.0) );
1501   }
1502   else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1503     Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1504     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1505   }
1506   else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1507     Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1508     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1509   }
1510   else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1511     Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.964, -3.132, -0.435};
1512     //Double_t p[10] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286, 0.96968, -2.68720, -0.831607};//same SDM piece as iCalib==2
1513     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5]))))) * (p[7] + exp(p[8]+p[9]*energy*2.0));
1514   }
1515   else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1516     Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1517     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1518   }
1519   else if(iCalib==10){// Jason played with test beam data
1520     Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1521     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1522   }
1523   else if(iCalib==11){// Jason played with test beam MC
1524     Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1525     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1526   }
1527
1528   return recalibfactor;
1529 }
1530
1531
1532 //________________________________________________________________________
1533 Double_t AliAnalysisTaskEMCALMesonGGSDMpPb::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1534 {
1535   // Get maximum energy of attached cell.
1536
1537   id = -1;
1538   AliVCaloCells *fVCells=NULL;
1539   if(fEsdEv)      fVCells = fEsdEv->GetEMCALCells();
1540   else if(fAodEv) fVCells = fAodEv->GetEMCALCells();
1541   if(!fVCells)
1542     return 0;
1543
1544   Double_t maxe = 0;
1545   Int_t ncells = cluster->GetNCells();
1546   for (Int_t i=0; i<ncells; i++) {
1547     Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1548     if (e>maxe) {
1549       maxe = e;
1550       id   = cluster->GetCellAbsId(i);
1551     }
1552   }
1553   return maxe;
1554 }
1555
1556
1557 //________________________________________________________________________
1558 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){
1559
1560   AliMCParticle *mcP  = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1561
1562   Int_t nPTracks= mcEvent->GetNumberOfPrimaries();
1563
1564   Int_t isPhysPrimary   = 1;
1565   Int_t ismHF           = 0;
1566   Int_t ismLongLivedOrK = 0;
1567
1568   if(mcP->GetMother()<0)//if it has no mother...
1569     return isPhysPrimary;
1570
1571   Int_t imTrack = mcP->GetMother();
1572   AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1573
1574   if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark...
1575     return isPhysPrimary;
1576
1577
1578   //############################################
1579   //get the PDG digits....
1580   int num = mcPm->PdgCode();
1581   int RevDigits[10] = {0};
1582   int nDigits = 0;
1583   while (num >= 1){
1584     RevDigits[nDigits++] = num%10;
1585     num = num / 10;
1586   }
1587   //##############################################
1588
1589
1590   if(RevDigits[3]>3)//Baryons
1591     ismHF = 1;
1592   else if(RevDigits[2]>3)//Mesons
1593     ismHF = 1;
1594
1595   ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1596
1597   if(!ismHF && ismLongLivedOrK)
1598     isPhysPrimary = 0;
1599   else{ // check grandmother, greatgrandmothers, etc...
1600     while(imTrack >= nPTracks){
1601
1602       if(mcPm->GetMother()<0)//if it has no mother...
1603         break;
1604
1605       if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark...
1606         return isPhysPrimary;
1607
1608       imTrack = mcPm->GetMother();
1609       mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1610
1611       //############################################
1612       //get the PDG digits....
1613       num = mcPm->PdgCode();
1614       for(int i=0; i<10; i++)  RevDigits[i] = 0;
1615       nDigits = 0;
1616       while (num >= 1){
1617         RevDigits[nDigits++] = num%10;
1618         num = num / 10;
1619       }
1620       //##############################################
1621       if(RevDigits[3]>3)//Baryons
1622         ismHF = 1;
1623       else if(RevDigits[2]>3)//Mesons
1624         ismHF = 1;
1625
1626       ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1627
1628       if(!ismHF && ismLongLivedOrK)
1629         isPhysPrimary = 0;
1630
1631     }//while( >=nPTracks)
1632   }
1633
1634   return isPhysPrimary;
1635 }
1636
1637
1638 //________________________________________________________________________
1639 Int_t AliAnalysisTaskEMCALMesonGGSDMpPb::IsLongLivedOrK(Int_t MyPDGcode){
1640
1641   Int_t MyFlag = 0;
1642
1643   if(
1644      (TMath::Abs(MyPDGcode) == 22  ) ||        // Photon
1645      (TMath::Abs(MyPDGcode) == 11  ) ||        // Electron
1646      (TMath::Abs(MyPDGcode) == 13  ) ||        // Muon(-)
1647      (TMath::Abs(MyPDGcode) == 211 ) ||        // Pion
1648      (TMath::Abs(MyPDGcode) == 321 ) ||        // Kaon
1649      (TMath::Abs(MyPDGcode) == 310 ) ||        // K0s
1650      (TMath::Abs(MyPDGcode) == 130 ) ||        // K0l
1651      (TMath::Abs(MyPDGcode) == 2212) ||        // Proton
1652      (TMath::Abs(MyPDGcode) == 2112) ||        // Neutron
1653      (TMath::Abs(MyPDGcode) == 3122) ||        // Lambda_0
1654      (TMath::Abs(MyPDGcode) == 3112) ||        // Sigma Minus
1655      (TMath::Abs(MyPDGcode) == 3222) ||        // Sigma Plus
1656      (TMath::Abs(MyPDGcode) == 3312) ||        // Xsi Minus
1657      (TMath::Abs(MyPDGcode) == 3322) ||        // Xsi
1658      (TMath::Abs(MyPDGcode) == 3334) ||        // Omega
1659      (TMath::Abs(MyPDGcode) == 12  ) ||        // Electron Neutrino
1660      (TMath::Abs(MyPDGcode) == 14  ) ||        // Muon Neutrino
1661      (TMath::Abs(MyPDGcode) == 16  )   )       // Tau Neutrino
1662     MyFlag = 1;
1663
1664   return MyFlag;
1665 }