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