PWGGA/EMCALTasks/AliAnalysisTaskSDMGammaMC.cxx

   1 #include "AliAnalysisTaskSDMGammaMC.h"
   2
   3 #include <vector>
   4 #include <Riostream.h>
   5 #include <TChain.h>
   6 #include <TTree.h>
   7 #include <TF1.h>
   8 #include <TH1F.h>
   9 #include <TH2F.h>
  10 #include <TH3F.h>
  11 #include <TH1D.h>
  12 #include <TH2D.h>
  13 #include <TH3D.h>
  14 #include <TCanvas.h>
  15 #include <TList.h>
  16 #include <TFile.h>
  17 #include <TLorentzVector.h>
  18 #include <TNtuple.h>
  19 #include <TRandom3.h>
  20
  21 #include "AliAnalysisTaskSE.h"
  22 #include "AliAnalysisManager.h"
  23 #include "AliStack.h"
  24 #include "AliESDtrackCuts.h"
  25 #include "AliESDEvent.h"
  26 #include "AliESDInputHandler.h"
  27 #include "AliAODEvent.h"
  28 #include "AliMCEvent.h"
  29 #include "AliEMCALGeometry.h"
  30 #include "AliInputEventHandler.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliAODInputHandler.h"
  33
  34 #include "AliEMCALRecoUtils.h"
  35 #include "AliExternalTrackParam.h"
  36
  37 // ROOT includes
  38 #include <TGeoManager.h>
  39 #include <TGeoMatrix.h>
  40 #include <TGeoBBox.h>
  41 #include <TH2F.h>
  42 #include <TArrayI.h>
  43 #include <TArrayF.h>
  44 #include <TObjArray.h>
  45
  46 // STEER includes
  47 #include "AliVCluster.h"
  48 #include "AliVCaloCells.h"
  49 #include "AliLog.h"
  50 #include "AliPID.h"
  51 #include "AliESDEvent.h"
  52 #include "AliAODEvent.h"
  53 #include "AliESDtrack.h"
  54 #include "AliAODTrack.h"
  55 #include "AliExternalTrackParam.h"
  56 #include "AliESDfriendTrack.h"
  57 #include "AliTrackerBase.h"
  58
  59 // EMCAL includes
  60 #include "AliEMCALRecoUtils.h"
  61 #include "AliEMCALGeometry.h"
  62 #include "AliTrackerBase.h"
  63 #include "AliEMCALCalibTimeDepCorrection.h" // Run dependent
  64 #include "AliEMCALPIDUtils.h"
  65
  66 using std::cout;
  67 using std::endl;
  68
  69 ClassImp(AliAnalysisTaskSDMGammaMC)
  70
  71 //________________________________________________________________________
  72 AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC() :
  73   AliAnalysisTaskSE(),
  74   fOutput(0),
  75   fMcMode(0),
  76   fRecalibrator(0),
  77   fPhimin(0),
  78   fPhimax(0),
  79   fEtamin(0),
  80   fEtamax(0),
  81   fTrackCuts(0),
  82   fEsdEv(0),
  83   fAodEv(0),
  84   h1_nClusters(0),
  85   h1_zvtx(0),
  86   h1_trigger(0),
  87   h1_E(0),
  88   h1_Phi(0),
  89   h2_PiMotherID(0),
  90   h2_GaMotherID(0),
  91   h3_gE_RecTruth(0),
  92   h3_gE_RecTruth_ncellscut(0),
  93   h1_Pi0TruthPt(0),
  94   h1_PriPi0TruthPt(0),
  95   h1_Pi0TruthPtEmcal(0),
  96   h1_PriPi0TruthPtEmcal(0),
  97   h2_Pi0TruthPhiEta(0),
  98   h2_PriPi0TruthPhiEta(0),
  99   h2_Pi0TruthPhiEtaEmcal(0),
 100   h2_PriPi0TruthPhiEtaEmcal(0),
 101   h1_TruthPhotonsEmcal(0),
 102   h2_TruthPhotonsPhiEta(0),
 103   h1_PhotonsEmcal(0),
 104   h1_PhotonsNCellsCut(0),
 105   h1_PhotonsTrackMatchCut(0),
 106   h1_PhotonsAllCut(0),
 107   h2_PhotonsPhiEtaIsEmcal(0),
 108   h1_dR_RealMC(0),
 109   h1_Eta(0),
 110   h1_Chi2(0),
 111   h1_nTrkMatch(0),
 112   h1_nCells(0),
 113   h1_ClusterDisp(0),
 114   h2_Ellipse(0),
 115   h2_EtaPt(0),
 116   h2_dphi_deta(0),
 117   h2_dphi_deta_mix(0),
 118   h2_DispRes(0),
 119   h2_cells_M02(0),
 120   TriggerList(0)
 121 {
 122   // Dummy constructor ALWAYS needed for I/O.
 123 }
 124
 125 //________________________________________________________________________
 126 AliAnalysisTaskSDMGammaMC::AliAnalysisTaskSDMGammaMC(const char *name) :
 127   AliAnalysisTaskSE(name),
 128   fOutput(0),
 129   fMcMode(0),
 130   fRecalibrator(0),
 131   fPhimin(0),
 132   fPhimax(0),
 133   fEtamin(0),
 134   fEtamax(0),
 135   fTrackCuts(0),
 136   fEsdEv(0),
 137   fAodEv(0),
 138   h1_nClusters(0),
 139   h1_zvtx(0),
 140   h1_trigger(0),
 141   h1_E(0),
 142   h1_Phi(0),
 143   h2_PiMotherID(0),
 144   h2_GaMotherID(0),
 145   h3_gE_RecTruth(0),
 146   h3_gE_RecTruth_ncellscut(0),
 147   h1_Pi0TruthPt(0),
 148   h1_PriPi0TruthPt(0),
 149   h1_Pi0TruthPtEmcal(0),
 150   h1_PriPi0TruthPtEmcal(0),
 151   h2_Pi0TruthPhiEta(0),
 152   h2_PriPi0TruthPhiEta(0),
 153   h2_Pi0TruthPhiEtaEmcal(0),
 154   h2_PriPi0TruthPhiEtaEmcal(0),
 155   h1_TruthPhotonsEmcal(0),
 156   h2_TruthPhotonsPhiEta(0),
 157   h1_PhotonsEmcal(0),
 158   h1_PhotonsNCellsCut(0),
 159   h1_PhotonsTrackMatchCut(0),
 160   h1_PhotonsAllCut(0),
 161   h2_PhotonsPhiEtaIsEmcal(0),
 162   h1_dR_RealMC(0),
 163   h1_Eta(0),
 164   h1_Chi2(0),
 165   h1_nTrkMatch(0),
 166   h1_nCells(0),
 167   h1_ClusterDisp(0),
 168   h2_Ellipse(0),
 169   h2_EtaPt(0),
 170   h2_dphi_deta(0),
 171   h2_dphi_deta_mix(0),
 172   h2_DispRes(0),
 173   h2_cells_M02(0),
 174   TriggerList(0)
 175 {
 176   // Constructor
 177   // Define input and output slots here (never in the dummy constructor)
 178   // Input slot #0 works with a TChain - it is connected to the default input container
 179   // Output slot #1 writes into a TH1 container
 180
 181
 182   DefineOutput(1, TList::Class());                                            // for output list
 183 }
 184
 185 //________________________________________________________________________
 186 AliAnalysisTaskSDMGammaMC::~AliAnalysisTaskSDMGammaMC()
 187 {
 188   // Destructor. Clean-up the output list, but not the histograms that are put inside
 189   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 190   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 191     delete fOutput;
 192   }
 193   delete fTrackCuts;
 194 }
 195
 196 //________________________________________________________________________
 197 void AliAnalysisTaskSDMGammaMC::UserCreateOutputObjects()
 198 {
 199   // Create histograms
 200   // Called once (on the worker node)
 201
 202   fOutput = new TList();
 203   fOutput->SetOwner();  // IMPORTANT!
 204
 205   fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 206
 207   cout << "__________AliAnalysisTaskSDMGammaMC: Input settings__________" << endl;
 208   cout << " fMcMode:             " << fMcMode   << endl;
 209   cout << " fRecalibrator:       " << fRecalibrator << endl;
 210   cout << " phi range:           " << fPhimin << ", " << fPhimax << endl;
 211   cout << " eta range:           " << fEtamin << ", " << fEtamax << endl;
 212   cout << " number of zvtx bins: " << zvtx_bins << endl;
 213   cout << " number of mult bins: " << mult_bins << endl;
 214   cout << " poolDepth:           " << poolDepth << endl;
 215   cout << endl;
 216
 217
 218   double TotalNBins = 0.0;
 219
 220   // Create histograms
 221   Int_t nClustersbins = 501;
 222   Float_t nClusterslow = -0.5, nClustersup = 500.5;
 223   h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup);
 224   h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt");
 225   h1_nClusters->GetYaxis()->SetTitle("counts");
 226   h1_nClusters->SetMarkerStyle(kFullCircle);
 227   TotalNBins+=nClustersbins;
 228
 229   Int_t nZvertexbins = 501;
 230   Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
 231   h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
 232   h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
 233   h1_zvtx->GetYaxis()->SetTitle("counts");
 234   h1_zvtx->SetMarkerStyle(kFullCircle);
 235   TotalNBins+=nZvertexbins;
 236
 237   h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
 238   TotalNBins+=1001;
 239
 240   Int_t ptbins = 2000;
 241   Float_t ptlow = 0.0, ptup = 20.0;
 242   Int_t Ebins = 1000;
 243   Float_t Elow = 0.0, Eup = 20.0;
 244   h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup);
 245   h1_E->GetXaxis()->SetTitle("E [GeV]");
 246   h1_E->GetYaxis()->SetTitle("counts");
 247   h1_E->SetMarkerStyle(kFullCircle);
 248   TotalNBins+=Ebins;
 249
 250   h1_Phi = new TH1F("h1_Phi", "phi distribution", 1000, -7, 7);
 251   h1_Phi->GetXaxis()->SetTitle("#phi [rad]");
 252   h1_Phi->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 253   h1_Phi->SetMarkerStyle(kFullCircle);
 254   TotalNBins+=1000;
 255
 256   h2_PiMotherID = new TH2F("h2_PiMotherID", "Mother ID for Truth Pi0's", 100001, -0.5,100000.5, 2, 0.5,2.5);
 257   h2_PiMotherID->GetXaxis()->SetTitle("#pi^{0} Mother Particle ID");
 258   h2_PiMotherID->GetYaxis()->SetTitle("primary or non-primary");
 259   h2_PiMotherID->GetZaxis()->SetTitle("counts");
 260   h2_PiMotherID->SetMarkerStyle(kFullCircle);
 261   TotalNBins+=2*100001;
 262
 263   h2_GaMotherID = new TH2F("h2_GaMotherID", "Mother ID for Truth #gamma's", 100001, -0.5,100000.5, 2, 0.5,2.5);
 264   h2_GaMotherID->GetXaxis()->SetTitle("#gamma Mother Particle ID");
 265   h2_GaMotherID->GetYaxis()->SetTitle("primary or non-primary");
 266   h2_GaMotherID->GetZaxis()->SetTitle("counts");
 267   h2_GaMotherID->SetMarkerStyle(kFullCircle);
 268   TotalNBins+=2*100001;
 269
 270   h3_gE_RecTruth = new TH3F("h3_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
 271   h3_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
 272   h3_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
 273   h3_gE_RecTruth->GetZaxis()->SetTitle("category");
 274   h3_gE_RecTruth->SetMarkerStyle(kFullCircle);
 275   TotalNBins+=Ebins*1000*4;
 276
 277   h3_gE_RecTruth_ncellscut = new TH3F("h3_gE_RecTruth_ncellscut", "#gamma E_{truth}/E_{clust} vs E_{clust} (for nCells>1)", Ebins,Elow,Eup, 1000,0,4, 4,0.5,4.5);
 278   h3_gE_RecTruth_ncellscut->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
 279   h3_gE_RecTruth_ncellscut->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
 280   h3_gE_RecTruth_ncellscut->GetZaxis()->SetTitle("category");
 281   h3_gE_RecTruth_ncellscut->SetMarkerStyle(kFullCircle);
 282   TotalNBins+=Ebins*1000*4;
 283
 284   h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup);
 285   h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 286   h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 287   h1_Pi0TruthPt->SetMarkerStyle(kFullCircle);
 288   TotalNBins+=ptbins;
 289
 290   h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup);
 291   h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 292   h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 293   h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle);
 294   TotalNBins+=ptbins;
 295
 296   h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup);
 297   h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 298   h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 299   h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 300   TotalNBins+=ptbins;
 301
 302   h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
 303   h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 304   h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 305   h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 306   TotalNBins+=ptbins;
 307
 308   h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 309   h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 310   h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 311   TotalNBins+=380*200;
 312
 313   h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 314   h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 315   h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 316   TotalNBins+=380*200;
 317
 318   h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
 319   h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 320   h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 321   TotalNBins+=380*150;
 322
 323   h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
 324   h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 325   h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 326   TotalNBins+=380*150;
 327
 328   h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
 329   h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 330   h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 331   h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle);
 332   TotalNBins+=ptbins;
 333
 334   h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
 335                                    "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
 336   h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 337   h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
 338   h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
 339   TotalNBins+=380*150;
 340
 341   h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
 342   h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 343   h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 344   h1_PhotonsEmcal->SetMarkerStyle(kFullCircle);
 345   TotalNBins+=ptbins;
 346
 347   h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup);
 348   h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 349   h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 350   h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle);
 351   TotalNBins+=ptbins;
 352
 353   h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup);
 354   h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 355   h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 356   h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle);
 357   TotalNBins+=ptbins;
 358
 359   h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup);
 360   h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 361   h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 362   h1_PhotonsAllCut->SetMarkerStyle(kFullCircle);
 363   TotalNBins+=ptbins;
 364
 365   h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
 366                                      "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 150,-1.5,1.5);
 367   h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
 368   h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
 369   TotalNBins+=380*150;
 370
 371   h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10);
 372   h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
 373   h1_dR_RealMC->GetYaxis()->SetTitle("N");
 374   h1_dR_RealMC->SetMarkerStyle(kFullCircle);
 375   TotalNBins+=2000;
 376
 377   Int_t etabins = 150;
 378   Float_t etalow = -1.5, etaup = 1.5;
 379   h1_Eta = new TH1F("h1_Eta","#eta distribution for reconstructed",etabins, etalow, etaup);
 380   h1_Eta->GetXaxis()->SetTitle("#eta");
 381   h1_Eta->GetYaxis()->SetTitle("counts");
 382   TotalNBins+=etabins;
 383
 384   Int_t chi2bins = 100;
 385   Float_t chi2low = -2, chi2up = 2;
 386   h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up);
 387   h1_Chi2->GetXaxis()->SetTitle("#chi^{2}");
 388   h1_Chi2->GetYaxis()->SetTitle("counts");
 389   TotalNBins+=chi2bins;
 390
 391   h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5);
 392   h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched");
 393   h1_nTrkMatch->GetYaxis()->SetTitle("counts");
 394   TotalNBins+=14;
 395
 396   h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3);
 397   h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
 398   h1_ClusterDisp->GetYaxis()->SetTitle("counts");
 399   TotalNBins+=1000;
 400
 401   h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1);
 402   h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()");
 403   h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()");
 404   h2_Ellipse->GetZaxis()->SetTitle("counts");
 405   TotalNBins+=500*500;
 406
 407   h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup);
 408   h2_EtaPt->GetXaxis()->SetTitle("E [GeV]");
 409   h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 410   TotalNBins+=etabins*ptbins;
 411
 412   h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0);
 413   h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi");
 414   h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta");
 415   TotalNBins+=349*400;
 416
 417   h2_dphi_deta_mix = new TH2F("h2_dphi_deta_mix","#Delta#phi vs #Delta#eta (mixed events)", 349,-1.5,5, 400,-2.0,2.0);
 418   h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi");
 419   h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta");
 420   TotalNBins+=349*400;
 421
 422   h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2);
 423   h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
 424   h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
 425   h2_DispRes->GetZaxis()->SetTitle("counts");
 426   TotalNBins+=500*500;
 427
 428   h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5);
 429   h2_cells_M02->GetXaxis()->SetTitle("nCells");
 430   h2_cells_M02->GetYaxis()->SetTitle("M02");
 431   h2_cells_M02->GetZaxis()->SetTitle("counts");
 432   TotalNBins+=204*500;
 433
 434   cout << endl << "Total number of bins in booked histograms:  " << TotalNBins << endl << endl;
 435
 436   //TFile *f = OpenFile(1);
 437   //TDirectory::TContext context(f);
 438
 439   fOutput->Add(h1_nClusters);
 440   fOutput->Add(h1_zvtx);
 441   fOutput->Add(h1_trigger);
 442   fOutput->Add(h1_E);
 443   fOutput->Add(h1_Phi);
 444   fOutput->Add(h2_PiMotherID);
 445   fOutput->Add(h2_GaMotherID);
 446   fOutput->Add(h3_gE_RecTruth);
 447   fOutput->Add(h3_gE_RecTruth_ncellscut);
 448   fOutput->Add(h1_Pi0TruthPt);
 449   fOutput->Add(h1_PriPi0TruthPt);
 450   fOutput->Add(h1_Pi0TruthPtEmcal);
 451   fOutput->Add(h1_PriPi0TruthPtEmcal);
 452   fOutput->Add(h2_Pi0TruthPhiEta);
 453   fOutput->Add(h2_PriPi0TruthPhiEta);
 454   fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
 455   fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
 456   fOutput->Add(h1_TruthPhotonsEmcal);
 457   fOutput->Add(h2_TruthPhotonsPhiEta);
 458   fOutput->Add(h1_PhotonsEmcal);
 459   fOutput->Add(h1_PhotonsNCellsCut);
 460   fOutput->Add(h1_PhotonsTrackMatchCut);
 461   fOutput->Add(h1_PhotonsAllCut);
 462   fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
 463   fOutput->Add(h1_dR_RealMC);
 464   fOutput->Add(h1_Eta);
 465   fOutput->Add(h1_Chi2);
 466   fOutput->Add(h1_nTrkMatch);
 467   fOutput->Add(h1_ClusterDisp);
 468   fOutput->Add(h2_Ellipse);
 469   fOutput->Add(h2_EtaPt);
 470   fOutput->Add(h2_dphi_deta);
 471   fOutput->Add(h2_dphi_deta_mix);
 472   fOutput->Add(h2_DispRes);
 473   fOutput->Add(h2_cells_M02);
 474
 475   // Post data for ALL output slots >0 here,
 476   // To get at least an empty histogram
 477   // 1 is the outputnumber of a certain weg of task 1
 478   PostData(1, fOutput);
 479 }
 480
 481 //________________________________________________________________________
 482 void AliAnalysisTaskSDMGammaMC::UserExec(Option_t *)
 483 {
 484   // Main loop Called for each event
 485
 486   AliMCEvent *mcEvent = MCEvent();
 487   Bool_t isMC = bool(mcEvent);//is this the right way to do this?
 488   if (!mcEvent){
 489     cout << "no MC event" << endl;
 490     return;
 491   }
 492
 493   TRandom3 randy; randy.SetSeed(0);
 494
 495   double recalScale = 1.0;
 496
 497   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 498
 499   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
 500   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
 501   if (!aodH && !esdH)  Printf("ERROR: Could not get ESD or AODInputHandler");
 502
 503   if(esdH)      fEsdEv = esdH->GetEvent();
 504   else if(aodH) fAodEv = aodH->GetEvent();
 505   else{
 506     AliFatal("Neither ESD nor AOD event found");
 507     return;
 508   }
 509
 510
 511   // get pointer to reconstructed event
 512   AliVEvent *event = InputEvent();
 513   if (!event){
 514     AliError("Pointer == 0, this can not happen!");  return;}
 515   //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
 516   //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 517   //if (!fEsdEv){
 518   //AliError("Cannot get the ESD event");  return;}
 519
 520   Int_t iTrigger = 0;
 521   if (fEsdEv)       iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
 522   else if (fAodEv)  iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
 523   //h1_trigger->Fill(iTrigger);
 524
 525   char saythis[500];
 526   Int_t iTriggerBin = 0;
 527   for(unsigned long j=0; j<TriggerList.size(); j++){
 528     if(iTrigger==TriggerList[j])
 529       iTriggerBin=j+1;
 530   }
 531   if(iTriggerBin==0){
 532     TriggerList.push_back(iTrigger);
 533     iTriggerBin=TriggerList.size();
 534   }
 535
 536   h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
 537   sprintf(saythis,"%d",iTrigger);
 538   h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
 539
 540   if(fEsdEv){
 541     TString trigClasses = fEsdEv->GetFiredTriggerClasses();
 542     // remove "fast cluster events":
 543     if (trigClasses.Contains("FAST")  && !trigClasses.Contains("ALL"))
 544       return;
 545   }
 546   else if(fAodEv){
 547     TString trigClasses = fAodEv->GetFiredTriggerClasses();
 548     // remove "fast cluster events":
 549     if (trigClasses.Contains("FAST")  && !trigClasses.Contains("ALL"))
 550       return;
 551   }
 552
 553   if (fEsdEv){
 554     if(!(fEsdEv->GetPrimaryVertex()->GetStatus()))   return;
 555   }
 556   //else if (fAodEv){
 557   //if(!(fAodEv->GetPrimaryVertex()->GetStatus()))   return;
 558   //}
 559
 560   Double_t vertDisp=0.0;
 561   Double_t vertZres=0.0;
 562   Bool_t vertIsfromZ=0;
 563   if (fEsdEv){
 564     vertDisp    = fEsdEv->GetPrimaryVertex()->GetDispersion();
 565     vertZres    = fEsdEv->GetPrimaryVertex()->GetZRes();
 566     vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
 567   }
 568   else if (fAodEv){
 569     vertDisp    = 0;
 570     vertZres    = 0;
 571     vertIsfromZ = 0;
 572   }
 573
 574   h2_DispRes->Fill(vertDisp, vertZres);
 575   // if vertex is from spd vertexZ, require more stringent cut
 576   if (vertIsfromZ) {
 577     if (vertDisp>0.02 ||  vertZres>0.25 )
 578       return; // bad vertex from VertexerZ
 579   }
 580
 581   Int_t nclusters=0;
 582   if(fEsdEv){
 583     //Int_t evtN      = fEsdEv->GetEventNumberInFile();
 584     //Int_t ntracks   = fEsdEv->GetNumberOfTracks();
 585     nclusters = fEsdEv->GetNumberOfCaloClusters();
 586   }
 587   else if(fAodEv){
 588     //Int_t evtN      = fAodEv->GetEventNumberInFile();
 589     //Int_t ntracks   = fAodEv->GetNumberOfTracks();
 590     nclusters = fAodEv->GetNumberOfCaloClusters();
 591   }
 592
 593   // EMCal cluster loop for reconstructed event
 594   //numberofclusters set above!
 595   Double_t vertZ=0.0;
 596   if (fEsdEv)       vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
 597   else if (fAodEv)  vertZ = fAodEv->GetPrimaryVertex()->GetZ();
 598
 599   h1_zvtx->Fill(vertZ);
 600   //zvertex cut:
 601   if(fabs(vertZ)>10.0)
 602     return;
 603
 604   h1_nClusters->Fill(nclusters);
 605
 606   //cout << iskip << " " << izvtx << " " << imult << endl;
 607   //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
 608
 609
 610
 611   //######################### ~~~~~~~~~~~ ##################################
 612   //######################### STARTING MC ##################################
 613   //######################### ~~~~~~~~~~~ ##################################
 614
 615   if(isMC){
 616     int isPrimary  = 0;
 617     int isK0sDecay = 0;
 618
 619     if (!mcEvent){
 620       cout << "no MC event" << endl;
 621       return;
 622     }
 623
 624     const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
 625     if (!evtVtx)
 626       return;
 627
 628     mcEvent->PreReadAll();
 629
 630     Int_t nTracksMC  = mcEvent->GetNumberOfTracks();
 631     Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
 632
 633     for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
 634       AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
 635       if (!mcP)
 636         continue;
 637
 638       if(iTrack<nPTracksMC)  isPrimary = 1;
 639       else                   isPrimary = 0;
 640
 641       if(mcP->PdgCode() == 22){
 642         if(isPrimary==1){
 643           if(mcP->GetMother()>-1)
 644             h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
 645           else
 646             h2_GaMotherID->Fill(0.0,1);
 647         }
 648         else
 649           h2_GaMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
 650       }
 651
 652       // it's a pion !!
 653       if(mcP->PdgCode() != 111)
 654         continue;
 655
 656       isK0sDecay = 0;
 657       if(mcP->GetMother()>-1){
 658         if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  310 ||
 659             ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310  )
 660           isK0sDecay = 1;
 661       }
 662
 663       // primary particle
 664       //Double_t dR_vtx = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
 665       //                            (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
 666       //if(dR_vtx <= 0.01)  isPrimary = 1;
 667       //else            isPrimary = 0;
 668
 669
 670       if(isPrimary==1){
 671         if(mcP->GetMother()>-1)
 672           h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 1);
 673         else
 674           h2_PiMotherID->Fill(0.0,1);
 675       }
 676       else
 677         h2_PiMotherID->Fill(( (AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()) )->PdgCode(), 2);
 678
 679       h1_Pi0TruthPt    ->Fill(mcP->Pt());
 680       h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
 681
 682       if(isPrimary==1){
 683         h1_PriPi0TruthPt    ->Fill(mcP->Pt());
 684         h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
 685       }
 686
 687       if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
 688         continue;
 689
 690       Int_t DecayPhotonLabel[2] = {mcP->GetFirstDaughter(),
 691                                    mcP->GetLastDaughter() };
 692
 693       if (DecayPhotonLabel[0]<0)  continue;
 694       if (DecayPhotonLabel[1]<0)  DecayPhotonLabel[1]=DecayPhotonLabel[0];
 695       if (DecayPhotonLabel[1]-DecayPhotonLabel[0] != 1)  continue;
 696
 697       bool bacc = true;
 698       bool binp = true;
 699       bool isConv[2] = {1,1};
 700       Int_t convIndices[2][2] = { {-1,-1},{-1,-1} };
 701       Double_t eta_d[2] = {0.0,0.0};
 702       Double_t phi_d[2] = {0.0,0.0};
 703       Int_t daughter_index = -1;
 704       for (Int_t iPhoton=DecayPhotonLabel[0];iPhoton<=DecayPhotonLabel[1];++iPhoton){
 705         if(iPhoton==DecayPhotonLabel[0]) daughter_index=0;
 706         else                             daughter_index=1;
 707         const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(iPhoton));
 708         eta_d[daughter_index] = dmc->Eta();
 709         phi_d[daughter_index] = dmc->Phi();
 710         if(!(dmc->PdgCode()==22))             binp = false;
 711         if(!(eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax &&
 712              phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax   ))   bacc = false;
 713
 714         if( ((TParticle*)dmc->Particle())->GetNDaughters() != 2 )  isConv[daughter_index] = 0;
 715         else{//if photon has 2 daughters.
 716
 717           Int_t dd1 = dmc->GetFirstDaughter();
 718           Int_t dd2 = dmc->GetLastDaughter();
 719           if (dd2-dd1 != 1)  cout << "How can this happen???? " << endl;
 720           const AliMCParticle *dd1mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd1));
 721           const AliMCParticle *dd2mc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dd2));
 722           if( dd1mc->PdgCode() != -dd2mc->PdgCode() )
 723             isConv[daughter_index] = 0;
 724           else if( TMath::Abs(dd1mc->PdgCode())!=11 )
 725             isConv[daughter_index] = 0;
 726           if(isConv[daughter_index]==1){//store the e+e- indices...
 727             convIndices[daughter_index][0] = dd1;
 728             convIndices[daughter_index][1] = dd2;
 729           }//if this photon converted.
 730         }//close else-if photon has 2 daughters.
 731       }//loop over 2 decay photons (iPhoton)
 732
 733       if(binp && bacc){// 2 Photons hit the EMCAL!
 734
 735         h1_Pi0TruthPtEmcal    ->Fill(mcP->Pt());
 736         h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
 737
 738         if(isPrimary==1){
 739           h1_PriPi0TruthPtEmcal    ->Fill(mcP->Pt());
 740           h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
 741         }
 742
 743         Int_t PhotonClusterMatch[2][3]  = { {0,-1,-1},
 744                                             {0,-1,-1} };
 745         Int_t PhotonElectronMatch[2][6] = { {0,-1,-1,-1,-1,-1},
 746                                             {0,-1,-1,-1,-1,-1} };
 747
 748         for(int iCluster=0; iCluster<nclusters; iCluster++) {
 749
 750           AliESDCaloCluster* esdCluster=NULL;
 751           AliAODCaloCluster* aodCluster=NULL;
 752           if (fEsdEv)       esdCluster = fEsdEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
 753           else if (fAodEv)  aodCluster = fAodEv->GetCaloCluster(iCluster); // pointer to EMCal cluster
 754
 755           Double_t clustMC_phi, clustMC_eta;
 756           if(fEsdEv){
 757             if(esdCluster->IsEMCAL()){
 758
 759               if(!isGoodEsdCluster(esdCluster))
 760                 continue;
 761
 762               Float_t pos[3] = {0,0,0};
 763               esdCluster->GetPosition(pos);
 764               TVector3 vpos(pos);
 765               h1_Phi->Fill(vpos.Phi());
 766               clustMC_phi = vpos.Phi();
 767               clustMC_eta = vpos.Eta();
 768
 769               Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
 770                                         (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
 771               h1_dR_RealMC->Fill(dR);
 772               //matches_pion_photon = 0;
 773               //if(dR<=0.04) matches_pion_photon = 1;
 774
 775               TArrayI *TruthLabelsA = esdCluster->GetLabelsArray();
 776               if(TruthLabelsA){
 777                 for(int itl=0; itl<TruthLabelsA->GetSize(); itl++){
 778
 779                   for(int iPhoton=0; iPhoton<2; iPhoton++){
 780                     if(TruthLabelsA->At(itl)==DecayPhotonLabel[iPhoton]){
 781                       PhotonClusterMatch[iPhoton][0] = 1;
 782                       PhotonClusterMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
 783                       PhotonClusterMatch[iPhoton][2] = iCluster;
 784                     }
 785                   }//loop over truth labels.
 786
 787                   AliMCParticle *elecCandidate = (AliMCParticle*)(mcEvent->GetTrack(TruthLabelsA->At(itl)));
 788                   if(TMath::Abs(elecCandidate->PdgCode())==11){//if we have an electron...
 789                     Int_t elecMother_index = elecCandidate->GetMother();
 790                     if(elecMother_index>1 && elecMother_index<nTracksMC){
 791                       AliMCParticle *elecMother   = (AliMCParticle*)(mcEvent->GetTrack(elecMother_index));
 792                       if( TMath::Abs(elecMother->PdgCode())==22 ){//if the e's mother is a photon...
 793                         Int_t elecGrandMother_index = elecMother->GetMother();
 794                         if(elecGrandMother_index==iTrack){//if the e's gMother is THE pi0 in question...
 795                           AliMCParticle *elecGrandMother = (AliMCParticle*)(mcEvent->GetTrack(elecGrandMother_index));
 796                           if( TMath::Abs(elecGrandMother->PdgCode())!=111 ) cout << "|| This can't happen!!  A pion is a pion is a pion is a pion... ||" << endl;
 797
 798                           for(int iPhoton=0; iPhoton<2; iPhoton++){
 799                             //if(convIndices[iPhoton][0]==elecMother_index){
 800                             if(convIndices[iPhoton][0]==TruthLabelsA->At(itl) || convIndices[iPhoton][1]==TruthLabelsA->At(itl)){
 801                               if(PhotonElectronMatch[iPhoton][1] == DecayPhotonLabel[iPhoton]) PhotonElectronMatch[iPhoton][0] = 2;
 802                               else                                                             PhotonElectronMatch[iPhoton][0] = 1;
 803                               PhotonElectronMatch[iPhoton][1] = DecayPhotonLabel[iPhoton];
 804                               if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][2] = iCluster;//first cluster
 805                               else                                    PhotonElectronMatch[iPhoton][3] = iCluster;//second cluster
 806                               if(PhotonElectronMatch[iPhoton][2]==-1) PhotonElectronMatch[iPhoton][4] = elecCandidate->PdgCode();
 807                               else                                    PhotonElectronMatch[iPhoton][5] = elecCandidate->PdgCode();
 808                             }
 809                           }//loop over both decay photons
 810
 811                         }//if it's THE pi0.
 812                       }//if we have a photon.
 813                     }//if we have an electron.
 814                   }//if the candidate has a real mother.
 815
 816                 }//itl (TruthLabel loop)
 817               }//if(TruthLabelsA exists)
 818             }//if(isEMCal)
 819           }//if(esdEv)
 820         }//loop over clusters.
 821
 822
 823         for(int iPhoton=0; iPhoton<2; iPhoton++){
 824           AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(DecayPhotonLabel[iPhoton]));
 825           if(!truthP)
 826             continue;
 827           if(PhotonClusterMatch[iPhoton][0]==1){
 828             AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonClusterMatch[iPhoton][2]);
 829             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
 830             h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
 831             if(esdCluster->GetNCells()>=2)
 832               h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 1);
 833           }
 834           else if(PhotonElectronMatch[iPhoton][0]==2 && PhotonElectronMatch[iPhoton][2] == PhotonElectronMatch[iPhoton][3]){//merged conv photon
 835             AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
 836             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
 837             h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
 838             if(esdCluster->GetNCells()>=2)
 839               h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 2);
 840           }
 841           else if(PhotonElectronMatch[iPhoton][0]==2){//non-merged conv photon (but both hit emcal)
 842             AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
 843             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
 844             h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
 845             if(esdCluster->GetNCells()>=2)
 846               h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
 847                                esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][3]);
 848             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
 849             h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
 850             if(esdCluster->GetNCells()>=2)
 851               h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 4);
 852           }
 853           else if(PhotonElectronMatch[iPhoton][0]==1){//non-merged conv photon (one missed emcal)
 854             AliESDCaloCluster *esdCluster = fEsdEv->GetCaloCluster(PhotonElectronMatch[iPhoton][2]);
 855             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
 856             h3_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
 857             if(esdCluster->GetNCells()>=2)
 858               h3_gE_RecTruth_ncellscut->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()), 3);
 859           }
 860         }//loop over decay photons (iPhoton).
 861
 862       }// 2 Photons pointed at the EMCAL!
 863     }//for(nTracksMC) ie. Truth Pion loop.
 864
 865   }//if(isMC)
 866
 867   //######################### ~~~~~~~~~~~~ ##################################
 868   //######################### DONE WITH MC ##################################
 869   //######################### ~~~~~~~~~~~~ ##################################
 870
 871
 872
 873   // NEW HISTO should be filled before this point, as PostData puts the
 874   // information for this iteration of the UserExec in the container
 875   PostData(1, fOutput);
 876   }
 877
 878 //________________________________________________________________________
 879 void AliAnalysisTaskSDMGammaMC::Terminate(Option_t *) //specify what you want to have done
 880 {
 881   // Called once at the end of the query.
 882
 883 }
 884
 885 //________________________________________________________________________
 886 Int_t AliAnalysisTaskSDMGammaMC::GetZvtxBin(Double_t vertZ)
 887 {
 888
 889   int izvtx = -1;
 890
 891   if     (vertZ<-35)
 892     izvtx=0;
 893   else if(vertZ<-30)
 894     izvtx=1;
 895   else if(vertZ<-25)
 896     izvtx=2;
 897   else if(vertZ<-20)
 898     izvtx=3;
 899   else if(vertZ<-15)
 900     izvtx=4;
 901   else if(vertZ<-10)
 902     izvtx=5;
 903   else if(vertZ< -5)
 904     izvtx=6;
 905   else if(vertZ<  0)
 906     izvtx=7;
 907   else if(vertZ<  5)
 908     izvtx=8;
 909   else if(vertZ< 10)
 910     izvtx=9;
 911   else if(vertZ< 15)
 912     izvtx=10;
 913   else if(vertZ< 20)
 914     izvtx=11;
 915   else if(vertZ< 25)
 916     izvtx=12;
 917   else if(vertZ< 30)
 918     izvtx=13;
 919   else if(vertZ< 35)
 920     izvtx=14;
 921   else
 922     izvtx=15;
 923
 924   return izvtx;
 925 }
 926
 927 //________________________________________________________________________
 928 Int_t AliAnalysisTaskSDMGammaMC::GetMultBin(Int_t mult){
 929
 930   int imult = -1;
 931
 932   if     (mult<2)
 933     imult=0;
 934   else if(mult<25)
 935     imult=mult-2;
 936   else
 937     imult=24;
 938
 939   return imult;
 940 }
 941
 942 //________________________________________________________________________
 943 Int_t AliAnalysisTaskSDMGammaMC::isGoodEsdCluster(AliESDCaloCluster* esdclust){
 944
 945   int pass = 1;
 946   int nMinCells  = 1;
 947   double MinE    = 0.4;
 948   //double MinErat = 0;
 949   //double MinEcc  = 0;
 950
 951   if (!esdclust)
 952     pass = 0;
 953   if (!esdclust->IsEMCAL())
 954     pass = 0;//removes ~70% of clusters.
 955   if (esdclust->E()<MinE)
 956     pass = 0;//does nothing
 957   if (esdclust->GetNCells()<nMinCells)
 958     pass = 0;//does nothing
 959   //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
 960   //pass = 0;
 961   //if (esdclust->Chi2()<MinEcc) // eccentricity cut
 962   //pass = 0;//this is always -1.
 963
 964   /*
 965   //This cuts out more than just 1 cell clusters
 966   //and drains the statistics badly.
 967   //haven't figured out what it does yet.
 968   if(esdclust->GetM20()<0.02)
 969   pass = 0;
 970   */
 971   //if(esdclust->GetM02()<0.1)
 972   //  pass = 0;
 973   //if(esdclust->GetM02()>0.5)
 974   //  pass = 0;
 975   //if(esdclust->GetNCells()<2)
 976   //  pass = 0;
 977
 978   Float_t pos[3] = {0,0,0};
 979   esdclust->GetPosition(pos);
 980   TVector3 clusterPosition(pos);
 981   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
 982      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
 983     pass = 0;
 984   clusterPosition.Delete();
 985
 986   //doing this by hand now...
 987   //if(!esdclust->GetNTracksMatched()==0)
 988   //pass = 0;
 989
 990   return pass;
 991 }
 992
 993 //________________________________________________________________________
 994 Int_t AliAnalysisTaskSDMGammaMC::isGoodAodCluster(AliAODCaloCluster* aodclust){
 995
 996   int pass = 1;
 997   int nMinCells  = 1;
 998   double MinE    = 0.4;
 999   //double MinErat = 0;
1000   //double MinEcc  = 0;
1001
1002   if (!aodclust)
1003     pass = 0;
1004   if (!aodclust->IsEMCAL())
1005     pass = 0;//removes ~70% of clusters.
1006   if (aodclust->E()<MinE)
1007     pass = 0;//does nothing
1008   if (aodclust->GetNCells()<nMinCells)
1009     pass = 0;//does nothing
1010   //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1011   //pass = 0;
1012   //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1013   //pass = 0;//this is always -1.
1014
1015   /*
1016   //This cuts out more than just 1 cell clusters
1017   //and drains the statistics badly.
1018   //haven't figured out what it does yet.
1019   if(aodclust->GetM20()<0.02)
1020   pass = 0;
1021   if(aodclust->GetM02()<0.02)
1022   pass = 0;
1023   */
1024   //if(aodclust->GetM02()<0.1)
1025   //  pass = 0;
1026   //if(aodclust->GetM02()>0.5)
1027   //  pass = 0;
1028   //if(aodclust->GetNCells()<2)
1029   //  pass = 0;
1030
1031   Float_t pos[3] = {0,0,0};
1032   aodclust->GetPosition(pos);
1033   TVector3 clusterPosition(pos);
1034   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1035      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
1036     pass = 0;
1037   clusterPosition.Delete();
1038
1039   //if(!aodclust->GetNTracksMatched()==0)
1040   //pass = 0;
1041
1042   return pass;
1043 }
1044
1045 //________________________________________________________________________
1046 Double_t AliAnalysisTaskSDMGammaMC::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1047
1048   double dphi = p1.Phi() - p2.Phi();
1049
1050   if(dphi<0.5*TMath::Pi())
1051     dphi = dphi + 2.0*TMath::Pi();
1052
1053   if(dphi>1.5*TMath::Pi())
1054     dphi = dphi - 2.0*TMath::Pi();
1055
1056   return dphi;
1057 }
1058
1059 //________________________________________________________________________
1060 Double_t AliAnalysisTaskSDMGammaMC::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1061
1062   double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1063
1064   return deta;
1065 }
1066
1067
1068 //________________________________________________________________________
1069 Double_t AliAnalysisTaskSDMGammaMC::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1070
1071   double recalibfactor = 0.0;
1072
1073   if(iCalib==0){// no recalibration!
1074     recalibfactor = 1.0;
1075   }
1076   else if(iCalib==1){// just a scale factor:
1077     recalibfactor = 0.984;
1078   }
1079   else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1080     Double_t p[3] = {0.96968, -2.68720, -0.831607};
1081     recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1082   }
1083   else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1084     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};
1085     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1086   }
1087   else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1088     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};
1089     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1090   }
1091   else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1092     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};
1093     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) );
1094   }
1095   else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1096     Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1097     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1098   }
1099   else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1100     Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1101     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1102   }
1103   else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1104     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};
1105     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));
1106   }
1107   else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1108     Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1109     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1110   }
1111   else if(iCalib==10){// Jason played with test beam data
1112     Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1113     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1114   }
1115   else if(iCalib==11){// Jason played with test beam MC
1116     Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1117     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1118   }
1119
1120   return recalibfactor;
1121 }
1122