PWGGA/EMCALTasks/AliAnalysisTaskEMCALMesonGGSDM.cxx

   1 #include "AliAnalysisTaskEMCALMesonGGSDM.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 #include "AliGenCocktailEventHeader.h"
  67
  68 using std::cout;
  69 using std::endl;
  70
  71 ClassImp(AliAnalysisTaskEMCALMesonGGSDM)
  72
  73 //________________________________________________________________________
  74 AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM() :
  75   AliAnalysisTaskSE(),
  76   fOutput(0),
  77   fMcMode(0),
  78   fMyMCType(0),
  79   fRecalibrator(0),
  80   fdRmin_ClustTrack(0),
  81   fPhimin(0),
  82   fPhimax(0),
  83   fEtamin(0),
  84   fEtamax(0),
  85   fTrackCuts(0),
  86   fEsdEv(0),
  87   fAodEv(0),
  88   h1_nClusters(0),
  89   h1_zvtx(0),
  90   h1_trigger(0),
  91   h1_M(0),
  92   h1_M_mix(0),
  93   h1_E(0),
  94   h2_PhiEtaCluster(0),
  95   h2_PhiEtaClusterCut(0),
  96   h2_PhiEtaMaxCell(0),
  97   h2_PhiEtaMaxCellCut(0),
  98   h1_dR_ClustTrk(0),
  99   h2_gE_RecTruth(0),
 100   h2_eop_E(0),
 101   h2_eop_pT(0),
 102   h2_E_time(0),
 103   h1_Pi0TruthPt(0),
 104   h1_K0Pi0TruthPt(0),
 105   h1_PriPi0TruthPt(0),
 106   h1_PhysPi0TruthPt(0),
 107   h1_Pi0TruthPtEmcal(0),
 108   h1_K0Pi0TruthPtEmcal(0),
 109   h1_PriPi0TruthPtEmcal(0),
 110   h1_PhysPi0TruthPtEmcal(0),
 111   h1_Pi0TruthPtPhi2piEta065(0),
 112   h1_K0Pi0TruthPtPhi2piEta065(0),
 113   h1_PriPi0TruthPtPhi2piEta065(0),
 114   h1_PhysPi0TruthPtPhi2piEta065(0),
 115   h1_Pi0TruthPtPhi2piEta1(0),
 116   h1_K0Pi0TruthPtPhi2piEta1(0),
 117   h1_PriPi0TruthPtPhi2piEta1(0),
 118   h1_PhysPi0TruthPtPhi2piEta1(0),
 119   h2_Pi0TruthPhiEta(0),
 120   h2_PriPi0TruthPhiEta(0),
 121   h2_Pi0TruthPhiEtaEmcal(0),
 122   h2_PriPi0TruthPhiEtaEmcal(0),
 123   h1_TruthPhotonsEmcal(0),
 124   h2_TruthPhotonsPhiEta(0),
 125   h1_PhotonsEmcal(0),
 126   h1_PhotonsNCellsCut(0),
 127   h1_PhotonsTrackMatchCut(0),
 128   h1_PhotonsAllCut(0),
 129   h2_PhotonsPhiEtaIsEmcal(0),
 130   h1_dR_RealMC(0),
 131   h2_Mpt_Pri(0),
 132   h2_Mpt_Sec(0),
 133   h3_MptR_Sec(0),
 134   h3_MptR_K0s(0),
 135   h3_MptR_Mat(0),
 136   h2_PtR_MatM(0),
 137   h2_Mpt_Pri_conv(0),
 138   h2_Mpt_Sec_conv(0),
 139   h3_MptR_Sec_conv(0),
 140   h3_MptR_K0s_conv(0),
 141   h3_MptR_Mat_conv(0),
 142   h1_eConversionR(0),
 143   h1_PriPi0Mother(0),
 144   h1_SecPi0Mother(0),
 145   h1_Chi2(0),
 146   h1_nTrkMatch(0),
 147   h1_nCells(0),
 148   h1_ClusterDisp(0),
 149   h2_Ellipse(0),
 150   h2_EtaPt(0),
 151 //h2_Mpt(0),
 152   h3_MptAsymm(0),
 153 //h2_Mpt_mix(0),
 154   h3_MptAsymm_mix(0),
 155   h2_dphi_deta(0),
 156   h2_dphi_deta_mix(0),
 157   h2_DispRes(0),
 158   h2_cells_M02(0),
 159   TriggerList(0),
 160   fHelperClass(0)
 161 {
 162   // Dummy constructor ALWAYS needed for I/O.
 163 }
 164
 165 //________________________________________________________________________
 166 AliAnalysisTaskEMCALMesonGGSDM::AliAnalysisTaskEMCALMesonGGSDM(const char *name) :
 167   AliAnalysisTaskSE(name),
 168   fOutput(0),
 169   fMcMode(0),
 170   fMyMCType(0),
 171   fRecalibrator(0),
 172   fdRmin_ClustTrack(0),
 173   fPhimin(0),
 174   fPhimax(0),
 175   fEtamin(0),
 176   fEtamax(0),
 177   fTrackCuts(0),
 178   fEsdEv(0),
 179   fAodEv(0),
 180   h1_nClusters(0),
 181   h1_zvtx(0),
 182   h1_trigger(0),
 183   h1_M(0),
 184   h1_M_mix(0),
 185   h1_E(0),
 186   h2_PhiEtaCluster(0),
 187   h2_PhiEtaClusterCut(0),
 188   h2_PhiEtaMaxCell(0),
 189   h2_PhiEtaMaxCellCut(0),
 190   h1_dR_ClustTrk(0),
 191   h2_gE_RecTruth(0),
 192   h2_eop_E(0),
 193   h2_eop_pT(0),
 194   h2_E_time(0),
 195   h1_Pi0TruthPt(0),
 196   h1_K0Pi0TruthPt(0),
 197   h1_PriPi0TruthPt(0),
 198   h1_PhysPi0TruthPt(0),
 199   h1_Pi0TruthPtEmcal(0),
 200   h1_K0Pi0TruthPtEmcal(0),
 201   h1_PriPi0TruthPtEmcal(0),
 202   h1_PhysPi0TruthPtEmcal(0),
 203   h1_Pi0TruthPtPhi2piEta065(0),
 204   h1_K0Pi0TruthPtPhi2piEta065(0),
 205   h1_PriPi0TruthPtPhi2piEta065(0),
 206   h1_PhysPi0TruthPtPhi2piEta065(0),
 207   h1_Pi0TruthPtPhi2piEta1(0),
 208   h1_K0Pi0TruthPtPhi2piEta1(0),
 209   h1_PriPi0TruthPtPhi2piEta1(0),
 210   h1_PhysPi0TruthPtPhi2piEta1(0),
 211   h2_Pi0TruthPhiEta(0),
 212   h2_PriPi0TruthPhiEta(0),
 213   h2_Pi0TruthPhiEtaEmcal(0),
 214   h2_PriPi0TruthPhiEtaEmcal(0),
 215   h1_TruthPhotonsEmcal(0),
 216   h2_TruthPhotonsPhiEta(0),
 217   h1_PhotonsEmcal(0),
 218   h1_PhotonsNCellsCut(0),
 219   h1_PhotonsTrackMatchCut(0),
 220   h1_PhotonsAllCut(0),
 221   h2_PhotonsPhiEtaIsEmcal(0),
 222   h1_dR_RealMC(0),
 223   h2_Mpt_Pri(0),
 224   h2_Mpt_Sec(0),
 225   h3_MptR_Sec(0),
 226   h3_MptR_K0s(0),
 227   h3_MptR_Mat(0),
 228   h2_PtR_MatM(0),
 229   h2_Mpt_Pri_conv(0),
 230   h2_Mpt_Sec_conv(0),
 231   h3_MptR_Sec_conv(0),
 232   h3_MptR_K0s_conv(0),
 233   h3_MptR_Mat_conv(0),
 234   h1_eConversionR(0),
 235   h1_PriPi0Mother(0),
 236   h1_SecPi0Mother(0),
 237   h1_Chi2(0),
 238   h1_nTrkMatch(0),
 239   h1_nCells(0),
 240   h1_ClusterDisp(0),
 241   h2_Ellipse(0),
 242   h2_EtaPt(0),
 243 //h2_Mpt(0),
 244   h3_MptAsymm(0),
 245 //h2_Mpt_mix(0),
 246   h3_MptAsymm_mix(0),
 247   h2_dphi_deta(0),
 248   h2_dphi_deta_mix(0),
 249   h2_DispRes(0),
 250   h2_cells_M02(0),
 251   TriggerList(0),
 252   fHelperClass(0)
 253 {
 254   // Constructor
 255   // Define input and output slots here (never in the dummy constructor)
 256   // Input slot #0 works with a TChain - it is connected to the default input container
 257   // Output slot #1 writes into a TH1 container
 258
 259   DefineOutput(1, TList::Class());                                            // for output list
 260 }
 261
 262 //________________________________________________________________________
 263 AliAnalysisTaskEMCALMesonGGSDM::~AliAnalysisTaskEMCALMesonGGSDM()
 264 {
 265   // Destructor. Clean-up the output list, but not the histograms that are put inside
 266   // (the list is owner and will clean-up these histograms). Protect in PROOF case.
 267   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {
 268     delete fOutput;
 269   }
 270   delete fTrackCuts;
 271 }
 272
 273 //________________________________________________________________________
 274 void AliAnalysisTaskEMCALMesonGGSDM::UserCreateOutputObjects()
 275 {
 276   // Create histograms
 277   // Called once (on the worker node)
 278
 279   fOutput = new TList();
 280   fOutput->SetOwner();  // IMPORTANT!
 281
 282   fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 283
 284   cout << "__________AliAnalysisTaskEMCALMesonGGSDM: Input settings__________" << endl;
 285   cout << " fMcMode:             " << fMcMode       << endl;
 286   cout << " fMyMCType:           " << fMyMCType     << endl;
 287   cout << " fRecalibrator:       " << fRecalibrator << endl;
 288   cout << " dRmin_ClustTrack:    " << fdRmin_ClustTrack << endl;
 289   cout << " phi range:           " << fPhimin << ", " << fPhimax << endl;
 290   cout << " eta range:           " << fEtamin << ", " << fEtamax << endl;
 291   cout << " number of zvtx bins: " << zvtx_bins     << endl;
 292   cout << " number of mult bins: " << mult_bins     << endl;
 293   cout << " poolDepth:           " << poolDepth     << endl;
 294   cout << endl;
 295
 296
 297   //AliAnalysisManager  *man = AliAnalysisManager::GetAnalysisManager();
 298   //AliInputEventHandler* inputHandler = (AliInputEventHandler*)(man->GetInputEventHandler());
 299   //fPIDResponse = (AliPIDResponse*)inputHandler->GetPIDResponse();
 300   //
 301   //fPIDCombined = new AliPIDCombined();
 302   //fPIDCombined->SetSelectedSpecies(AliPID::kSPECIES);
 303   //fPIDCombined->SetDetectorMask(AliPIDResponse::kDetEMCAL);
 304   //fPIDCombined->SetEnablePriors(kFALSE);
 305
 306   double TotalNBins = 0.0;
 307
 308   // Create histograms
 309   Int_t nClustersbins = 501;
 310   Float_t nClusterslow = -0.5, nClustersup = 500.5;
 311   h1_nClusters = new TH1F("h1_nClusters", "# of clusters", nClustersbins, nClusterslow, nClustersup);
 312   h1_nClusters->GetXaxis()->SetTitle("number of clusters/evt");
 313   h1_nClusters->GetYaxis()->SetTitle("counts");
 314   h1_nClusters->SetMarkerStyle(kFullCircle);
 315   TotalNBins+=nClustersbins;
 316
 317   Int_t nZvertexbins = 501;
 318   Float_t Zvertexlow = -50.0, Zvertexup = 50.0;
 319   h1_zvtx = new TH1F("h1_zvtx", "# of clusters", nZvertexbins, Zvertexlow, Zvertexup);
 320   h1_zvtx->GetXaxis()->SetTitle("z_{vertex}");
 321   h1_zvtx->GetYaxis()->SetTitle("counts");
 322   h1_zvtx->SetMarkerStyle(kFullCircle);
 323   TotalNBins+=nZvertexbins;
 324
 325   h1_trigger = new TH1F("h1_trigger", "trigger number returned", 1001,-0.5,1000.5);
 326   TotalNBins+=1001;
 327
 328   Int_t Mbins = 3000;
 329   Float_t Mlow = 0.0, Mup = 3.0;
 330   h1_M = new TH1F("h1_M", "Invariant Mass", Mbins, Mlow, Mup);
 331   h1_M->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 332   h1_M->GetYaxis()->SetTitle("counts");
 333   h1_M->SetMarkerStyle(kFullCircle);
 334   TotalNBins+=Mbins;
 335
 336   h1_M_mix = new TH1F("h1_M_mix", "Invariant Mass (mixed events)", Mbins, Mlow, Mup);
 337   h1_M_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 338   h1_M_mix->GetYaxis()->SetTitle("counts");
 339   h1_M_mix->SetMarkerStyle(kFullCircle);
 340   TotalNBins+=Mbins;
 341
 342   Int_t ptbins = 2000;
 343   Float_t ptlow = 0.0, ptup = 20.0;
 344   Int_t Ebins = 1000;
 345   Float_t Elow = 0.0, Eup = 20.0;
 346   h1_E = new TH1F("h1_E", "Cluster Energy in EMCal", Ebins, Elow, Eup);
 347   h1_E->GetXaxis()->SetTitle("E [GeV]");
 348   h1_E->GetYaxis()->SetTitle("counts");
 349   h1_E->SetMarkerStyle(kFullCircle);
 350   TotalNBins+=Ebins;
 351
 352   h2_PhiEtaCluster = new TH2F("h2_PhiEtaCluster", "cluster phi vs eta", 400,1.362,3.178, 300,-0.728,0.728);
 353   h2_PhiEtaCluster->GetXaxis()->SetTitle("#phi [rad]");
 354   h2_PhiEtaCluster->GetYaxis()->SetTitle("#eta");
 355   h2_PhiEtaCluster->GetZaxis()->SetTitle("hits");
 356   h2_PhiEtaCluster->SetMarkerStyle(kFullCircle);
 357   TotalNBins+=400*300;
 358
 359   h2_PhiEtaClusterCut = new TH2F("h2_PhiEtaClusterCut", "cluster phi vs eta (after cuts)", 400,1.362,3.178, 300,-0.728,0.728);
 360   h2_PhiEtaClusterCut->GetXaxis()->SetTitle("#phi [rad]");
 361   h2_PhiEtaClusterCut->GetYaxis()->SetTitle("#eta");
 362   h2_PhiEtaClusterCut->GetZaxis()->SetTitle("hits");
 363   h2_PhiEtaClusterCut->SetMarkerStyle(kFullCircle);
 364   TotalNBins+=400*300;
 365
 366 // eta binning
 367   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};
 368
 369   // phi binning
 370   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};
 371
 372   h2_PhiEtaMaxCell = new TH2F("h2_PhiEtaMaxCell", "maxcell phi vs eta", 124,PhiBins, 96,EtaBins);
 373   h2_PhiEtaMaxCell->GetXaxis()->SetTitle("#phi [rad]");
 374   h2_PhiEtaMaxCell->GetYaxis()->SetTitle("#eta");
 375   h2_PhiEtaMaxCell->GetZaxis()->SetTitle("hits");
 376   h2_PhiEtaMaxCell->SetMarkerStyle(kFullCircle);
 377   TotalNBins+=96*124;
 378
 379   h2_PhiEtaMaxCellCut = new TH2F("h2_PhiEtaMaxCellCut", "maxcell phi vs eta (after cuts)", 124,PhiBins, 96,EtaBins);
 380   h2_PhiEtaMaxCellCut->GetXaxis()->SetTitle("#phi [rad]");
 381   h2_PhiEtaMaxCellCut->GetYaxis()->SetTitle("#eta");
 382   h2_PhiEtaMaxCellCut->GetZaxis()->SetTitle("hits");
 383   h2_PhiEtaMaxCellCut->SetMarkerStyle(kFullCircle);
 384   TotalNBins+=96*124;
 385
 386   h1_dR_ClustTrk = new TH1F("h1_dR_ClustTrk", "Cluster-Track matching", 5000, -0.01, 5);
 387   h1_dR_ClustTrk->GetXaxis()->SetTitle("dR [sqrt(d#phi^{2}+d#eta^{2})]");
 388   h1_dR_ClustTrk->GetYaxis()->SetTitle("N");
 389   h1_dR_ClustTrk->SetMarkerStyle(kFullCircle);
 390   TotalNBins+=5000;
 391
 392   h2_gE_RecTruth = new TH2F("h2_gE_RecTruth", "#gamma E_{truth}/E_{clust} vs E_{clust}", Ebins,Elow,Eup, 500,0,2);
 393   h2_gE_RecTruth->GetXaxis()->SetTitle("E^{rec}_{clust} [GeV]");
 394   h2_gE_RecTruth->GetYaxis()->SetTitle("E^{rec}_{clust}/E^{truth}_{#gamma}");
 395   h2_gE_RecTruth->GetZaxis()->SetTitle("counts");
 396   h2_gE_RecTruth->SetMarkerStyle(kFullCircle);
 397   TotalNBins+=Ebins*500;
 398
 399   h2_eop_E = new TH2F("h2_eop_E","E/p vs E (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
 400   h2_eop_E->GetXaxis()->SetTitle("cluster Energy [GeV]");
 401   h2_eop_E->GetYaxis()->SetTitle("E/p");
 402   TotalNBins+=Ebins*1200;
 403
 404   h2_eop_pT = new TH2F("h2_eop_pT","E/p vs p_{T} (using built-in track matching)", Ebins, Elow, Eup, 1200,0,3);
 405   h2_eop_pT->GetXaxis()->SetTitle("cluster Energy [GeV]");
 406   h2_eop_pT->GetYaxis()->SetTitle("E/p");
 407   TotalNBins+=Ebins*1200;
 408
 409   h2_E_time = new TH2F("h2_E_time","cluster energy vs time", Ebins, Elow, Eup, 1000,-1e-6,1e-6);
 410   h2_E_time->GetXaxis()->SetTitle("cluster Energy [GeV]");
 411   h2_E_time->GetYaxis()->SetTitle("time [s]");
 412   TotalNBins+=Ebins*1000;
 413
 414   h1_Pi0TruthPt = new TH1F("h1_Pi0TruthPt", "P_{T} distribution for Truth Pi0's", ptbins, ptlow, ptup);
 415   h1_Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 416   h1_Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 417   h1_Pi0TruthPt->SetMarkerStyle(kFullCircle);
 418   TotalNBins+=ptbins;
 419
 420   h1_K0Pi0TruthPt = new TH1F("h1_K0Pi0TruthPt", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays", ptbins, ptlow, ptup);
 421   h1_K0Pi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 422   h1_K0Pi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 423   h1_K0Pi0TruthPt->SetMarkerStyle(kFullCircle);
 424   TotalNBins+=ptbins;
 425
 426   h1_PriPi0TruthPt = new TH1F("h1_PriPi0TruthPt", "P_{T} distribution for Truth Primary Pi0's", ptbins, ptlow, ptup);
 427   h1_PriPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 428   h1_PriPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 429   h1_PriPi0TruthPt->SetMarkerStyle(kFullCircle);
 430   TotalNBins+=ptbins;
 431
 432   h1_PhysPi0TruthPt = new TH1F("h1_PhysPi0TruthPt", "P_{T} distribution for Truth Physical Primary Pi0's", ptbins, ptlow, ptup);
 433   h1_PhysPi0TruthPt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 434   h1_PhysPi0TruthPt->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 435   h1_PhysPi0TruthPt->SetMarkerStyle(kFullCircle);
 436   TotalNBins+=ptbins;
 437
 438   h1_Pi0TruthPtEmcal = new TH1F("h1_Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's (hit EMCal)", ptbins, ptlow, ptup);
 439   h1_Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 440   h1_Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 441   h1_Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 442   TotalNBins+=ptbins;
 443
 444   h1_K0Pi0TruthPtEmcal = new TH1F("h1_K0Pi0TruthPtEmcal", "P_{T} distribution for Truth Pi0's from K^{0}_{s} decays (hit EMCal)", ptbins, ptlow, ptup);
 445   h1_K0Pi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 446   h1_K0Pi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 447   h1_K0Pi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 448   TotalNBins+=ptbins;
 449
 450   h1_PriPi0TruthPtEmcal = new TH1F("h1_PriPi0TruthPtEmcal", "P_{T} distribution for Truth Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
 451   h1_PriPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 452   h1_PriPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 453   h1_PriPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 454   TotalNBins+=ptbins;
 455
 456   h1_PhysPi0TruthPtEmcal = new TH1F("h1_PhysPi0TruthPtEmcal", "P_{T} distribution for Truth Physical Primary Pi0's (hit EMCal)", ptbins, ptlow, ptup);
 457   h1_PhysPi0TruthPtEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 458   h1_PhysPi0TruthPtEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 459   h1_PhysPi0TruthPtEmcal->SetMarkerStyle(kFullCircle);
 460   TotalNBins+=ptbins;
 461
 462   h1_Pi0TruthPtPhi2piEta065 = new TH1F("h1_Pi0TruthPtPhi2piEta065",
 463                                        "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 464   h1_Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 465   h1_Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 466   h1_Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
 467   TotalNBins+=ptbins;
 468
 469   h1_K0Pi0TruthPtPhi2piEta065 = new TH1F("h1_K0Pi0TruthPtPhi2piEta065",
 470                                          "P_{T} for Truth Pi0's (K^{0}_{s} decays) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 471   h1_K0Pi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 472   h1_K0Pi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 473   h1_K0Pi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
 474   TotalNBins+=ptbins;
 475
 476   h1_PriPi0TruthPtPhi2piEta065 = new TH1F("h1_PriPi0TruthPtPhi2piEta065",
 477                                           "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 478   h1_PriPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 479   h1_PriPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 480   h1_PriPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
 481   TotalNBins+=ptbins;
 482
 483   h1_PhysPi0TruthPtPhi2piEta065 = new TH1F("h1_PhysPi0TruthPtPhi2piEta065",
 484                                            "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<0.65 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 485   h1_PhysPi0TruthPtPhi2piEta065->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 486   h1_PhysPi0TruthPtPhi2piEta065->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 487   h1_PhysPi0TruthPtPhi2piEta065->SetMarkerStyle(kFullCircle);
 488   TotalNBins+=ptbins;
 489
 490   h1_Pi0TruthPtPhi2piEta1 = new TH1F("h1_Pi0TruthPtPhi2piEta1",
 491                                      "P_{T} for Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 492   h1_Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 493   h1_Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 494   h1_Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
 495   TotalNBins+=ptbins;
 496
 497   h1_K0Pi0TruthPtPhi2piEta1 = new TH1F("h1_K0Pi0TruthPtPhi2piEta1",
 498                                        "P_{T} for Truth Pi0's (k^{0}_{s} decays) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 499   h1_K0Pi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 500   h1_K0Pi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 501   h1_K0Pi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
 502   TotalNBins+=ptbins;
 503
 504   h1_PriPi0TruthPtPhi2piEta1 = new TH1F("h1_PriPi0TruthPtPhi2piEta1",
 505                                         "P_{T} for Primary Truth Pi0's [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 506   h1_PriPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 507   h1_PriPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 508   h1_PriPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
 509   TotalNBins+=ptbins;
 510
 511   h1_PhysPi0TruthPtPhi2piEta1 = new TH1F("h1_PhysPi0TruthPtPhi2piEta1",
 512                                      "P_{T} for Truth Pi0's (not from material) [|#eta_{#pi^{0}}|<1.0 && 0<#phi_{#pi^{0}}<2#pi]", ptbins, ptlow, ptup);
 513   h1_PhysPi0TruthPtPhi2piEta1->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 514   h1_PhysPi0TruthPtPhi2piEta1->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 515   h1_PhysPi0TruthPtPhi2piEta1->SetMarkerStyle(kFullCircle);
 516   TotalNBins+=ptbins;
 517
 518   h2_Pi0TruthPhiEta = new TH2F("h2_Pi0TruthPhiEta","Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 519   h2_Pi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 520   h2_Pi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 521   TotalNBins+=380*200;
 522
 523   h2_PriPi0TruthPhiEta = new TH2F("h2_PriPi0TruthPhiEta","Primary Pi0Truth Phi vs Eta ", 380,-0.02,6.30, 200,-10,10);
 524   h2_PriPi0TruthPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 525   h2_PriPi0TruthPhiEta->GetYaxis()->SetTitle("#eta ");
 526   TotalNBins+=380*200;
 527
 528   h2_Pi0TruthPhiEtaEmcal = new TH2F("h2_Pi0TruthPhiEtaEmcal","Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
 529   h2_Pi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 530   h2_Pi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 531   TotalNBins+=380*150;
 532
 533   h2_PriPi0TruthPhiEtaEmcal = new TH2F("h2_PriPi0TruthPhiEtaEmcal","Primary Pi0Truth Phi vs Eta (in EMCal)", 380,-0.02,6.30, 150,-5,5);
 534   h2_PriPi0TruthPhiEtaEmcal->GetXaxis()->SetTitle("#phi [rad]");
 535   h2_PriPi0TruthPhiEtaEmcal->GetYaxis()->SetTitle("#eta ");
 536   TotalNBins+=380*150;
 537
 538   h1_TruthPhotonsEmcal = new TH1F("h1_TruthPhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
 539   h1_TruthPhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 540   h1_TruthPhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 541   h1_TruthPhotonsEmcal->SetMarkerStyle(kFullCircle);
 542     TotalNBins+=ptbins;
 543
 544   h2_TruthPhotonsPhiEta = new TH2F("h2_TruthPhotonsPhiEta",
 545                                    "Truth Photons Phi vs Eta (pointed at emcal)", 380,-0.02,6.30, 150,-1.5,1.5);
 546   h2_TruthPhotonsPhiEta->GetXaxis()->SetTitle("#phi [rad]");
 547   h2_TruthPhotonsPhiEta->GetYaxis()->SetTitle("#eta ");
 548   h2_TruthPhotonsPhiEta->SetMarkerStyle(kFullCircle);
 549   TotalNBins+=380*150;
 550
 551   h1_PhotonsEmcal = new TH1F("h1_PhotonsEmcal", "P_{T} distribution for photons (in EMCal)", ptbins, ptlow, ptup);
 552   h1_PhotonsEmcal->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 553   h1_PhotonsEmcal->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 554   h1_PhotonsEmcal->SetMarkerStyle(kFullCircle);
 555   TotalNBins+=ptbins;
 556
 557   h1_PhotonsNCellsCut = new TH1F("h1_PhotonsNCellsCut", "P_{T} distribution for #gamma's that survive NCells cut", ptbins, ptlow, ptup);
 558   h1_PhotonsNCellsCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 559   h1_PhotonsNCellsCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 560   h1_PhotonsNCellsCut->SetMarkerStyle(kFullCircle);
 561   TotalNBins+=ptbins;
 562
 563   h1_PhotonsTrackMatchCut = new TH1F("h1_PhotonsTrackMatchCut", "P_{T} distribution for #gamma's that survive TrackMatch cut", ptbins, ptlow, ptup);
 564   h1_PhotonsTrackMatchCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 565   h1_PhotonsTrackMatchCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 566   h1_PhotonsTrackMatchCut->SetMarkerStyle(kFullCircle);
 567   TotalNBins+=ptbins;
 568
 569   h1_PhotonsAllCut = new TH1F("h1_PhotonsAllCut", "P_{T} distribution for #gamma's that survive All cut", ptbins, ptlow, ptup);
 570   h1_PhotonsAllCut->GetXaxis()->SetTitle("P_{T} (GeV/c)");
 571   h1_PhotonsAllCut->GetYaxis()->SetTitle("dN/dP_{T} (c/GeV)");
 572   h1_PhotonsAllCut->SetMarkerStyle(kFullCircle);
 573   TotalNBins+=ptbins;
 574
 575   h2_PhotonsPhiEtaIsEmcal = new TH2F("h2_PhotonsPhiEtaIsEmcal",
 576                                      "Photons Phi vs Eta (IsEMCAL()==1)", 380,-0.02,6.30, 100,-1.0,1.0);
 577   h2_PhotonsPhiEtaIsEmcal->GetXaxis()->SetTitle("#phi [rad]");
 578   h2_PhotonsPhiEtaIsEmcal->GetYaxis()->SetTitle("#eta ");
 579   TotalNBins+=380*100;
 580
 581   h1_dR_RealMC = new TH1F("h1_dR_RealMC", "P_{T} distribution for #gamma's that survive All cut", 2000, -0.01, 10);
 582   h1_dR_RealMC->GetXaxis()->SetTitle("dR sqrt(dx^{2}+dy^{2})");
 583   h1_dR_RealMC->GetYaxis()->SetTitle("N");
 584   h1_dR_RealMC->SetMarkerStyle(kFullCircle);
 585   TotalNBins+=2000;
 586
 587   h2_Mpt_Pri = new TH2F("h2_Mpt_Pri", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
 588   h2_Mpt_Pri->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 589   h2_Mpt_Pri->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 590   h2_Mpt_Pri->SetMarkerStyle(kFullCircle);
 591   TotalNBins+=Mbins*ptbins;
 592
 593   h2_Mpt_Sec = new TH2F("h2_Mpt_Sec", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
 594   h2_Mpt_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 595   h2_Mpt_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 596   h2_Mpt_Sec->SetMarkerStyle(kFullCircle);
 597   TotalNBins+=Mbins*ptbins;
 598
 599   h3_MptR_Sec = new TH3F("h3_MptR_Sec", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600);
 600   h3_MptR_Sec->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 601   h3_MptR_Sec->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 602   h3_MptR_Sec->GetZaxis()->SetTitle("production radius [cm]");
 603   h3_MptR_Sec->SetMarkerStyle(kFullCircle);
 604   TotalNBins+=500*100*300;
 605
 606   h3_MptR_K0s = new TH3F("h3_MptR_K0s", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600);
 607   h3_MptR_K0s->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 608   h3_MptR_K0s->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 609   h3_MptR_K0s->GetZaxis()->SetTitle("production radius [cm]");
 610   h3_MptR_K0s->SetMarkerStyle(kFullCircle);
 611   TotalNBins+=500*100*300;
 612
 613   h3_MptR_Mat = new TH3F("h3_MptR_Mat", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600);
 614   h3_MptR_Mat->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 615   h3_MptR_Mat->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 616   h3_MptR_Mat->GetZaxis()->SetTitle("production radius [cm]");
 617   h3_MptR_Mat->SetMarkerStyle(kFullCircle);
 618   TotalNBins+=500*100*300;
 619
 620   h2_PtR_MatM = new TH2F("h2_PtR_MatM", "pT vs production radius for merged material pions (pi mass assumed)", 100,0,20, 300,0,600);
 621   h2_PtR_MatM->GetXaxis()->SetTitle("p_{T} [GeV/c]");
 622   h2_PtR_MatM->GetYaxis()->SetTitle("production radius [cm]");
 623   h2_PtR_MatM->SetMarkerStyle(kFullCircle);
 624   TotalNBins+=100*300;
 625
 626   h2_Mpt_Pri_conv = new TH2F("h2_Mpt_Pri_conv", "mass vs pT for primary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
 627   h2_Mpt_Pri_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 628   h2_Mpt_Pri_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 629   h2_Mpt_Pri_conv->SetMarkerStyle(kFullCircle);
 630   TotalNBins+=Mbins*ptbins;
 631
 632   h2_Mpt_Sec_conv = new TH2F("h2_Mpt_Sec_conv", "mass vs pT for secondary pions", Mbins, Mlow, Mup, ptbins, ptlow, ptup);
 633   h2_Mpt_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 634   h2_Mpt_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 635   h2_Mpt_Sec_conv->SetMarkerStyle(kFullCircle);
 636   TotalNBins+=Mbins*ptbins;
 637
 638   h3_MptR_Sec_conv = new TH3F("h3_MptR_Sec_conv", "mass vs pT vs production radius for secondary pions", 500,0,0.5, 100,0,20, 300,0,600);
 639   h3_MptR_Sec_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 640   h3_MptR_Sec_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 641   h3_MptR_Sec_conv->GetZaxis()->SetTitle("production radius [cm]");
 642   h3_MptR_Sec_conv->SetMarkerStyle(kFullCircle);
 643   TotalNBins+=500*100*300;
 644
 645   h3_MptR_K0s_conv = new TH3F("h3_MptR_K0s_conv", "mass vs pT vs production radius for K0s pions", 500,0,0.5, 100,0,20, 300,0,600);
 646   h3_MptR_K0s_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 647   h3_MptR_K0s_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 648   h3_MptR_K0s_conv->GetZaxis()->SetTitle("production radius [cm]");
 649   h3_MptR_K0s_conv->SetMarkerStyle(kFullCircle);
 650   TotalNBins+=500*100*300;
 651
 652   h3_MptR_Mat_conv = new TH3F("h3_MptR_Mat_conv", "mass vs pT vs production radius for material pions", 500,0,0.5, 100,0,20, 300,0,600);
 653   h3_MptR_Mat_conv->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 654   h3_MptR_Mat_conv->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 655   h3_MptR_Mat_conv->GetZaxis()->SetTitle("production radius [cm]");
 656   h3_MptR_Mat_conv->SetMarkerStyle(kFullCircle);
 657   TotalNBins+=500*100*300;
 658
 659   h1_eConversionR = new TH1F("h1_eConversionR", "conversion point (radius)", 600,0,600);
 660   h1_eConversionR->GetXaxis()->SetTitle("production radius [cm]");
 661   h1_eConversionR->SetMarkerStyle(kFullCircle);
 662   TotalNBins+=600;
 663
 664   h1_PriPi0Mother = new TH1F("h1_PriPi0Mother", "primary pi0 mother ID", 12001,-6000.5,6000.5);
 665   h1_PriPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID");
 666   h1_PriPi0Mother->SetMarkerStyle(kFullCircle);
 667   TotalNBins+=12001;
 668
 669   h1_SecPi0Mother = new TH1F("h1_SecPi0Mother", "secondray pi0 mother ID", 12001,-6000.5,6000.5);
 670   h1_SecPi0Mother->GetXaxis()->SetTitle("#pi^{0} mother ID");
 671   h1_SecPi0Mother->SetMarkerStyle(kFullCircle);
 672   TotalNBins+=12001;
 673
 674   Int_t chi2bins = 100;
 675   Float_t chi2low = -2, chi2up = 2;
 676   h1_Chi2 = new TH1F("h1_Chi2","#chi^{2} distribution for reconstructed",chi2bins, chi2low, chi2up);
 677   h1_Chi2->GetXaxis()->SetTitle("#chi^{2}");
 678   h1_Chi2->GetYaxis()->SetTitle("counts");
 679   TotalNBins+=chi2bins;
 680
 681   h1_nTrkMatch = new TH1F("h1_nTrkMatch","number of matched tracks",14, -1.5, 5.5);
 682   h1_nTrkMatch->GetXaxis()->SetTitle("nTracksMatched");
 683   h1_nTrkMatch->GetYaxis()->SetTitle("counts");
 684   TotalNBins+=14;
 685
 686   h1_ClusterDisp = new TH1F("h1_ClusterDisp","Dispersion of CaloCluster",1000, -1, 3);
 687   h1_ClusterDisp->GetXaxis()->SetTitle("cluster->GetClusterDisp()");
 688   h1_ClusterDisp->GetYaxis()->SetTitle("counts");
 689   TotalNBins+=1000;
 690
 691   h2_Ellipse = new TH2F("h2_Ellipse","Ellipse axis M20 vs M02",500, -0.01, 1, 500, -0.01, 1);
 692   h2_Ellipse->GetXaxis()->SetTitle("cluster->GetM20()");
 693   h2_Ellipse->GetYaxis()->SetTitle("cluster->GetM02()");
 694   h2_Ellipse->GetZaxis()->SetTitle("counts");
 695   TotalNBins+=500*500;
 696
 697   Int_t etabins = 150;
 698   Float_t etalow = -1.5, etaup = 1.5;
 699   h2_EtaPt = new TH2F("h2_EtaPt","Cluster Energy vs ",etabins, etalow, etaup, ptbins, ptlow, ptup);
 700   h2_EtaPt->GetXaxis()->SetTitle("E [GeV]");
 701   h2_EtaPt->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 702   TotalNBins+=etabins*ptbins;
 703
 704   h3_MptAsymm = new TH3F("h3_MptAsymm","mass vs p_{T} vs Asymm cut",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
 705   h3_MptAsymm->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 706   h3_MptAsymm->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 707   h3_MptAsymm->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
 708   TotalNBins+=Mbins*ptbins*3.0;
 709
 710   h3_MptAsymm_mix = new TH3F("h3_MptAsymm_mix","mass vs p_{T} vs Asymm cut (mixed events)",Mbins,Mlow,Mup, ptbins,ptlow,ptup, 3,0.5,3.5);
 711   h3_MptAsymm_mix->GetXaxis()->SetTitle("mass [GeV/c^{2}]");
 712   h3_MptAsymm_mix->GetYaxis()->SetTitle("p_{T} [GeV/c]");
 713   h3_MptAsymm_mix->GetZaxis()->SetTitle("Asymmetry Cut (edges: 0.0, 0.1, 0.7, 1.0)");
 714   TotalNBins+=Mbins*ptbins*3.0;
 715
 716   h2_dphi_deta = new TH2F("h2_dphi_deta","#Delta#phi vs #Delta#eta", 349,-1.5,5, 400,-2.0,2.0);
 717   h2_dphi_deta->GetXaxis()->SetTitle("#Delta#phi");
 718   h2_dphi_deta->GetYaxis()->SetTitle("#Delta#eta");
 719   TotalNBins+=349*400;
 720
 721   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);
 722   h2_dphi_deta_mix->GetXaxis()->SetTitle("#Delta#phi");
 723   h2_dphi_deta_mix->GetYaxis()->SetTitle("#Delta#eta");
 724   TotalNBins+=349*400;
 725
 726   h2_DispRes = new TH2F("h2_DispRes", "zvtx info", 500,-0.01,1, 500,-0.1,2);
 727   h2_DispRes->GetXaxis()->SetTitle("EvtVtx->GetDispersion()");
 728   h2_DispRes->GetYaxis()->SetTitle("EvtVtx->GetZRes()");
 729   h2_DispRes->GetZaxis()->SetTitle("counts");
 730   TotalNBins+=500*500;
 731
 732   h2_cells_M02 = new TH2F("h2_cells_M02", "nCells vs M02", 204,-1.5,100.5, 500,-1,1.5);
 733   h2_cells_M02->GetXaxis()->SetTitle("nCells");
 734   h2_cells_M02->GetYaxis()->SetTitle("M02");
 735   h2_cells_M02->GetZaxis()->SetTitle("counts");
 736   TotalNBins+=204*500;
 737
 738   cout << endl << "Total number of bins in booked histograms:  " << TotalNBins << endl << endl;
 739
 740   // Initialize helper class (for vertex selection & pile up correction)
 741   fHelperClass = new AliAnalysisUtils();
 742
 743   //TFile *f = OpenFile(1);
 744   //TDirectory::TContext context(f);
 745
 746   fOutput->Add(h1_nClusters);
 747   fOutput->Add(h1_zvtx);
 748   fOutput->Add(h1_trigger);
 749   fOutput->Add(h1_M);
 750   fOutput->Add(h1_M_mix);
 751   fOutput->Add(h1_E);
 752   fOutput->Add(h2_PhiEtaCluster);
 753   fOutput->Add(h2_PhiEtaClusterCut);
 754   fOutput->Add(h2_PhiEtaMaxCell);
 755   fOutput->Add(h2_PhiEtaMaxCellCut);
 756   fOutput->Add(h1_dR_ClustTrk);
 757   fOutput->Add(h2_gE_RecTruth);
 758   fOutput->Add(h2_eop_E);
 759   fOutput->Add(h2_eop_pT);
 760   fOutput->Add(h2_E_time);
 761   fOutput->Add(h1_Pi0TruthPt);
 762   fOutput->Add(h1_K0Pi0TruthPt);
 763   fOutput->Add(h1_PriPi0TruthPt);
 764   fOutput->Add(h1_PhysPi0TruthPt);
 765   fOutput->Add(h1_Pi0TruthPtEmcal);
 766   fOutput->Add(h1_K0Pi0TruthPtEmcal);
 767   fOutput->Add(h1_PriPi0TruthPtEmcal);
 768   fOutput->Add(h1_PhysPi0TruthPtEmcal);
 769   fOutput->Add(h1_Pi0TruthPtPhi2piEta065);
 770   fOutput->Add(h1_K0Pi0TruthPtPhi2piEta065);
 771   fOutput->Add(h1_PriPi0TruthPtPhi2piEta065);
 772   fOutput->Add(h1_PhysPi0TruthPtPhi2piEta065);
 773   fOutput->Add(h1_Pi0TruthPtPhi2piEta1);
 774   fOutput->Add(h1_K0Pi0TruthPtPhi2piEta1);
 775   fOutput->Add(h1_PriPi0TruthPtPhi2piEta1);
 776   fOutput->Add(h1_PhysPi0TruthPtPhi2piEta1);
 777   fOutput->Add(h2_Pi0TruthPhiEta);
 778   fOutput->Add(h2_PriPi0TruthPhiEta);
 779   fOutput->Add(h2_Pi0TruthPhiEtaEmcal);
 780   fOutput->Add(h2_PriPi0TruthPhiEtaEmcal);
 781   fOutput->Add(h1_TruthPhotonsEmcal);
 782   fOutput->Add(h2_TruthPhotonsPhiEta);
 783   fOutput->Add(h1_PhotonsEmcal);
 784   fOutput->Add(h1_PhotonsNCellsCut);
 785   fOutput->Add(h1_PhotonsTrackMatchCut);
 786   fOutput->Add(h1_PhotonsAllCut);
 787   fOutput->Add(h2_PhotonsPhiEtaIsEmcal);
 788   fOutput->Add(h1_dR_RealMC);
 789   fOutput->Add(h2_Mpt_Pri);
 790   fOutput->Add(h2_Mpt_Sec);
 791   fOutput->Add(h3_MptR_Sec);
 792   fOutput->Add(h3_MptR_K0s);
 793   fOutput->Add(h3_MptR_Mat);
 794   fOutput->Add(h2_PtR_MatM);
 795   fOutput->Add(h2_Mpt_Pri_conv);
 796   fOutput->Add(h2_Mpt_Sec_conv);
 797   fOutput->Add(h3_MptR_Sec_conv);
 798   fOutput->Add(h3_MptR_K0s_conv);
 799   fOutput->Add(h3_MptR_Mat_conv);
 800   fOutput->Add(h1_eConversionR);
 801   fOutput->Add(h1_PriPi0Mother);
 802   fOutput->Add(h1_SecPi0Mother);
 803   fOutput->Add(h1_Chi2);
 804   fOutput->Add(h1_nTrkMatch);
 805   fOutput->Add(h1_ClusterDisp);
 806   fOutput->Add(h2_Ellipse);
 807   fOutput->Add(h2_EtaPt);
 808   fOutput->Add(h3_MptAsymm);
 809   fOutput->Add(h3_MptAsymm_mix);
 810   fOutput->Add(h2_dphi_deta);
 811   fOutput->Add(h2_dphi_deta_mix);
 812   fOutput->Add(h2_DispRes);
 813   fOutput->Add(h2_cells_M02);
 814
 815   // Post data for ALL output slots >0 here,
 816   // To get at least an empty histogram
 817   // 1 is the outputnumber of a certain weg of task 1
 818   PostData(1, fOutput);
 819 }
 820
 821 //________________________________________________________________________
 822 void AliAnalysisTaskEMCALMesonGGSDM::UserExec(Option_t *)
 823 {
 824   // Main loop Called for each event
 825
 826   AliMCEvent *mcEvent = MCEvent();
 827   Bool_t isMC = bool(mcEvent);//is this the right way to do this?
 828
 829   TRandom3 randy; randy.SetSeed(0);
 830   unsigned int iskip = -1;
 831   TLorentzVector ParentMix;
 832
 833   double recalScale = 1.0;
 834
 835   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 836
 837   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (am->GetInputEventHandler());
 838   AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (am->GetInputEventHandler());
 839   if (!aodH && !esdH)  Printf("ERROR: Could not get ESD or AODInputHandler");
 840
 841   if(esdH)      fEsdEv = esdH->GetEvent();
 842   else if(aodH) fAodEv = aodH->GetEvent();
 843   else{
 844     AliFatal("Neither ESD nor AOD event found");
 845     return;
 846   }
 847
 848   // get pointer to reconstructed event
 849   AliVEvent *event = InputEvent();
 850   if (!event){
 851     AliError("Pointer == 0, this can not happen!");  return;}
 852   //AliESDEvent* fEsdEv = dynamic_cast<AliESDEvent*>(event);
 853   //AliAODEvent* aod = dynamic_cast<AliAODEvent*>(event);
 854   //if (!fEsdEv){
 855   //AliError("Cannot get the ESD event");  return;}
 856
 857   /*
 858   fHelperClass->SetCutOnZVertexSPD(kFALSE);//does the zvtx have to match the spd vertex?
 859   fHelperClass->SetMaxVtxZ(1.0e6);//i set this myself later..
 860   // simply makes sure that there is at least 1 contributer to the zvtx determination.
 861   // this should only remove the *extra* events at zvtx==0.
 862   if(!fHelperClass->IsVertexSelected2013pA(event))
 863     return;
 864   */
 865
 866   Int_t iTrigger = 0;
 867   if (fEsdEv)       iTrigger = fEsdEv->GetHeader()->GetL0TriggerInputs();
 868   else if (fAodEv)  iTrigger = fAodEv->GetHeader()->GetL0TriggerInputs();
 869
 870   char saythis[500];
 871   Int_t iTriggerBin = 0;
 872   for(unsigned long j=0; j<TriggerList.size(); j++){
 873     if(iTrigger==TriggerList[j])
 874       iTriggerBin=j+1;
 875   }
 876   if(iTriggerBin==0){
 877     TriggerList.push_back(iTrigger);
 878     iTriggerBin=TriggerList.size();
 879   }
 880
 881   h1_trigger->SetBinContent(iTriggerBin, h1_trigger->GetBinContent(iTriggerBin)+1);
 882   sprintf(saythis,"%d",iTrigger);
 883   h1_trigger->GetXaxis()->SetBinLabel(iTriggerBin, saythis);
 884
 885   if(fEsdEv){
 886     TString trigClasses = fEsdEv->GetFiredTriggerClasses();
 887     // remove "fast cluster events":
 888     if (trigClasses.Contains("FAST")  && !trigClasses.Contains("ALL"))
 889       return;
 890   }
 891   else if(fAodEv){
 892     TString trigClasses = fAodEv->GetFiredTriggerClasses();
 893     // remove "fast cluster events":
 894     if (trigClasses.Contains("FAST")  && !trigClasses.Contains("ALL"))
 895       return;
 896   }
 897
 898   if (fEsdEv){
 899     if(!(fEsdEv->GetPrimaryVertex()->GetStatus()))   return;
 900   }
 901   //else if (fAodEv){
 902   //if(!(fAodEv->GetPrimaryVertex()->GetStatus()))   return;
 903   //}
 904
 905   Double_t vertDisp=0.0;
 906   Double_t vertZres=0.0;
 907   Bool_t vertIsfromZ=0;
 908   if (fEsdEv){
 909     vertDisp    = fEsdEv->GetPrimaryVertex()->GetDispersion();
 910     vertZres    = fEsdEv->GetPrimaryVertex()->GetZRes();
 911     vertIsfromZ = fEsdEv->GetPrimaryVertex()->IsFromVertexerZ();
 912   }
 913   else if (fAodEv){
 914     vertDisp    = 0;
 915     vertZres    = 0;
 916     vertIsfromZ = 0;
 917   }
 918
 919   h2_DispRes->Fill(vertDisp, vertZres);
 920   // if vertex is from spd vertexZ, require more stringent cut
 921   if (vertIsfromZ) {
 922     if (vertDisp>0.02 ||  vertZres>0.25 )
 923       return; // bad vertex from VertexerZ
 924   }
 925
 926
 927   Int_t nclusters=0;
 928   if(fEsdEv){
 929     //Int_t evtN      = fEsdEv->GetEventNumberInFile();
 930     //Int_t ntracks   = fEsdEv->GetNumberOfTracks();
 931     nclusters = fEsdEv->GetNumberOfCaloClusters();
 932   }
 933   else if(fAodEv){
 934     //Int_t evtN      = fAodEv->GetEventNumberInFile();
 935     //Int_t ntracks   = fAodEv->GetNumberOfTracks();
 936     nclusters = fAodEv->GetNumberOfCaloClusters();
 937   }
 938
 939   // EMCal cluster loop for reconstructed event
 940   //numberofclusters set above!
 941   TLorentzVector Photon1, Photon2, Parent;
 942   Double_t vertex[3];
 943   Double_t E1=0.0;
 944   Double_t E2=0.0;
 945   Double_t vertZ=0.0;
 946   if (fEsdEv)       vertZ = fEsdEv->GetPrimaryVertex()->GetZ();
 947   else if (fAodEv)  vertZ = fAodEv->GetPrimaryVertex()->GetZ();
 948
 949   h1_zvtx->Fill(vertZ);
 950   //zvertex cut:
 951   if(fabs(vertZ)>10.0)
 952     return;
 953
 954   h1_nClusters->Fill(nclusters);
 955
 956   int izvtx = GetZvtxBin(vertZ);
 957   int imult = GetMultBin(nclusters);
 958
 959   //cout << iskip << " " << izvtx << " " << imult << endl;
 960   //cout << "GetNumberOfVertices(): " << fAodEv->GetNumberOfVertices() << endl;
 961
 962
 963
 964   //######################### ~~~~~~~~~~~ ##################################
 965   //######################### STARTING MC ##################################
 966   //######################### ~~~~~~~~~~~ ##################################
 967
 968   if(isMC){
 969     int isPrimary     = 0;
 970     int isMaterialSec    = 0;
 971     int isK0sDecay    = 0;
 972
 973     if (!mcEvent){
 974       cout << "no MC event" << endl;
 975       return;
 976     }
 977
 978     const AliVVertex *evtVtx = mcEvent->GetPrimaryVertex();
 979     if (!evtVtx)
 980       return;
 981
 982     mcEvent->PreReadAll();
 983
 984     Int_t nTracksMC  = mcEvent->GetNumberOfTracks();
 985     Int_t nPTracksMC = mcEvent->GetNumberOfPrimaries();
 986     Int_t MyEMCPion = -1;
 987     for (Int_t iTrack = 0; iTrack<nTracksMC; ++iTrack) {
 988       AliMCParticle *mcP = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
 989       if (!mcP)
 990         continue;
 991
 992       if(iTrack<nPTracksMC)  isPrimary = 1;
 993       else                   isPrimary = 0;
 994
 995       isK0sDecay = 0;
 996       if(mcP->GetMother()>-1){
 997         if( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  310 ||
 998             ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -310  )
 999           isK0sDecay = 1;
1000       }
1001
1002       // it's a pion !!
1003       if(mcP->PdgCode() != 111)
1004         continue;
1005
1006       MyEMCPion = 0;
1007       if(strcmp(fMyMCType,"ANY")==0)
1008         MyEMCPion = 1;
1009       else if(IsMyMCHeaderType(iTrack, fMyMCType, mcEvent)){
1010         //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << "   i: " << iTrack << "    pdg: " << mcP->PdgCode() << "   pT: " << mcP->Pt() << endl;
1011         //cout << "iTrack: " << iTrack << "   nPrimaryMC: " << nPTracksMC << "   nTracksMC: " << nTracksMC << endl;
1012         MyEMCPion = 1;
1013       }
1014
1015       //if(MyEMCPion)
1016       //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << "   i: " << iTrack << "    pdg: " << mcP->PdgCode() << "   pT: " << mcP->Pt() << endl;
1017       //cout << "evtN: " << fEsdEv->GetEventNumberInFile() << "   i: " << iTrack << "    pdg: " << static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack+1))->PdgCode() << "   pT: " << mcP->Pt() << endl;
1018
1019       if(MyEMCPion!=1 && isPrimary==1)
1020         continue;
1021
1022
1023       if(isPrimary==1 && mcP->GetMother()>-1)
1024         h1_PriPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() );
1025       else if(isPrimary==0 && mcP->GetMother()>-1)
1026         h1_SecPi0Mother->Fill( ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() );
1027
1028       Int_t daughter[2] = {-1,-1};
1029       daughter[0] = mcP->GetFirstDaughter();
1030       daughter[1] = mcP->GetLastDaughter();
1031
1032       if (daughter[0]<0)  continue;
1033       if (daughter[1]<0)  daughter[1]=daughter[0];
1034       if (daughter[1]-daughter[0] != 1)  continue;
1035
1036       Int_t eIndexofConvertedPhoton[2] = {-1,-1};
1037
1038       bool bacc = true;
1039       bool binp = true;
1040       Double_t eta_d[2] = {0.0,0.0};
1041       Double_t phi_d[2] = {0.0,0.0};
1042       for (Int_t daughter_index=0; daughter_index<2; daughter_index++){
1043         const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index]));
1044         eta_d[daughter_index] = dmc->Eta();
1045         phi_d[daughter_index] = dmc->Phi();
1046         if(!(dmc->PdgCode()==22))         binp = false;
1047         if(!(dmc->PdgCode()==22 &&
1048              eta_d[daughter_index]>fEtamin && eta_d[daughter_index]<fEtamax &&
1049              phi_d[daughter_index]>fPhimin && phi_d[daughter_index]<fPhimax))   bacc = false;
1050
1051         if(dmc->GetFirstDaughter()>0 && dmc->GetLastDaughter()>0) {
1052           // get the photons's daughters...
1053           const AliMCParticle *dmcd1 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetFirstDaughter()));
1054           const AliMCParticle *dmcd2 = static_cast<const AliMCParticle *>(mcEvent->GetTrack(dmc->GetLastDaughter()));
1055           Double_t productionR1 = TMath::Sqrt(dmcd1->Xv()*dmcd1->Xv() + dmcd1->Yv()*dmcd1->Yv());
1056           if(bacc)  h1_eConversionR->Fill(productionR1);
1057           // check if this is a conversion...
1058           if( (dmcd1->PdgCode()== -1.0*dmcd2->PdgCode()) &&
1059               (dmcd1->PdgCode()==11 || dmcd1->PdgCode()==-11) &&
1060               productionR1<440.0){
1061             //find the conv e with highest energy, assign it to be that photon decay product.
1062             if( dmcd1->E() > dmcd2->E() )
1063               eIndexofConvertedPhoton[daughter_index] = dmc->GetFirstDaughter();
1064             else
1065               eIndexofConvertedPhoton[daughter_index] = dmc->GetLastDaughter();
1066           }
1067         }
1068       }
1069
1070       if(binp!=true)
1071         continue;
1072
1073       // primary particle
1074       //Double_t dR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1075       //                          (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1076       //if(dR <= 0.01)  isPrimary = 1;
1077       //else            isPrimary = 0;
1078
1079       isMaterialSec = 0;
1080       if(isPrimary!=1){
1081         if(((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  2212 || //proton
1082            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2212 || //anti-proton
1083            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  2112 || //neutron
1084            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -2112 || //anti-neutron
1085            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  321  || //K+
1086            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -321  || //K-
1087            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() ==  211  || //pi+
1088            ((AliMCParticle*)mcEvent->GetTrack(mcP->GetMother()))->PdgCode() == -211     //pi-
1089            )
1090           isMaterialSec = 1;
1091       }
1092
1093       h1_Pi0TruthPt                  ->Fill(mcP->Pt());
1094       if(isK0sDecay)  h1_K0Pi0TruthPt->Fill(mcP->Pt());
1095       h2_Pi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
1096
1097       if(isPrimary==1){
1098         h1_PriPi0TruthPt    ->Fill(mcP->Pt());
1099         h2_PriPi0TruthPhiEta->Fill(mcP->Phi(),mcP->Eta());
1100       }
1101       if(isPrimary!=1 && isMaterialSec!=1)   h1_PhysPi0TruthPt->Fill(mcP->Pt());
1102
1103       if(mcP->Eta()<-1.0 || mcP->Eta()>1.0)
1104         continue;
1105
1106       h1_Pi0TruthPtPhi2piEta1         ->Fill(mcP->Pt());
1107       if(isPrimary==1)
1108         h1_PriPi0TruthPtPhi2piEta1    ->Fill(mcP->Pt());
1109       if(isK0sDecay)
1110         h1_K0Pi0TruthPtPhi2piEta1     ->Fill(mcP->Pt());
1111       if(isPrimary!=1 && isMaterialSec!=1)
1112         h1_PhysPi0TruthPtPhi2piEta1   ->Fill(mcP->Pt());
1113
1114       if(mcP->Eta()>fEtamin && mcP->Eta()<fEtamax){
1115         h1_Pi0TruthPtPhi2piEta065       ->Fill(mcP->Pt());
1116         if(isPrimary==1)
1117           h1_PriPi0TruthPtPhi2piEta065  ->Fill(mcP->Pt());
1118         if(isK0sDecay)
1119           h1_K0Pi0TruthPtPhi2piEta065   ->Fill(mcP->Pt());
1120         if(isPrimary!=1 && isMaterialSec!=1)
1121           h1_PhysPi0TruthPtPhi2piEta065 ->Fill(mcP->Pt());
1122       }
1123
1124
1125       if(binp && bacc){// 2 Photons hit the EMCAL!
1126
1127         Int_t Nfoundphotons = 0;
1128         Int_t iFoundphotons[10] = {0,0,0,0,0,0,0,0,0,0};
1129         Int_t Nfoundelectrons = 0;
1130         Int_t iFoundelectrons[10] = {0,0,0,0,0,0,0,0,0,0};
1131         for (Int_t daughter_index=0; daughter_index<2; daughter_index++){//both truth photons. (also includes conversions..)
1132
1133           const AliMCParticle *dmc = static_cast<const AliMCParticle *>(mcEvent->GetTrack(daughter[daughter_index]));
1134
1135           h1_TruthPhotonsEmcal->Fill(dmc->Pt());
1136           h2_TruthPhotonsPhiEta->Fill(dmc->Phi(),dmc->Eta());
1137
1138           for(int i=0; i<nclusters; i++) {
1139             Bool_t matches_pion_photon = 0;
1140
1141             AliESDCaloCluster* esdCluster=NULL;
1142             AliAODCaloCluster* aodCluster=NULL;
1143             if (fEsdEv)       esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1144             else if (fAodEv)  aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1145
1146             Double_t clustMC_phi, clustMC_eta;
1147
1148             if(fEsdEv){
1149
1150               if(esdCluster->IsEMCAL()){
1151
1152                 Float_t pos[3] = {0,0,0};
1153                 esdCluster->GetPosition(pos);
1154                 TVector3 vpos(pos);
1155                 //h1_Phi->Fill(vpos.Phi());
1156                 clustMC_phi = vpos.Phi();
1157                 clustMC_eta = vpos.Eta();
1158
1159                 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
1160                                           (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
1161                 h1_dR_RealMC->Fill(dR);
1162                 matches_pion_photon = 0;
1163                 //if(dR<=0.04) matches_pion_photon = 1;
1164
1165
1166                 TArrayI *TruthLabelsA = esdCluster->GetLabelsArray();
1167                 if(TruthLabelsA){
1168                   Int_t trackindex = TruthLabelsA->At(0);
1169                   if( trackindex==daughter[daughter_index] ){
1170                     matches_pion_photon = 1;
1171                     iFoundphotons[Nfoundphotons] = i;
1172                     Nfoundphotons++;
1173                   }
1174                   else if( trackindex==eIndexofConvertedPhoton[daughter_index] ){
1175                     iFoundelectrons[Nfoundelectrons] = i;
1176                     Nfoundelectrons++;
1177                   }
1178                   AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex));
1179
1180                   if(matches_pion_photon){
1181
1182                     h1_PhotonsEmcal->Fill(esdCluster->E());
1183                     h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1184                     if(esdCluster->GetNCells()>=2)
1185                       h1_PhotonsNCellsCut->Fill(esdCluster->E());
1186                     if(esdCluster->GetNTracksMatched()==0)
1187                       h1_PhotonsTrackMatchCut->Fill(esdCluster->E());
1188                     if(esdCluster->GetNCells()>=2 && esdCluster->GetNTracksMatched()==0)
1189                       h1_PhotonsAllCut->Fill(esdCluster->E());
1190
1191                     if(esdCluster->GetNCells()>=2){
1192                       recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1193                       h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()));
1194                     }
1195                   }//if(matches_pion_photon)
1196
1197                 }//if Truthlabels exists
1198                 vpos.Delete();
1199               }//if(IsEMCAL())
1200
1201             }//if(fEsdEv)
1202             else if(fAodEv){
1203
1204               if(aodCluster->IsEMCAL()){
1205
1206                 Float_t pos[3] = {0,0,0};
1207                 aodCluster->GetPosition(pos);
1208                 TVector3 vpos(pos);
1209                 //h1_Phi->Fill(vpos.Phi());
1210                 clustMC_phi = vpos.Phi();
1211                 clustMC_eta = vpos.Eta();
1212
1213                 Double_t dR = TMath::Sqrt((eta_d[daughter_index]-clustMC_eta)*(eta_d[daughter_index]-clustMC_eta) +
1214                                           (phi_d[daughter_index]-clustMC_phi)*(phi_d[daughter_index]-clustMC_phi));
1215                 h1_dR_RealMC->Fill(dR);
1216                 matches_pion_photon = 0;
1217                 if(dR<=0.04) matches_pion_photon = 1;
1218
1219                 //TArrayI *TruthLabelsA = aodCluster->GetLabelsArray();
1220                 //if(TruthLabelsA){
1221                 //  Int_t trackindex = TruthLabelsA->At(0);
1222                 //  if( trackindex==daughter[daughter_index] )
1223                 //    matches_pion_photon = 1;
1224                 //  AliMCParticle *truthP = (AliMCParticle*)(mcEvent->GetTrack(trackindex));
1225
1226                 if(matches_pion_photon){
1227
1228                   h1_PhotonsEmcal->Fill(aodCluster->E());
1229                   h2_PhotonsPhiEtaIsEmcal->Fill(clustMC_phi,clustMC_eta);
1230                   if(aodCluster->GetNCells()>=2)
1231                     h1_PhotonsNCellsCut->Fill(aodCluster->E());
1232                   if(aodCluster->GetNTracksMatched()==0)
1233                     h1_PhotonsTrackMatchCut->Fill(aodCluster->E());
1234                   if(aodCluster->GetNCells()>=2 && aodCluster->GetNTracksMatched()==0)
1235                     h1_PhotonsAllCut->Fill(aodCluster->E());
1236
1237                   //if(aodCluster->GetNCells()>=2){
1238                   //recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1239                   //h2_gE_RecTruth->Fill(recalScale*esdCluster->E(), truthP->E()/(recalScale*esdCluster->E()));
1240                   //}
1241                 }//if(matches_pion_photon)
1242
1243                 //}//if Truthlabels exists
1244                 vpos.Delete();
1245               }//if(IsEMCAL())
1246
1247             }//if(fAodEv)
1248
1249           }//loop over nclusters.
1250
1251         }//both truth photons.
1252
1253         if(Nfoundphotons>1){
1254           AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1255           AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundphotons[1]); // pointer to EMCal cluster
1256
1257           if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1258
1259             recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1260             E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1261             fEsdEv->GetVertex()->GetXYZ(vertex);
1262             esdCluster1->GetMomentum(Photon1,vertex);
1263             Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1264             Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1265             Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1266
1267             recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1268             E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1269             fEsdEv->GetVertex()->GetXYZ(vertex);
1270             esdCluster2->GetMomentum(Photon2,vertex);
1271             Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1272             Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1273             Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1274
1275             Parent =  TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1276
1277             //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1278             Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1279                                                (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1280             if(isPrimary==1){
1281               //cout << "Primary production vertex: " << productionR << endl;
1282               h2_Mpt_Pri->Fill(Parent.M(),Parent.Pt());
1283             }
1284             else{
1285               //cout << "Secondary production vertex: " << productionR << endl;
1286               h2_Mpt_Sec ->Fill(Parent.M(),Parent.Pt());
1287               h3_MptR_Sec->Fill(Parent.M(),Parent.Pt(),productionR);
1288               if(isK0sDecay)
1289                 h3_MptR_K0s->Fill(Parent.M(),Parent.Pt(),productionR);
1290               if(isMaterialSec)
1291                 h3_MptR_Mat->Fill(Parent.M(),Parent.Pt(),productionR);
1292             }
1293           }//both good clusters
1294         }//found 2 photons.
1295         else if(Nfoundphotons==1){
1296           int mergedPion = 0;
1297           AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1298
1299           TArrayI *TruthLabelsA = esdCluster1->GetLabelsArray();
1300           if(TruthLabelsA){
1301             if(TruthLabelsA->GetSize()>1){
1302               Int_t trackindex[2];
1303               trackindex[0] = TruthLabelsA->At(0);
1304               trackindex[1] = TruthLabelsA->At(1);
1305               if( (trackindex[0]==daughter[0] && trackindex[1]==daughter[1]) ||
1306                   (trackindex[0]==daughter[1] && trackindex[1]==daughter[0]) ){
1307                 mergedPion = 1;
1308               }
1309               if(mergedPion==1){
1310                 recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1311                 E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1312                 Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1313                                                    (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1314                 h2_PtR_MatM->Fill(E1,productionR);
1315               }//if merged pion.
1316             }//truthlabel.size > 1
1317           }//if truthlabels
1318         }// Nfoundphotons==1
1319
1320         if(Nfoundphotons==1 && Nfoundelectrons==1){
1321           AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundphotons[0]); // pointer to EMCal cluster
1322           AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster
1323
1324           if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1325
1326             recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1327             E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1328             fEsdEv->GetVertex()->GetXYZ(vertex);
1329             esdCluster1->GetMomentum(Photon1,vertex);
1330             Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1331             Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1332             Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1333
1334             recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1335             E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1336             fEsdEv->GetVertex()->GetXYZ(vertex);
1337             esdCluster2->GetMomentum(Photon2,vertex);
1338             Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1339             Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1340             Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1341
1342             Parent =  TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1343
1344             //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1345             Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1346                                                (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1347             if(isPrimary==1){
1348               //cout << "Primary production vertex: " << productionR << endl;
1349               h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt());
1350             }
1351             else{
1352               //cout << "Secondary production vertex: " << productionR << endl;
1353               h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt());
1354               h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1355               if(isK0sDecay)
1356                 h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1357               if(isMaterialSec)
1358                 h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1359             }
1360           }//both good clusters
1361         }// Nfoundphotons==1 && Nfoundelectrons==1
1362         else if(Nfoundelectrons==2){
1363           AliESDCaloCluster* esdCluster1 = fEsdEv->GetCaloCluster(iFoundelectrons[0]); // pointer to EMCal cluster
1364           AliESDCaloCluster* esdCluster2 = fEsdEv->GetCaloCluster(iFoundelectrons[1]); // pointer to EMCal cluster
1365
1366           if( isGoodEsdCluster(esdCluster1) && isGoodEsdCluster(esdCluster2) ){
1367
1368             recalScale = PrivateEnergyRecal(esdCluster1->E(), fRecalibrator);
1369             E1 = esdCluster1->E()*recalScale;// TOTAL HACK - JJ
1370             fEsdEv->GetVertex()->GetXYZ(vertex);
1371             esdCluster1->GetMomentum(Photon1,vertex);
1372             Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1373             Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1374             Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1375
1376             recalScale = PrivateEnergyRecal(esdCluster2->E(), fRecalibrator);
1377             E2 = esdCluster2->E()*recalScale;// TOTAL HACK - JJ
1378             fEsdEv->GetVertex()->GetXYZ(vertex);
1379             esdCluster2->GetMomentum(Photon2,vertex);
1380             Photon2.SetPx(Photon2.Px()*recalScale);// TOTAL HACK - JJ
1381             Photon2.SetPy(Photon2.Py()*recalScale);// TOTAL HACK - JJ
1382             Photon2.SetPz(Photon2.Pz()*recalScale);// TOTAL HACK - JJ
1383
1384             Parent =  TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) + TLorentzVector(Photon2.Px(),Photon2.Py(),Photon2.Pz(),E2);
1385
1386             //double productionR = TMath::Sqrt( mcP->Xv()*mcP->Xv() + mcP->Yv()*mcP->Yv() );
1387             Double_t productionR = TMath::Sqrt((mcP->Xv()-evtVtx->GetX())*(mcP->Xv()-evtVtx->GetX()) +
1388                                                (mcP->Yv()-evtVtx->GetY())*(mcP->Yv()-evtVtx->GetY()));
1389             if(isPrimary==1){
1390               //cout << "Primary production vertex: " << productionR << endl;
1391               h2_Mpt_Pri_conv->Fill(Parent.M(),Parent.Pt());
1392             }
1393             else{
1394               //cout << "Secondary production vertex: " << productionR << endl;
1395               h2_Mpt_Sec_conv ->Fill(Parent.M(),Parent.Pt());
1396               h3_MptR_Sec_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1397               if(isK0sDecay)
1398                 h3_MptR_K0s_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1399               if(isMaterialSec)
1400                 h3_MptR_Mat_conv->Fill(Parent.M(),Parent.Pt(),productionR);
1401             }
1402           }//both good clusters
1403         }// Nfoundelectrons==2
1404
1405         h1_Pi0TruthPtEmcal    ->Fill(mcP->Pt());
1406         if(isK0sDecay)
1407           h1_K0Pi0TruthPtEmcal    ->Fill(mcP->Pt());
1408         h2_Pi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1409
1410         if(isPrimary==1){
1411           h1_PriPi0TruthPtEmcal    ->Fill(mcP->Pt());
1412           h2_PriPi0TruthPhiEtaEmcal->Fill(mcP->Phi(),mcP->Eta());
1413         }
1414         if(isPrimary!=1 && isMaterialSec!=1)   h1_PhysPi0TruthPtEmcal->Fill(mcP->Pt());
1415
1416
1417       }// 2 Photons hit the EMCAL!
1418
1419     }//for(nTracksMC)
1420
1421   }//if(isMC)
1422
1423   //######################### ~~~~~~~~~~~~ ##################################
1424   //######################### DONE WITH MC ##################################
1425   //######################### ~~~~~~~~~~~~ ##################################
1426
1427
1428   for(int i=0; i<nclusters; i++) {
1429
1430     AliESDCaloCluster* esdCluster=NULL;
1431     AliAODCaloCluster* aodCluster=NULL;
1432     if (fEsdEv)       esdCluster = fEsdEv->GetCaloCluster(i); // pointer to EMCal cluster
1433     else if (fAodEv)  aodCluster = fAodEv->GetCaloCluster(i); // pointer to EMCal cluster
1434     if(!esdCluster && !aodCluster) {
1435       AliError(Form("ERROR: Could not retrieve any (ESD or AOD) Cluster %d",i));
1436       continue;
1437     }
1438
1439     if(fEsdEv){
1440
1441       recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1442
1443       //uncomment this to do the track matching (1 of 3 lines, esd part)!!
1444       //Bool_t MatchesToTrack = 0;
1445       if(esdCluster->IsEMCAL()){
1446
1447         Float_t pos[3] = {0,0,0};
1448         Short_t maxCellID = -1;
1449         Float_t celleta, cellphi;
1450         esdCluster->GetPosition(pos);
1451         TVector3 clusterPosition(pos);
1452         h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1453         GetMaxCellEnergy(esdCluster, maxCellID);
1454         AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1455         fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1456         h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1457
1458         // _______________Track loop for reconstructed event_____________
1459         for(Int_t itrk = 0; itrk < fEsdEv->GetNumberOfTracks(); itrk++) {
1460           AliESDtrack* esdTrack = fEsdEv->GetTrack(itrk); // pointer to reconstructed to track
1461           if(!esdTrack) {
1462             AliError(Form("ERROR: Could not retrieve any (ESD) track %d",itrk));
1463             continue;
1464           }
1465
1466           Double_t posTrk[3] = {0,0,0};
1467           esdTrack->GetXYZ(posTrk);
1468           TVector3 vposTrk(posTrk);
1469
1470           Double_t fMass          = 0.139;
1471           Double_t fStepSurface   = 20.;
1472           Float_t etaproj, phiproj, pttrackproj;
1473
1474           AliExternalTrackParam *trackParam =  const_cast<AliExternalTrackParam*>(esdTrack->GetInnerParam());
1475           if(!trackParam) continue;
1476           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1477
1478           double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1479                                     (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1480
1481           h1_dR_ClustTrk->Fill(dR_clusttrk);
1482
1483           //uncomment this to do the track matching (2 of 3 lines, esd part)!!
1484           //if(dR_clusttrk<fdRmin_ClustTrack)
1485           //MatchesToTrack = 1;
1486
1487         }//_____________________________nTracks__________________________
1488
1489         h2_cells_M02  ->Fill(esdCluster->GetNCells(),esdCluster->GetM02());
1490         h2_Ellipse    ->Fill(esdCluster->GetM20(),esdCluster->GetM02());
1491         h1_Chi2       ->Fill(esdCluster->Chi2());//always -1.
1492         h1_nTrkMatch  ->Fill(esdCluster->GetNTracksMatched());
1493         h1_ClusterDisp->Fill(esdCluster->GetDispersion());
1494         h2_E_time     ->Fill(esdCluster->E(),esdCluster->GetTOF());
1495
1496         TArrayI *TrackLabels = esdCluster->GetTracksMatched();
1497         if(TrackLabels){
1498           if(TrackLabels->GetSize()>0){
1499             Int_t trackindex = TrackLabels->At(0);
1500             AliESDtrack* matchingT = fEsdEv->GetTrack(trackindex); // pointer to reconstructed to track
1501
1502             recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1503             h2_eop_E ->Fill(esdCluster->E()*recalScale, esdCluster->E()*recalScale/matchingT->P());
1504             h2_eop_pT->Fill(matchingT->Pt(),            esdCluster->E()*recalScale/matchingT->P());
1505           }
1506         }
1507
1508         //uncomment this to do the track matching (2 of 3 lines, esd part)!!
1509         //if(isGoodEsdCluster(esdCluster) && !MatchesToTrack){
1510         if(isGoodEsdCluster(esdCluster)){
1511           recalScale = PrivateEnergyRecal(esdCluster->E(), fRecalibrator);
1512           E1 = esdCluster->E()*recalScale;// TOTAL HACK - JJ
1513           fEsdEv->GetVertex()->GetXYZ(vertex);
1514           esdCluster->GetMomentum(Photon1,vertex);
1515           Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1516           Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1517           Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1518           Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1519           h1_E->Fill(E1);
1520           h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1521           h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1522         }
1523         clusterPosition.Delete();
1524       }//if(esdCluster->isEMCAL())
1525     }//if(fEsdEv)
1526     else if(fAodEv){
1527
1528       recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1529
1530       //uncomment this to do the track matching (1 of 3 lines, aod part)!!
1531       //Bool_t MatchesToTrack = 0;
1532       if(aodCluster->IsEMCAL()){
1533
1534         Float_t pos[3] = {0,0,0};
1535         Short_t maxCellID = -1;
1536         Float_t celleta, cellphi;
1537         aodCluster->GetPosition(pos);
1538         TVector3 clusterPosition(pos);
1539         h2_PhiEtaCluster->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1540         GetMaxCellEnergy(aodCluster, maxCellID);
1541         AliEMCALGeometry *fGeom = AliEMCALGeometry::GetInstance();
1542         fGeom->EtaPhiFromIndex(maxCellID,celleta,cellphi);
1543         h2_PhiEtaMaxCell->Fill(cellphi,celleta);
1544
1545         // _______________Track loop for reconstructed event_____________
1546         for(Int_t itrk = 0; itrk < fAodEv->GetNumberOfTracks(); itrk++) {
1547           AliAODTrack* aodTrack = fAodEv->GetTrack(itrk); // pointer to reconstructed to track
1548           if(!aodTrack) {
1549             AliError(Form("ERROR: Could not retrieve any (AOD) track %d",itrk));
1550             continue;
1551           }
1552
1553           Double_t posTrk[3] = {0,0,0};
1554           Double_t momTrk[3] = {0,0,0};
1555           aodTrack->GetXYZ(posTrk);
1556           aodTrack->GetPxPyPz(momTrk);
1557           //TVector3 vposTrk(posTrk);
1558
1559           //####################################################################################################
1560           //
1561           // commented all this stuff just to satisfy aliroot warnings.
1562           // but I may need it again if I want to do the track matching for aods.
1563           /*
1564           Double_t fMass          = 0.139;
1565           Double_t fStepSurface   = 20.;
1566           Float_t etaproj=0.0;
1567           Float_t phiproj=0.0;
1568           Float_t pttrackproj=0.0;
1569
1570           Double_t cv[21] = {0.0};
1571           aodTrack->GetCovarianceXYZPxPyPz(cv);
1572           AliExternalTrackParam *trackParam = new AliExternalTrackParam(posTrk,momTrk,cv,aodTrack->Charge());
1573           //AliExternalTrackParam emcalParam(*trackParam);
1574           //AliExternalTrackParam *trackParam =  const_cast<AliExternalTrackParam*>(aodTrack->GetInnerParam());
1575           if(!trackParam) continue;
1576           ////AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(trackParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1577           //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(&emcalParam, 440., fMass, fStepSurface, etaproj, phiproj, pttrackproj);
1578           delete trackParam;
1579
1580           //Constantin's implementation... gives funny result.
1581           //AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(aodTrack,440.0);
1582           //phiproj = aodTrack->GetTrackPhiOnEMCal();
1583           //etaproj = aodTrack->GetTrackPhiOnEMCal();
1584
1585           double dR_clusttrk = sqrt((phiproj-clusterPosition.Phi())*(phiproj-clusterPosition.Phi()) +
1586                                     (etaproj-clusterPosition.Eta())*(etaproj-clusterPosition.Eta()) );
1587
1588           h1_dR_ClustTrk->Fill(dR_clusttrk);
1589           */
1590           //####################################################################################################
1591
1592           //uncomment this to do the track matching (2 of 3 lines, aod part)!!
1593           //if(dR_clusttrk<fdRmin_ClustTrack)
1594           //MatchesToTrack = 1;
1595
1596
1597         }//_____________________________nTracks__________________________
1598
1599         h2_cells_M02  ->Fill(aodCluster->GetNCells(),aodCluster->GetM02());
1600         h2_Ellipse    ->Fill(aodCluster->GetM20(),aodCluster->GetM02());
1601         h1_Chi2       ->Fill(aodCluster->Chi2());//always -1.
1602         h1_nTrkMatch  ->Fill(aodCluster->GetNTracksMatched());
1603         h1_ClusterDisp->Fill(aodCluster->GetDispersion());
1604         h2_E_time     ->Fill(aodCluster->E(),aodCluster->GetTOF());
1605
1606         // #################################################
1607         // track matching eop histograms are handled here...
1608         // #################################################
1609
1610         //uncomment this to do the track matching (3 of 3 lines, aod part)!!
1611         //if(isGoodAodCluster(aodCluster) && !MatchesToTrack){
1612         if(isGoodAodCluster(aodCluster)){
1613           recalScale = PrivateEnergyRecal(aodCluster->E(), fRecalibrator);
1614           E1 = aodCluster->E()*recalScale;// TOTAL HACK - JJ
1615           fAodEv->GetVertex(0)->GetXYZ(vertex);
1616           aodCluster->GetMomentum(Photon1,vertex);
1617           Photon1.SetPx(Photon1.Px()*recalScale);// TOTAL HACK - JJ
1618           Photon1.SetPy(Photon1.Py()*recalScale);// TOTAL HACK - JJ
1619           Photon1.SetPz(Photon1.Pz()*recalScale);// TOTAL HACK - JJ
1620           Photons[0][izvtx][imult].push_back( TLorentzVector(Photon1.Px(),Photon1.Py(),Photon1.Pz(),E1) );
1621           h1_E->Fill(E1);
1622           h2_PhiEtaClusterCut->Fill(clusterPosition.Phi(),clusterPosition.Eta());
1623           h2_PhiEtaMaxCellCut->Fill(cellphi,celleta);
1624         }
1625         clusterPosition.Delete();
1626       }//if(aodCluster->IsEMCAL())
1627     }//if(fAodEv)
1628
1629   }//loop over nclusters.
1630
1631   //Make same event pions...
1632   for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1633     for(unsigned int j=i+1; j<Photons[0][izvtx][imult].size(); j++){
1634       Parent = Photons[0][izvtx][imult][i] + Photons[0][izvtx][imult][j];
1635       Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1636       Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[0][izvtx][imult][j]);
1637       Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[0][izvtx][imult][j].Pt())/
1638                               (Photons[0][izvtx][imult][i].Pt()+Photons[0][izvtx][imult][j].Pt());
1639       Int_t asymCut = 0;
1640       if     (pairasym<0.1)  asymCut = 1;
1641       else if(pairasym<0.7)  asymCut = 2;
1642       else                   asymCut = 3;
1643
1644       h1_M        ->Fill(Parent.M());
1645       h3_MptAsymm ->Fill(Parent.M(),Parent.Pt(),asymCut);
1646       h2_dphi_deta->Fill(deltaphi,deltaeta);
1647     }
1648   }
1649
1650   //Make mixed event...
1651   for(unsigned int i=0; i<Photons[0][izvtx][imult].size(); i++){
1652     for(unsigned int ipool=1; ipool<poolDepth; ipool++){
1653       for(unsigned int j=0; j<Photons[ipool][izvtx][imult].size(); j++){
1654         iskip = randy.Integer(Photons[0][izvtx][imult].size());
1655         if(j==iskip) continue;
1656         Parent = Photons[0][izvtx][imult][i]+Photons[ipool][izvtx][imult][j];
1657         Double_t deltaphi = getDeltaPhi(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1658         Double_t deltaeta = getDeltaEta(Photons[0][izvtx][imult][i],Photons[ipool][izvtx][imult][j]);
1659         Double_t pairasym = fabs(Photons[0][izvtx][imult][i].Pt()-Photons[ipool][izvtx][imult][j].Pt())/
1660                                 (Photons[0][izvtx][imult][i].Pt()+Photons[ipool][izvtx][imult][j].Pt());
1661         Int_t asymCut = 0;
1662         if     (pairasym<0.1)  asymCut = 1;
1663         else if(pairasym<0.7)  asymCut = 2;
1664         else                   asymCut = 3;
1665
1666         h1_M_mix        ->Fill(Parent.M());
1667         h3_MptAsymm_mix ->Fill(Parent.M(),Parent.Pt(),asymCut);
1668         h2_dphi_deta_mix->Fill(deltaphi,deltaeta);
1669       }
1670     }
1671   }
1672
1673   for(int ipool=poolDepth-1; ipool>0; ipool--){
1674     Photons[ipool][izvtx][imult].clear();
1675     for(unsigned int i=0; i<Photons[ipool-1][izvtx][imult].size(); i++)
1676       Photons[ipool][izvtx][imult].push_back(Photons[ipool-1][izvtx][imult][i]);
1677   }
1678   Photons[0][izvtx][imult].clear();
1679
1680
1681
1682   // NEW HISTO should be filled before this point, as PostData puts the
1683   // information for this iteration of the UserExec in the container
1684   PostData(1, fOutput);
1685   }
1686
1687 //________________________________________________________________________
1688 void AliAnalysisTaskEMCALMesonGGSDM::Terminate(Option_t *) //specify what you want to have done
1689 {
1690   // Called once at the end of the query.
1691
1692 }
1693
1694 //________________________________________________________________________
1695 Int_t AliAnalysisTaskEMCALMesonGGSDM::GetZvtxBin(Double_t vertZ)
1696 {
1697
1698   int izvtx = -1;
1699
1700   if     (vertZ<-35)
1701     izvtx=0;
1702   else if(vertZ<-30)
1703     izvtx=1;
1704   else if(vertZ<-25)
1705     izvtx=2;
1706   else if(vertZ<-20)
1707     izvtx=3;
1708   else if(vertZ<-15)
1709     izvtx=4;
1710   else if(vertZ<-10)
1711     izvtx=5;
1712   else if(vertZ< -5)
1713     izvtx=6;
1714   else if(vertZ<  0)
1715     izvtx=7;
1716   else if(vertZ<  5)
1717     izvtx=8;
1718   else if(vertZ< 10)
1719     izvtx=9;
1720   else if(vertZ< 15)
1721     izvtx=10;
1722   else if(vertZ< 20)
1723     izvtx=11;
1724   else if(vertZ< 25)
1725     izvtx=12;
1726   else if(vertZ< 30)
1727     izvtx=13;
1728   else if(vertZ< 35)
1729     izvtx=14;
1730   else
1731     izvtx=15;
1732
1733   return izvtx;
1734 }
1735
1736 //________________________________________________________________________
1737 Int_t AliAnalysisTaskEMCALMesonGGSDM::GetMultBin(Int_t mult){
1738
1739   int imult = -1;
1740
1741   if     (mult<2)
1742     imult=0;
1743   else if(mult<25)
1744     imult=mult-2;
1745   else
1746     imult=24;
1747
1748   return imult;
1749 }
1750
1751 //________________________________________________________________________
1752 Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodEsdCluster(AliESDCaloCluster* esdclust){
1753
1754   int pass = 1;
1755   int nMinCells  = 2;
1756   double MinE    = 0.4;
1757   //double MinErat = 0;
1758   //double MinEcc  = 0;
1759
1760   if (!esdclust)
1761     pass = 0;
1762   if (!esdclust->IsEMCAL())
1763     pass = 0;
1764   if (esdclust->E()<MinE)
1765     pass = 0;
1766   if (esdclust->GetNCells()<nMinCells)
1767     pass = 0;
1768   //if (GetMaxCellEnergy(esdclust)/esdclust->E()<MinErat)
1769   //pass = 0;
1770   //if (esdclust->Chi2()<MinEcc) // eccentricity cut
1771   //pass = 0;//this is always -1.
1772
1773   //if(esdclust->GetM02()<0.1)
1774   //  pass = 0;
1775   //if(esdclust->GetM02()>0.5)
1776   //  pass = 0;
1777
1778   Float_t pos[3] = {0,0,0};
1779   esdclust->GetPosition(pos);
1780   TVector3 clusterPosition(pos);
1781   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1782      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
1783     pass = 0;
1784   clusterPosition.Delete();
1785
1786   //DOING THIS BY HAND NOW...
1787   //if(!esdclust->GetNTracksMatched()==0)
1788   //pass = 0;
1789
1790   return pass;
1791 }
1792
1793 //________________________________________________________________________
1794 Int_t AliAnalysisTaskEMCALMesonGGSDM::isGoodAodCluster(AliAODCaloCluster* aodclust){
1795
1796   int pass = 1;
1797   int nMinCells  = 2;
1798   double MinE    = 0.4;
1799   //double MinErat = 0;
1800   //double MinEcc  = 0;
1801
1802   if (!aodclust)
1803     pass = 0;
1804   if (!aodclust->IsEMCAL())
1805     pass = 0;
1806   if (aodclust->E()<MinE)
1807     pass = 0;
1808   if (aodclust->GetNCells()<nMinCells)
1809     pass = 0;
1810   //if (GetMaxCellEnergy(aodclust)/aodclust->E()<MinErat)
1811   //pass = 0;
1812   //if (aodclust->Chi2()<MinEcc) // eccentricity cut
1813   //pass = 0;//this is always -1.
1814
1815   //if(aodclust->GetM02()<0.1)
1816   //pass = 0;
1817   //if(aodclust->GetM02()>0.5)
1818   //pass = 0;
1819
1820   Float_t pos[3] = {0,0,0};
1821   aodclust->GetPosition(pos);
1822   TVector3 clusterPosition(pos);
1823   if(clusterPosition.Eta()<fEtamin || clusterPosition.Eta()>fEtamax ||
1824      clusterPosition.Phi()<fPhimin || clusterPosition.Phi()>fPhimax  )
1825     pass = 0;
1826   clusterPosition.Delete();
1827
1828   //DOING THIS BY HAND NOW...
1829   //if(!aodclust->GetNTracksMatched()==0)
1830   //pass = 0;
1831
1832   return pass;
1833 }
1834
1835 //________________________________________________________________________
1836 Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaPhi(TLorentzVector p1, TLorentzVector p2){
1837
1838   double dphi = p1.Phi() - p2.Phi();
1839
1840   if(dphi<0.5*TMath::Pi())
1841     dphi = dphi + 2.0*TMath::Pi();
1842
1843   if(dphi>1.5*TMath::Pi())
1844     dphi = dphi - 2.0*TMath::Pi();
1845
1846   return dphi;
1847 }
1848
1849 //________________________________________________________________________
1850 Double_t AliAnalysisTaskEMCALMesonGGSDM::getDeltaEta(TLorentzVector p1, TLorentzVector p2){
1851
1852   double deta = p1.PseudoRapidity() - p2.PseudoRapidity();
1853
1854   return deta;
1855 }
1856
1857
1858 //________________________________________________________________________
1859 Double_t AliAnalysisTaskEMCALMesonGGSDM::PrivateEnergyRecal(Double_t energy, Int_t iCalib){
1860
1861   double recalibfactor = 0.0;
1862
1863   if(iCalib==0){// no recalibration!
1864     recalibfactor = 1.0;
1865   }
1866   else if(iCalib==1){// just a scale factor:
1867     recalibfactor = 0.984;
1868   }
1869   else if(iCalib==2){// Symmetric Decay Fit - corrects data to uncorrected MC.
1870     Double_t p[3] = {0.96968, -2.68720, -0.831607};
1871     recalibfactor = p[0] + exp(p[1] + p[2]*energy*2.0);
1872   }
1873   else if(iCalib==3){// Jason's fit to the LHC12f1a MC single photons - 04 Aug 2013 (call it kPi0MCv4??)
1874     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};
1875     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1876   }
1877   else if(iCalib==4){// Jason's fit to the test beam data - 04 Aug 2013(call it kBTCv3??)
1878     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};
1879     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1880   }
1881   else if(iCalib==5){// Based on kSDM/kTBCv3 (call it kPi0MCv4??)
1882     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};
1883     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) );
1884   }
1885   else if(iCalib==6){// kBeamTestCorrectedv2 - in AliROOT!
1886     Double_t p[7] = {9.83504e-01, 2.10106e-01, 8.97274e-01, 8.29064e-02, 1.52299e+02, 3.15028e+01, 0.968};
1887     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1888   }
1889   else if(iCalib==7){// kPi0MCv3 - in AliROOT!
1890     Double_t p[7] = {9.81039e-01, 1.13508e-01, 1.00173e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.0};
1891     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1892   }
1893   else if(iCalib==8){// Jason's fit to the noNL MC/data- based on kSDM and kPi0MCv5 - 28 Oct 2013 (call it... ??)
1894     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};
1895     //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
1896     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));
1897   }
1898   else if(iCalib==9){// Jason's fit to the LHC12f1a/b MC single photons (above 400MeV), including conversions - 28 Oct 2013 (call it kPi0MCv5??)
1899     Double_t p[7] = {1.0, 6.64778e-02, 1.57000e+00, 9.67998e-02, 2.19381e+02, 6.31604e+01, 1.01286};
1900     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1901   }
1902   else if(iCalib==10){// Jason played with test beam data
1903     Double_t p[7] = {1.0, 0.237767, 0.651203, 0.183741, 155.427, 17.0335, 0.987054};
1904     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1905   }
1906   else if(iCalib==11){// Jason played with test beam MC
1907     Double_t p[7] = {1.0, 0.0797873, 1.68322, 0.0806098, 244.586, 116.938, 1.00437};
1908     recalibfactor = ((p[6])/(p[0]*(1./(1.+p[1]*exp(-energy/p[2]))*1./(1.+p[3]*exp((energy-p[4])/p[5])))));
1909   }
1910
1911   return recalibfactor;
1912 }
1913
1914
1915 //________________________________________________________________________
1916 Double_t AliAnalysisTaskEMCALMesonGGSDM::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
1917 {
1918   // Get maximum energy of attached cell.
1919
1920   id = -1;
1921   AliVCaloCells *fVCells=NULL;
1922   if(fEsdEv)      fVCells = fEsdEv->GetEMCALCells();
1923   else if(fAodEv) fVCells = fAodEv->GetEMCALCells();
1924   if(!fVCells)
1925     return 0;
1926
1927   Double_t maxe = 0;
1928   Int_t ncells = cluster->GetNCells();
1929   for (Int_t i=0; i<ncells; i++) {
1930     Double_t e = fVCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
1931     if (e>maxe) {
1932       maxe = e;
1933       id   = cluster->GetCellAbsId(i);
1934     }
1935   }
1936   return maxe;
1937 }
1938
1939
1940 //________________________________________________________________________
1941 Int_t AliAnalysisTaskEMCALMesonGGSDM::IsPhysPrimJ(AliMCEvent *mcEvent, Int_t iTrack){
1942
1943   AliMCParticle *mcP  = static_cast<AliMCParticle*>(mcEvent->GetTrack(iTrack));
1944
1945   Int_t nPTracks= mcEvent->GetNumberOfPrimaries();
1946
1947   Int_t isPhysPrimary   = 1;
1948   Int_t ismHF           = 0;
1949   Int_t ismLongLivedOrK = 0;
1950
1951   if(mcP->GetMother()<0)//if it has no mother...
1952     return isPhysPrimary;
1953
1954   Int_t imTrack = mcP->GetMother();
1955   AliMCParticle *mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1956
1957   if( TMath::Abs(mcPm->PdgCode())<10 )//if mother is a single quark...
1958     return isPhysPrimary;
1959
1960
1961   //############################################
1962   //get the PDG digits....
1963   int num = mcPm->PdgCode();
1964   int RevDigits[10] = {0};
1965   int nDigits = 0;
1966   while (num >= 1){
1967     RevDigits[nDigits++] = num%10;
1968     num = num / 10;
1969   }
1970   //##############################################
1971
1972
1973   if(RevDigits[3]>3)//Baryons
1974     ismHF = 1;
1975   else if(RevDigits[2]>3)//Mesons
1976     ismHF = 1;
1977
1978   ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
1979
1980   if(!ismHF && ismLongLivedOrK)
1981     isPhysPrimary = 0;
1982   else{ // check grandmother, greatgrandmothers, etc...
1983     while(imTrack >= nPTracks){
1984
1985       if(mcPm->GetMother()<0)//if it has no mother...
1986         break;
1987
1988       if( TMath::Abs(mcPm->PdgCode()<10) )//if mother is a single quark...
1989         return isPhysPrimary;
1990
1991       imTrack = mcPm->GetMother();
1992       mcPm = static_cast<AliMCParticle*>(mcEvent->GetTrack(imTrack));
1993
1994       //############################################
1995       //get the PDG digits....
1996       num = mcPm->PdgCode();
1997       for(int i=0; i<10; i++)  RevDigits[i] = 0;
1998       nDigits = 0;
1999       while (num >= 1){
2000         RevDigits[nDigits++] = num%10;
2001         num = num / 10;
2002       }
2003       //##############################################
2004       if(RevDigits[3]>3)//Baryons
2005         ismHF = 1;
2006       else if(RevDigits[2]>3)//Mesons
2007         ismHF = 1;
2008
2009       ismLongLivedOrK = IsLongLivedOrK(mcPm->PdgCode());
2010
2011       if(!ismHF && ismLongLivedOrK)
2012         isPhysPrimary = 0;
2013
2014     }//while( >=nPTracks)
2015   }
2016
2017   return isPhysPrimary;
2018 }
2019
2020
2021 //________________________________________________________________________
2022 Int_t AliAnalysisTaskEMCALMesonGGSDM::IsLongLivedOrK(Int_t MyPDGcode){
2023
2024   Int_t MyFlag = 0;
2025
2026   if(
2027      (TMath::Abs(MyPDGcode) == 22  ) ||        // Photon
2028      (TMath::Abs(MyPDGcode) == 11  ) ||        // Electron
2029      (TMath::Abs(MyPDGcode) == 13  ) ||        // Muon(-)
2030      (TMath::Abs(MyPDGcode) == 211 ) ||        // Pion
2031      (TMath::Abs(MyPDGcode) == 321 ) ||        // Kaon
2032      (TMath::Abs(MyPDGcode) == 310 ) ||        // K0s
2033      (TMath::Abs(MyPDGcode) == 130 ) ||        // K0l
2034      (TMath::Abs(MyPDGcode) == 2212) ||        // Proton
2035      (TMath::Abs(MyPDGcode) == 2112) ||        // Neutron
2036      (TMath::Abs(MyPDGcode) == 3122) ||        // Lambda_0
2037      (TMath::Abs(MyPDGcode) == 3112) ||        // Sigma Minus
2038      (TMath::Abs(MyPDGcode) == 3222) ||        // Sigma Plus
2039      (TMath::Abs(MyPDGcode) == 3312) ||        // Xsi Minus
2040      (TMath::Abs(MyPDGcode) == 3322) ||        // Xsi
2041      (TMath::Abs(MyPDGcode) == 3334) ||        // Omega
2042      (TMath::Abs(MyPDGcode) == 12  ) ||        // Electron Neutrino
2043      (TMath::Abs(MyPDGcode) == 14  ) ||        // Muon Neutrino
2044      (TMath::Abs(MyPDGcode) == 16  )   )       // Tau Neutrino
2045     MyFlag = 1;
2046
2047   return MyFlag;
2048 }
2049
2050
2051 //________________________________________________________________________
2052 Int_t AliAnalysisTaskEMCALMesonGGSDM::IsMyMCHeaderType(Int_t iTrack, char *MyType, AliMCEvent *mcEvent) const
2053 {
2054
2055   Int_t isMyType = 0;
2056
2057   AliGenCocktailEventHeader *cocktail = dynamic_cast<AliGenCocktailEventHeader *>(mcEvent->GenEventHeader());
2058   if(!cocktail)
2059     return 0;
2060
2061   TList *genHeaders = cocktail->GetHeaders();
2062
2063   Int_t nGenerators = genHeaders->GetEntries();
2064   Int_t indexMyType = -1;
2065   Int_t startParticle=0;
2066
2067   for(Int_t igen = 0; igen < nGenerators; igen++){
2068     AliGenEventHeader* eventHeader2 = (AliGenEventHeader*)genHeaders->At(igen) ;
2069     TString name = eventHeader2->GetName();
2070     startParticle += eventHeader2->NProduced();
2071     //cout << name << endl;
2072     if (name.Contains(MyType,TString::kIgnoreCase)){
2073       indexMyType = igen;
2074       startParticle -= eventHeader2->NProduced();
2075       break;
2076     }
2077   }
2078
2079   AliGenEventHeader *addedPi0Header = (AliGenEventHeader*)genHeaders->At(indexMyType);
2080   Int_t ipi0min = startParticle;
2081   Int_t ipi0max = ipi0min+addedPi0Header->NProduced()-1;
2082   if(iTrack >= ipi0min && iTrack <= ipi0max)
2083     isMyType = 1;
2084
2085   return isMyType;
2086 }
2087
2088