]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/GammaConv/AliAnalysisTaskGammaConvCalo.cxx
85fe4c8cef85784c74d48e681a696ff77d802726
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskGammaConvCalo.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                                                                                                *
4  * Author: Baldo Sahlmueller, Friederike Bock                                             *
5  * Version 1.0                                                                                                                    *
6  *                                                                                                                                                *
7  *                                                                        *
8  * Permission to use, copy, modify and distribute this software and its   *
9  * documentation strictly for non-commercial purposes is hereby granted   *
10  * without fee, provided that the above copyright notice appears in all   *
11  * copies and that both the copyright notice and this permission notice   *
12  * appear in the supporting documentation. The authors make no claims     *
13  * about the suitability of this software for any purpose. It is                  *
14  * provided "as is" without express or implied warranty.                          *
15  **************************************************************************/
16
17 //////////////////////////////////////////////////////////////////
18 //----------------------------------------------------------------
19 // Class used to do analysis on conversion photons + calo photons
20 //----------------------------------------------------------------
21 //////////////////////////////////////////////////////////////////
22 #include "TChain.h"
23 #include "TTree.h"
24 #include "TBranch.h"
25 #include "TFile.h"
26 #include "TH1F.h"
27 #include "TH2F.h"
28 #include "TH3F.h"
29 #include "THnSparse.h"
30 #include "TCanvas.h"
31 #include "TNtuple.h"
32 #include "AliAnalysisTask.h"
33 #include "AliAnalysisManager.h"
34 #include "AliESDEvent.h"
35 #include "AliESDInputHandler.h"
36 #include "AliMCEventHandler.h"
37 #include "AliMCEvent.h"
38 #include "AliMCParticle.h"
39 #include "AliCentrality.h"
40 #include "AliESDVZERO.h"
41 #include "AliESDpid.h"
42 #include "AliAnalysisTaskGammaConvCalo.h"
43 #include "AliVParticle.h"
44 #include "AliESDtrack.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliKFVertex.h"
47 #include "AliV0ReaderV1.h"
48 #include "AliGenCocktailEventHeader.h"
49 #include "AliConversionAODBGHandlerRP.h"
50 #include "AliAODMCParticle.h"
51 #include "AliAODMCHeader.h"
52 #include "AliEventplane.h"
53 #include "AliAnalysisTaskEMCALClusterizeFast.h"
54 #include "AliAODEvent.h"
55 #include "AliAODInputHandler.h"
56 #include "AliESDEvent.h"
57 #include "AliESDInputHandler.h"
58 #include "AliInputEventHandler.h"
59
60 ClassImp(AliAnalysisTaskGammaConvCalo)
61
62 //________________________________________________________________________
63 AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(): AliAnalysisTaskSE(),
64         fV0Reader(NULL),
65         fBGHandler(NULL),
66         fBGHandlerRP(NULL),
67         fBGClusHandler(NULL),
68         fBGClusHandlerRP(NULL),
69         fInputEvent(NULL),
70         fMCEvent(NULL),
71         fMCStack(NULL),
72         fCutFolder(NULL),
73         fESDList(NULL),
74         fBackList(NULL),
75         fMotherList(NULL),
76         fPhotonDCAList(NULL),
77         fMesonDCAList(NULL),
78         fTrueList(NULL),
79         fMCList(NULL),
80         fHeaderNameList(NULL),
81         fTagOutputList(NULL),
82         fOutputContainer(NULL),
83         fReaderGammas(NULL),
84         fGammaCandidates(NULL),
85         fClusterCandidates(NULL),
86         fCutArray(NULL),
87         fConversionCuts(NULL),
88         fClusterCutArray(NULL),
89         fCaloPhotonCuts(NULL),
90         fMesonCutArray(NULL),
91         fMesonCuts(NULL),
92         fHistoConvGammaPt(NULL),
93         fHistoConvGammaR(NULL),
94         fHistoConvGammaEta(NULL),
95         fTreeConvGammaPtDcazCat(NULL),
96         fPtGamma(0),
97         fDCAzPhoton(0),
98         fRConvPhoton(0),
99         fEtaPhoton(0),
100         fCharCatPhoton(0),
101         fCharPhotonMCInfo(0),
102         fHistoMotherInvMassPt(NULL),
103         fSparseMotherInvMassPtZM(NULL),
104         fHistoMotherBackInvMassPt(NULL),
105         fSparseMotherBackInvMassPtZM(NULL),
106         fHistoMotherInvMassEalpha(NULL),
107         fHistoMotherPi0PtY(NULL),
108         fHistoMotherEtaPtY(NULL),
109         fHistoMotherPi0PtAlpha(NULL),
110         fHistoMotherEtaPtAlpha(NULL),
111         fHistoMotherPi0PtOpenAngle(NULL),
112         fHistoMotherEtaPtOpenAngle(NULL),
113         fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
114         fInvMass(0),
115         fPt(0),
116         fDCAzGammaMin(0),
117         fDCAzGammaMax(0),
118         fCharFlag(0),
119         fCharMesonMCInfo(0),
120         fHistoConvGammaUntagged(NULL),
121         fHistoConvGammaTagged(NULL),
122         fHistoConvGammaPi0Tagged(NULL),
123         fHistoConvGammaEtaTagged(NULL),
124         fHistoPhotonPairAll(NULL),
125         fHistoPhotonPairAllGam(NULL),
126         fHistoClusGammaPt(NULL),
127         fHistoMCHeaders(NULL),
128         fHistoMCAllGammaPt(NULL),
129         fHistoMCDecayGammaPi0Pt(NULL),
130         fHistoMCDecayGammaRhoPt(NULL),
131         fHistoMCDecayGammaEtaPt(NULL),
132         fHistoMCDecayGammaOmegaPt(NULL),
133         fHistoMCDecayGammaEtapPt(NULL),
134         fHistoMCDecayGammaPhiPt(NULL),
135         fHistoMCDecayGammaSigmaPt(NULL),
136         fHistoMCConvGammaPt(NULL),
137         fHistoMCConvGammaR(NULL),
138         fHistoMCConvGammaEta(NULL),
139         fHistoMCPi0Pt(NULL),
140         fHistoMCPi0WOWeightPt(NULL),
141         fHistoMCEtaPt(NULL),
142         fHistoMCEtaWOWeightPt(NULL),
143         fHistoMCPi0InAccPt(NULL),
144         fHistoMCEtaInAccPt(NULL),
145         fHistoMCPi0PtY(NULL),
146         fHistoMCEtaPtY(NULL),
147         fHistoMCK0sPt(NULL),
148         fHistoMCK0sWOWeightPt(NULL),
149         fHistoMCK0sPtY(NULL),
150         fHistoMCSecPi0PtvsSource(NULL),
151         fHistoMCSecPi0Source(NULL),
152         fHistoMCSecEtaPt(NULL),
153         fHistoMCSecEtaSource(NULL),
154         fHistoTrueMotherInvMassPt(NULL),
155         fHistoTruePrimaryMotherInvMassPt(NULL),
156         fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL),
157         fProfileTruePrimaryMotherWeightsInvMassPt(NULL),
158         fHistoTruePrimaryPi0MCPtResolPt(NULL),
159         fHistoTruePrimaryEtaMCPtResolPt(NULL),
160         fHistoTrueSecondaryMotherInvMassPt(NULL),
161         fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL),
162         fHistoTrueK0sWithPi0DaughterMCPt(NULL),
163         fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL),
164         fHistoTrueEtaWithPi0DaughterMCPt(NULL),
165         fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL),
166         fHistoTrueLambdaWithPi0DaughterMCPt(NULL),
167         fHistoTrueBckGGInvMassPt(NULL),
168         fHistoTrueBckContInvMassPt(NULL),
169         fHistoTruePi0PtY(NULL),
170         fHistoTrueEtaPtY(NULL),
171         fHistoTruePi0PtAlpha(NULL),
172         fHistoTrueEtaPtAlpha(NULL),
173         fHistoTruePi0PtOpenAngle(NULL),
174         fHistoTrueEtaPtOpenAngle(NULL),
175         fHistoTrueMotherDalitzInvMassPt(NULL),
176         fHistoTrueConvGammaPt(NULL),
177         fHistoTrueConvPi0GammaPt(NULL),
178         fHistoTrueConvGammaEta(NULL),
179         fHistoCombinatorialPt(NULL),
180         fHistoTruePrimaryConvGammaPt(NULL),
181         fHistoTruePrimaryConvGammaESDPtMCPt(NULL),
182         fHistoTrueSecondaryConvGammaPt(NULL),
183         fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL),
184         fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
185         fHistoTrueDalitzPsiPairDeltaPhi(NULL),
186         fHistoTrueGammaPsiPairDeltaPhi(NULL),
187         fHistoTrueClusGammaPt(NULL),
188         fHistoTruePrimaryClusGammaPt(NULL),
189         fHistoTruePrimaryClusGammaESDPtMCPt(NULL),
190         fHistoNEvents(NULL),
191         fHistoNGoodESDTracks(NULL),
192         fHistoNGammaCandidates(NULL),
193         fHistoNGoodESDTracksVsNGammaCanditates(NULL),
194         fHistoNV0Tracks(NULL),
195         fProfileEtaShift(NULL),
196         fEventPlaneAngle(-100),
197         fRandom(0),
198         fNGammaCandidates(0),
199         fUnsmearedPx(NULL),
200         fUnsmearedPy(NULL),
201         fUnsmearedPz(NULL),
202         fUnsmearedE(NULL),
203         fMCStackPos(NULL),
204         fMCStackNeg(NULL),
205         fESDArrayPos(NULL),
206         fESDArrayNeg(NULL),
207         fnCuts(0),
208         fiCut(0),
209         fMoveParticleAccordingToVertex(kTRUE),
210         fIsHeavyIon(0),
211         fDoMesonAnalysis(kTRUE),
212         fDoMesonQA(0),
213         fDoPhotonQA(0),
214         fIsFromMBHeader(kTRUE),
215         fIsMC(kFALSE),
216         fMinE(0.1),
217         fNminCells(2),
218         fEMCm02cut(0.5)
219 {
220   
221 }
222
223 //________________________________________________________________________
224 AliAnalysisTaskGammaConvCalo::AliAnalysisTaskGammaConvCalo(const char *name):
225         AliAnalysisTaskSE(name),
226         fV0Reader(NULL),
227         fBGHandler(NULL),
228         fBGHandlerRP(NULL),
229         fBGClusHandler(NULL),
230         fBGClusHandlerRP(NULL),
231         fInputEvent(NULL),
232         fMCEvent(NULL),
233         fMCStack(NULL),
234         fCutFolder(NULL),
235         fESDList(NULL),
236         fBackList(NULL),
237         fMotherList(NULL),
238         fPhotonDCAList(NULL),
239         fMesonDCAList(NULL),
240         fTrueList(NULL),
241         fMCList(NULL),
242         fHeaderNameList(NULL),
243         fTagOutputList(NULL),
244         fOutputContainer(0),
245         fReaderGammas(NULL),
246         fGammaCandidates(NULL),
247         fClusterCandidates(NULL),
248         fCutArray(NULL),
249         fConversionCuts(NULL),
250         fClusterCutArray(NULL),
251         fCaloPhotonCuts(NULL),
252         fMesonCutArray(NULL),
253         fMesonCuts(NULL),
254         fHistoConvGammaPt(NULL),
255         fHistoConvGammaR(NULL),
256         fHistoConvGammaEta(NULL),
257         fTreeConvGammaPtDcazCat(NULL),
258         fPtGamma(0),
259         fDCAzPhoton(0),
260         fRConvPhoton(0),
261         fEtaPhoton(0),
262         fCharCatPhoton(0),
263         fCharPhotonMCInfo(0),
264         fHistoMotherInvMassPt(NULL),
265         fSparseMotherInvMassPtZM(NULL),
266         fHistoMotherBackInvMassPt(NULL),
267         fSparseMotherBackInvMassPtZM(NULL),
268         fHistoMotherInvMassEalpha(NULL),
269         fHistoMotherPi0PtY(NULL),
270         fHistoMotherEtaPtY(NULL),
271         fHistoMotherPi0PtAlpha(NULL),
272         fHistoMotherEtaPtAlpha(NULL),
273         fHistoMotherPi0PtOpenAngle(NULL),
274         fHistoMotherEtaPtOpenAngle(NULL),
275         fTreeMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
276         fInvMass(0),
277         fPt(0),
278         fDCAzGammaMin(0),
279         fDCAzGammaMax(0),
280         fCharFlag(0),
281         fCharMesonMCInfo(0),
282         fHistoConvGammaUntagged(NULL),
283         fHistoConvGammaTagged(NULL),
284         fHistoConvGammaPi0Tagged(NULL),
285         fHistoConvGammaEtaTagged(NULL),
286         fHistoPhotonPairAll(NULL),
287         fHistoPhotonPairAllGam(NULL),
288         fHistoClusGammaPt(NULL),
289         fHistoMCHeaders(NULL),
290         fHistoMCAllGammaPt(NULL),
291         fHistoMCDecayGammaPi0Pt(NULL),
292         fHistoMCDecayGammaRhoPt(NULL),
293         fHistoMCDecayGammaEtaPt(NULL),
294         fHistoMCDecayGammaOmegaPt(NULL),
295         fHistoMCDecayGammaEtapPt(NULL),
296         fHistoMCDecayGammaPhiPt(NULL),
297         fHistoMCDecayGammaSigmaPt(NULL),
298         fHistoMCConvGammaPt(NULL),
299         fHistoMCConvGammaR(NULL),
300         fHistoMCConvGammaEta(NULL),
301         fHistoMCPi0Pt(NULL),
302         fHistoMCPi0WOWeightPt(NULL),
303         fHistoMCEtaPt(NULL),
304         fHistoMCEtaWOWeightPt(NULL),
305         fHistoMCPi0InAccPt(NULL),
306         fHistoMCEtaInAccPt(NULL),
307         fHistoMCPi0PtY(NULL),
308         fHistoMCEtaPtY(NULL),
309         fHistoMCK0sPt(NULL),
310         fHistoMCK0sWOWeightPt(NULL),
311         fHistoMCK0sPtY(NULL),
312         fHistoMCSecPi0PtvsSource(NULL),
313         fHistoMCSecPi0Source(NULL),
314         fHistoMCSecEtaPt(NULL),
315         fHistoMCSecEtaSource(NULL),
316         fHistoTrueMotherInvMassPt(NULL),
317         fHistoTruePrimaryMotherInvMassPt(NULL),
318         fHistoTruePrimaryMotherW0WeightingInvMassPt(NULL),
319         fProfileTruePrimaryMotherWeightsInvMassPt(NULL),
320         fHistoTruePrimaryPi0MCPtResolPt(NULL),
321         fHistoTruePrimaryEtaMCPtResolPt(NULL),
322         fHistoTrueSecondaryMotherInvMassPt(NULL),
323         fHistoTrueSecondaryMotherFromK0sInvMassPt(NULL),
324         fHistoTrueK0sWithPi0DaughterMCPt(NULL),
325         fHistoTrueSecondaryMotherFromEtaInvMassPt(NULL),
326         fHistoTrueEtaWithPi0DaughterMCPt(NULL),
327         fHistoTrueSecondaryMotherFromLambdaInvMassPt(NULL),
328         fHistoTrueLambdaWithPi0DaughterMCPt(NULL),
329         fHistoTrueBckGGInvMassPt(NULL),
330         fHistoTrueBckContInvMassPt(NULL),
331         fHistoTruePi0PtY(NULL),
332         fHistoTrueEtaPtY(NULL),
333         fHistoTruePi0PtAlpha(NULL),
334         fHistoTrueEtaPtAlpha(NULL),
335         fHistoTruePi0PtOpenAngle(NULL),
336         fHistoTrueEtaPtOpenAngle(NULL),
337         fHistoTrueMotherDalitzInvMassPt(NULL),
338         fHistoTrueConvGammaPt(NULL),
339         fHistoTrueConvPi0GammaPt(NULL),
340         fHistoTrueConvGammaEta(NULL),
341         fHistoCombinatorialPt(NULL),
342         fHistoTruePrimaryConvGammaPt(NULL),
343         fHistoTruePrimaryConvGammaESDPtMCPt(NULL),
344         fHistoTrueSecondaryConvGammaPt(NULL),
345         fHistoTrueSecondaryConvGammaFromXFromK0sPt(NULL),
346         fHistoTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
347         fHistoTrueDalitzPsiPairDeltaPhi(NULL),
348         fHistoTrueGammaPsiPairDeltaPhi(NULL),
349         fHistoTrueClusGammaPt(NULL),
350         fHistoTruePrimaryClusGammaPt(NULL),
351         fHistoTruePrimaryClusGammaESDPtMCPt(NULL),
352         fHistoNEvents(NULL),
353         fHistoNGoodESDTracks(NULL),
354         fHistoNGammaCandidates(NULL),
355         fHistoNGoodESDTracksVsNGammaCanditates(NULL),
356         fHistoNV0Tracks(NULL),
357         fProfileEtaShift(NULL),
358         fEventPlaneAngle(-100),
359         fRandom(0),
360         fNGammaCandidates(0),
361         fUnsmearedPx(NULL),
362         fUnsmearedPy(NULL),
363         fUnsmearedPz(NULL),
364         fUnsmearedE(NULL),
365         fMCStackPos(NULL),
366         fMCStackNeg(NULL),
367         fESDArrayPos(NULL),
368         fESDArrayNeg(NULL),
369         fnCuts(0),
370         fiCut(0),
371         fMoveParticleAccordingToVertex(kTRUE),
372         fIsHeavyIon(0),
373         fDoMesonAnalysis(kTRUE),
374         fDoMesonQA(0),
375         fDoPhotonQA(0),
376         fIsFromMBHeader(kTRUE),
377         fIsMC(kFALSE),
378         fMinE(0.1),
379         fNminCells(2),
380         fEMCm02cut(0.5)
381 {
382   // Define output slots here
383   DefineOutput(1, TList::Class());
384 }
385
386 AliAnalysisTaskGammaConvCalo::~AliAnalysisTaskGammaConvCalo()
387 {
388         if(fGammaCandidates){
389                 delete fGammaCandidates;
390                 fGammaCandidates = 0x0;
391         }
392         if(fClusterCandidates){
393                 delete fClusterCandidates;
394                 fClusterCandidates = 0x0;
395         }
396         if(fBGHandler){
397                 delete[] fBGHandler;
398                 fBGHandler = 0x0;
399         }
400         if(fBGHandlerRP){
401                 delete[] fBGHandlerRP;
402                 fBGHandlerRP = 0x0;
403         }
404         if(fBGClusHandler){
405                 delete[] fBGClusHandler;
406                 fBGClusHandler = 0x0;
407         }
408         if(fBGClusHandlerRP){
409                 delete[] fBGClusHandlerRP;
410                 fBGClusHandlerRP = 0x0;
411         }
412 }
413 //___________________________________________________________
414 void AliAnalysisTaskGammaConvCalo::InitBack(){
415         
416         const Int_t nDim = 4;
417         Int_t nBins[nDim] = {800,250,7,4};
418         Double_t xMin[nDim] = {0,0, 0,0};
419         Double_t xMax[nDim] = {0.8,25,7,4};
420         
421         fSparseMotherInvMassPtZM = new THnSparseF*[fnCuts];
422         fSparseMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
423         
424         fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
425         fBGHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
426
427         fBGClusHandler = new AliGammaConversionAODBGHandler*[fnCuts];
428         fBGClusHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
429         
430         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
431                 if (((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
432                         TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
433                         TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
434                         TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
435                         
436                         Int_t collisionSystem = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(0,1));
437                         Int_t centMin = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(1,1));
438                         Int_t centMax = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(2,1));
439                         
440                         if(collisionSystem == 1 || collisionSystem == 2 ||
441                                 collisionSystem == 5 || collisionSystem == 8 ||
442                                 collisionSystem == 9){
443                                 centMin = centMin*10;
444                                 centMax = centMax*10;
445                                 if(centMax ==0 && centMax!=centMin) centMax=100;
446                         } else if(collisionSystem == 3 || collisionSystem == 6){
447                                 centMin = centMin*5;
448                                 centMax = centMax*5;
449                         } else if(collisionSystem == 4 || collisionSystem == 7){
450                                 centMin = ((centMin*5)+45);
451                                 centMax = ((centMax*5)+45);
452                         }
453                         
454                         fBackList[iCut] = new TList();
455                         fBackList[iCut]->SetName(Form("%s_%s_%s Back histograms",cutstring.Data(),cutstringCalo.Data(), cutstringMeson.Data()));
456                         fBackList[iCut]->SetOwner(kTRUE);
457                         fCutFolder[iCut]->Add(fBackList[iCut]);
458                         
459                         fSparseMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
460                         fBackList[iCut]->Add(fSparseMotherBackInvMassPtZM[iCut]);
461                         
462                         fMotherList[iCut] = new TList();
463                         fMotherList[iCut]->SetName(Form("%s_%s_%s Mother histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
464                         fMotherList[iCut]->SetOwner(kTRUE);
465                         fCutFolder[iCut]->Add(fMotherList[iCut]);
466                         
467                         fSparseMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
468                         fMotherList[iCut]->Add(fSparseMotherInvMassPtZM[iCut]);
469                         
470                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
471                                 fBGHandler[iCut] = new AliGammaConversionAODBGHandler(
472                                                                                                                                         collisionSystem,centMin,centMax,
473                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
474                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity());
475                                 fBGClusHandler[iCut] = new AliGammaConversionAODBGHandler(
476                                                                                                                                         collisionSystem,centMin,centMax,
477                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
478                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity());
479                                 fBGHandlerRP[iCut] = NULL;
480                         } else{
481                                 fBGHandlerRP[iCut] = new AliConversionAODBGHandlerRP(
482                                                                                                                                         ((AliConversionCuts*)fCutArray->At(fiCut))->IsHeavyIon(),
483                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
484                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
485                                 fBGClusHandlerRP[iCut] = new AliConversionAODBGHandlerRP(
486                                                                                                                                         ((AliConversionCuts*)fCutArray->At(fiCut))->IsHeavyIon(),
487                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
488                                                                                                                                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
489                                 fBGHandler[iCut] = NULL;
490                         }
491                 }
492         }
493 }
494 //________________________________________________________________________
495 void AliAnalysisTaskGammaConvCalo::UserCreateOutputObjects(){
496   
497         // Create histograms
498         if(fOutputContainer != NULL){
499                 delete fOutputContainer;
500                 fOutputContainer = NULL;
501         }
502         if(fOutputContainer == NULL){
503                 fOutputContainer = new TList();
504                 fOutputContainer->SetOwner(kTRUE);
505         }
506   
507         // Array of current cut's gammas
508         fGammaCandidates = new TList();
509         fClusterCandidates = new TList();
510   
511         fCutFolder = new TList*[fnCuts];
512         fESDList = new TList*[fnCuts];
513         fBackList = new TList*[fnCuts];
514         fMotherList = new TList*[fnCuts];
515         fHistoNEvents = new TH1I*[fnCuts];
516         fHistoNGoodESDTracks = new TH1I*[fnCuts];
517         fHistoNGammaCandidates = new TH1I*[fnCuts];
518         fHistoNGoodESDTracksVsNGammaCanditates = new TH2F*[fnCuts];
519         fHistoNV0Tracks = new TH1I*[fnCuts];
520         fProfileEtaShift = new TProfile*[fnCuts];
521         fHistoConvGammaPt = new TH1F*[fnCuts];
522   
523         if (fDoPhotonQA == 2){
524                 fPhotonDCAList = new TList*[fnCuts];
525                 fTreeConvGammaPtDcazCat = new TTree*[fnCuts];
526         }
527         if (fDoPhotonQA > 0){
528                 fHistoConvGammaR = new TH1F*[fnCuts];
529                 fHistoConvGammaEta = new TH1F*[fnCuts];
530         }
531         
532         if(fDoMesonAnalysis){
533                 fHistoMotherInvMassPt = new TH2F*[fnCuts];
534                 fHistoMotherBackInvMassPt = new TH2F*[fnCuts];
535                 fHistoMotherInvMassEalpha = new TH2F*[fnCuts];
536                 if (fDoMesonQA == 2){
537                         fMesonDCAList = new TList*[fnCuts];
538                         fTreeMesonsInvMassPtDcazMinDcazMaxFlag = new TTree*[fnCuts];
539                 }
540                 if (fDoMesonQA > 0){
541                         fHistoMotherPi0PtY =  new TH2F*[fnCuts];
542                         fHistoMotherEtaPtY =  new TH2F*[fnCuts];
543                         fHistoMotherPi0PtAlpha =  new TH2F*[fnCuts];
544                         fHistoMotherEtaPtAlpha =  new TH2F*[fnCuts];
545                         fHistoMotherPi0PtOpenAngle =  new TH2F*[fnCuts];
546                         fHistoMotherEtaPtOpenAngle =  new TH2F*[fnCuts];
547                 }
548         }
549         fTagOutputList = new TList*[fnCuts];
550         
551         fHistoConvGammaUntagged = new TH1F*[fnCuts];
552         fHistoConvGammaTagged = new TH1F*[fnCuts];
553         fHistoConvGammaPi0Tagged = new TH1F*[fnCuts];
554         fHistoConvGammaEtaTagged = new TH1F*[fnCuts];
555         fHistoPhotonPairAll = new TH2F*[fnCuts];
556         fHistoPhotonPairAllGam = new TH2F*[fnCuts];
557         
558         fHistoClusGammaPt = new TH1F*[fnCuts];
559         
560         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
561
562                 TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
563                 TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
564                 TString cutstringMeson = "NoMesonCut";
565                 if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
566     
567                 fCutFolder[iCut] = new TList();
568                 fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
569                 fCutFolder[iCut]->SetOwner(kTRUE);
570                 fOutputContainer->Add(fCutFolder[iCut]);
571                 fESDList[iCut] = new TList();
572                 fESDList[iCut]->SetName(Form("%s_%s_%s ESD histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
573                 fESDList[iCut]->SetOwner(kTRUE);
574                 fCutFolder[iCut]->Add(fESDList[iCut]);
575     
576                 fHistoNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
577                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
578                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
579                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
580                 if (((AliConversionCuts*)fCutArray->At(iCut))->IsSpecialTrigger() == 4 ){
581                         TString TriggerNames = "Not Trigger: ";
582                         TriggerNames = TriggerNames+ ( (AliConversionCuts*)fCutArray->At(iCut))->GetSpecialTriggerName();
583                         fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
584                 } else {
585                         fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
586                 }
587                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
588                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
589                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
590                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
591                 fHistoNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
592                 fESDList[iCut]->Add(fHistoNEvents[iCut]);
593                 
594                 if(fIsHeavyIon == 1) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",4000,0,4000);
595                 else if(fIsHeavyIon == 2) fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",400,0,400);
596                 else fHistoNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
597                 fESDList[iCut]->Add(fHistoNGoodESDTracks[iCut]);
598                 if(fIsHeavyIon == 1) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",100,0,100);
599                 else if(fIsHeavyIon == 2) fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
600                 else fHistoNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
601                 fESDList[iCut]->Add(fHistoNGammaCandidates[iCut]);
602                 if(fIsHeavyIon == 1) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",4000,0,4000,100,0,100);
603                 else if(fIsHeavyIon == 2) fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",400,0,400,50,0,50);
604                 else fHistoNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",200,0,200,50,0,50);
605                 fESDList[iCut]->Add(fHistoNGoodESDTracksVsNGammaCanditates[iCut]);
606     
607                 
608                 if(fIsHeavyIon == 1) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",30000,0,30000);
609                 else if(fIsHeavyIon == 2) fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2500,0,2500);
610                 else fHistoNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",1500,0,1500);
611                 fESDList[iCut]->Add(fHistoNV0Tracks[iCut]);
612                 fProfileEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
613                 fESDList[iCut]->Add(fProfileEtaShift[iCut]);
614                 fHistoConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
615                 fESDList[iCut]->Add(fHistoConvGammaPt[iCut]);
616     
617                 if (fDoPhotonQA == 2){
618                         fPhotonDCAList[iCut] = new TList();
619                         fPhotonDCAList[iCut]->SetName(Form("%s_%s_%s Photon DCA tree",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
620                         fPhotonDCAList[iCut]->SetOwner(kTRUE);
621                         fCutFolder[iCut]->Add(fPhotonDCAList[iCut]);
622       
623                         fTreeConvGammaPtDcazCat[iCut] = new TTree("ESD_ConvGamma_Pt_Dcaz_R_Eta","ESD_ConvGamma_Pt_Dcaz_R_Eta_Cat");
624                         fTreeConvGammaPtDcazCat[iCut]->Branch("Pt",&fPtGamma,"fPtGamma/F");
625                         fTreeConvGammaPtDcazCat[iCut]->Branch("DcaZPhoton",&fDCAzPhoton,"fDCAzPhoton/F");
626       //          fTreeConvGammaPtDcazCat[iCut]->Branch("R",&fRConvPhoton,"fRConvPhoton/F");
627       //          fTreeConvGammaPtDcazCat[iCut]->Branch("Eta",&fEtaPhoton,"fEtaPhoton/F");
628                         
629                         fTreeConvGammaPtDcazCat[iCut]->Branch("cat",&fCharCatPhoton,"fCharCatPhoton/b");
630                         if(fIsMC){
631                                 fTreeConvGammaPtDcazCat[iCut]->Branch("photonMCInfo",&fCharPhotonMCInfo,"fCharPhotonMCInfo/b");
632                         }
633                         fPhotonDCAList[iCut]->Add(fTreeConvGammaPtDcazCat[iCut]);
634                 }
635     
636                 if (fDoPhotonQA > 0){
637                         fHistoConvGammaR[iCut] = new TH1F("ESD_ConvGamma_R","ESD_ConvGamma_R",800,0,200);
638                         fESDList[iCut]->Add(fHistoConvGammaR[iCut]);
639                         fHistoConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",2000,-2,2);
640                         fESDList[iCut]->Add(fHistoConvGammaEta[iCut]);
641                 }
642
643                 fTagOutputList[iCut] = new TList();
644                 fTagOutputList[iCut]->SetName(Form("Tagging Output Cut Number %s_%s_%s",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
645                 fTagOutputList[iCut]->SetOwner(1);
646                 fCutFolder[iCut]->Add(fTagOutputList[iCut]);
647                 
648                 const Int_t nptbins = 200;
649                 const Double_t ptmin = 0.;
650                 const Double_t ptmax = 20.;
651                 
652                 const Int_t nmbins = 180;
653                 const Double_t mmin = 0.;
654                 const Double_t mmax = 0.9;
655                 
656                 // photon candidates
657                 // this is maybe not necessary ...
658                         
659                 fHistoConvGammaUntagged[iCut] = new TH1F("ConvGammaUntagged","",nptbins,ptmin,ptmax);
660                 fHistoConvGammaUntagged[iCut]->SetXTitle("p_{T} (GeV/c)");
661                 fTagOutputList[iCut]->Add(fHistoConvGammaUntagged[iCut]);
662                 
663                 fHistoConvGammaTagged[iCut] = new TH1F("ConvGammaTagged","",nptbins,ptmin,ptmax);
664                 fHistoConvGammaTagged[iCut]->SetXTitle("p_{T} (GeV/c)");
665                 fTagOutputList[iCut]->Add(fHistoConvGammaTagged[iCut]);
666                 
667                 fHistoConvGammaPi0Tagged[iCut] = new TH1F("ConvGammaPi0Tagged","",nptbins,ptmin,ptmax);
668                 fHistoConvGammaPi0Tagged[iCut]->SetXTitle("p_{T} (GeV/c)");
669                 fTagOutputList[iCut]->Add(fHistoConvGammaPi0Tagged[iCut]);
670                 
671                 fHistoConvGammaEtaTagged[iCut] = new TH1F("ConvGammaEtaTagged","",nptbins,ptmin,ptmax);
672                 fHistoConvGammaEtaTagged[iCut]->SetXTitle("p_{T} (GeV/c)");
673                 fTagOutputList[iCut]->Add(fHistoConvGammaEtaTagged[iCut]);
674                 
675                 // pairs
676                 fHistoPhotonPairAll[iCut] = new TH2F("PhotonPairAll","",nmbins,mmin,mmax,nptbins,ptmin,ptmax);
677                 fHistoPhotonPairAll[iCut]->SetXTitle("M_{inv} (GeV/cc)");
678                 fHistoPhotonPairAll[iCut]->SetYTitle("p_{T} (GeV/c)");
679                 fTagOutputList[iCut]->Add(fHistoPhotonPairAll[iCut]);
680                 
681                 fHistoPhotonPairAllGam[iCut] = new TH2F("PhotonPairAllGammaConvPt","",nmbins,mmin,mmax,nptbins,ptmin,ptmax);
682                 fHistoPhotonPairAllGam[iCut]->SetXTitle("M_{inv} (GeV/cc)");
683                 fHistoPhotonPairAllGam[iCut]->SetYTitle("#gamma^{conv} p_{T} (GeV/c)");
684                 fTagOutputList[iCut]->Add(fHistoPhotonPairAllGam[iCut]);
685         
686                 fHistoClusGammaPt[iCut] = new TH1F("ClusGamma_Pt","ClusGamma_Pt",250,0,25);
687                 fTagOutputList[iCut]->Add(fHistoClusGammaPt[iCut]);
688
689                 
690                 if(fDoMesonAnalysis){
691                         fHistoMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",800,0,0.8,250,0,25);
692                         fESDList[iCut]->Add(fHistoMotherInvMassPt[iCut]);
693                         fHistoMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",800,0,0.8,250,0,25);
694                         fESDList[iCut]->Add(fHistoMotherBackInvMassPt[iCut]);
695                         fHistoMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",800,0,0.8,250,0,25);
696                         fESDList[iCut]->Add(fHistoMotherInvMassEalpha[iCut]);
697                         if (fDoMesonQA == 2){
698                                 fMesonDCAList[iCut] = new TList();
699                                 fMesonDCAList[iCut]->SetName(Form("%s_%s_%s Meson DCA tree",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
700                                 fMesonDCAList[iCut]->SetOwner(kTRUE);
701                                 fCutFolder[iCut]->Add(fMesonDCAList[iCut]);
702                                 
703                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut] = new TTree("ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag","ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag");
704                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("InvMass",&fInvMass,"fInvMass/F");
705                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("Pt",&fPt,"fPt/F");
706                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMin",&fDCAzGammaMin,"fDCAzGammaMin/F");
707                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMax",&fDCAzGammaMax,"fDCAzGammaMax/F");
708                                 fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("kind",&fCharFlag,"fCharFlag/b");
709                                 if(fIsMC){
710                                         fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("mesonMCInfo",&fCharMesonMCInfo,"fCharMesonMCInfo/b");
711                                 }
712                                 fMesonDCAList[iCut]->Add(fTreeMesonsInvMassPtDcazMinDcazMaxFlag[iCut]);
713                                 
714                         }
715                         if (fDoMesonQA > 0 ){
716                                 fHistoMotherPi0PtY[iCut] = new TH2F("ESD_MotherPi0_Pt_Y","ESD_MotherPi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
717                                 SetLogBinningXTH2(fHistoMotherPi0PtY[iCut]);
718                                 fESDList[iCut]->Add(fHistoMotherPi0PtY[iCut]);
719                                 fHistoMotherEtaPtY[iCut] = new TH2F("ESD_MotherEta_Pt_Y","ESD_MotherEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
720                                 SetLogBinningXTH2(fHistoMotherEtaPtY[iCut]);
721                                 fESDList[iCut]->Add(fHistoMotherEtaPtY[iCut]);
722                                 fHistoMotherPi0PtAlpha[iCut] = new TH2F("ESD_MotherPi0_Pt_Alpha","ESD_MotherPi0_Pt_Alpha",150,0.03,15.,100,0,1);
723                                 SetLogBinningXTH2(fHistoMotherPi0PtAlpha[iCut]);
724                                 fESDList[iCut]->Add(fHistoMotherPi0PtAlpha[iCut]);
725                                 fHistoMotherEtaPtAlpha[iCut] = new TH2F("ESD_MotherEta_Pt_Alpha","ESD_MotherEta_Pt_Alpha",150,0.03,15.,100,0,1);
726                                 SetLogBinningXTH2(fHistoMotherEtaPtAlpha[iCut]);
727                                 fESDList[iCut]->Add(fHistoMotherEtaPtAlpha[iCut]);
728                                 fHistoMotherPi0PtOpenAngle[iCut] = new TH2F("ESD_MotherPi0_Pt_OpenAngle","ESD_MotherPi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
729                                 SetLogBinningXTH2(fHistoMotherPi0PtOpenAngle[iCut]);
730                                 fESDList[iCut]->Add(fHistoMotherPi0PtOpenAngle[iCut]);
731                                 fHistoMotherEtaPtOpenAngle[iCut] = new TH2F("ESD_MotherEta_Pt_OpenAngle","ESD_MotherEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
732                                 SetLogBinningXTH2(fHistoMotherEtaPtOpenAngle[iCut]);
733                                 fESDList[iCut]->Add(fHistoMotherEtaPtOpenAngle[iCut]);
734                         }
735                 }    
736         }
737         if(fDoMesonAnalysis){
738                 InitBack(); // Init Background Handler
739         }
740   
741         if(fIsMC){
742                 // MC Histogramms
743                 fMCList = new TList*[fnCuts];
744                 // True Histogramms
745                 fTrueList = new TList*[fnCuts];
746                 // Selected Header List
747                 fHeaderNameList = new TList*[fnCuts];
748                 fHistoMCHeaders = new TH1I*[fnCuts];
749                 fHistoMCAllGammaPt = new TH1F*[fnCuts];
750                 fHistoMCDecayGammaPi0Pt = new TH1F*[fnCuts];
751                 fHistoMCDecayGammaRhoPt = new TH1F*[fnCuts];
752                 fHistoMCDecayGammaEtaPt = new TH1F*[fnCuts];
753                 fHistoMCDecayGammaOmegaPt = new TH1F*[fnCuts];
754                 fHistoMCDecayGammaEtapPt = new TH1F*[fnCuts];
755                 fHistoMCDecayGammaPhiPt = new TH1F*[fnCuts];
756                 fHistoMCDecayGammaSigmaPt = new TH1F*[fnCuts];
757                 fHistoMCConvGammaPt = new TH1F*[fnCuts];
758                 fHistoTrueConvGammaPt = new TH1F*[fnCuts];
759                 fHistoTrueConvPi0GammaPt = new TH1F*[fnCuts];
760     
761                 fHistoCombinatorialPt = new TH2F*[fnCuts];
762                 fHistoTruePrimaryConvGammaPt = new TH1F*[fnCuts];
763                 fHistoTruePrimaryConvGammaESDPtMCPt = new TH2F*[fnCuts];
764                 fHistoTrueSecondaryConvGammaPt = new TH1F*[fnCuts];
765                 fHistoTrueSecondaryConvGammaFromXFromK0sPt = new TH1F*[fnCuts];
766                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt = new TH1F*[fnCuts];
767     
768                 fHistoTrueDalitzPsiPairDeltaPhi= new TH2F*[fnCuts];
769                 fHistoTrueGammaPsiPairDeltaPhi= new TH2F*[fnCuts];
770
771                 fHistoTrueClusGammaPt = new TH1F*[fnCuts];
772                 fHistoTruePrimaryClusGammaPt = new TH1F*[fnCuts];
773                 fHistoTruePrimaryClusGammaESDPtMCPt = new TH2F*[fnCuts];
774
775                 if (fDoPhotonQA > 0){
776                         fHistoMCConvGammaR = new TH1F*[fnCuts];
777                         fHistoMCConvGammaEta = new TH1F*[fnCuts];
778                         fHistoTrueConvGammaEta = new TH1F*[fnCuts];
779                 }
780     
781                 if(fDoMesonAnalysis){
782                         fHistoMCPi0Pt = new TH1F*[fnCuts];
783                         fHistoMCPi0WOWeightPt = new TH1F*[fnCuts];
784                         fHistoMCEtaPt = new TH1F*[fnCuts];
785                         fHistoMCEtaWOWeightPt = new TH1F*[fnCuts];
786                         fHistoMCPi0InAccPt = new TH1F*[fnCuts];
787                         fHistoMCEtaInAccPt = new TH1F*[fnCuts];
788       
789                         fHistoTrueMotherInvMassPt = new TH2F*[fnCuts];
790                         fHistoTruePrimaryMotherInvMassPt = new TH2F*[fnCuts];
791                         fHistoTruePrimaryMotherW0WeightingInvMassPt = new TH2F*[fnCuts];
792                         fProfileTruePrimaryMotherWeightsInvMassPt = new TProfile2D*[fnCuts];
793                         fHistoTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
794                         fHistoTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
795                         fHistoTrueSecondaryMotherFromEtaInvMassPt = new TH2F*[fnCuts];
796                         fHistoTrueSecondaryMotherFromLambdaInvMassPt = new TH2F*[fnCuts];
797                         fHistoTrueMotherDalitzInvMassPt = new TH2F*[fnCuts];
798                         if (fDoMesonQA > 0){
799                                 fHistoMCPi0PtY = new TH2F*[fnCuts];
800                                 fHistoMCEtaPtY = new TH2F*[fnCuts];
801                                 fHistoMCK0sPt = new TH1F*[fnCuts];
802                                 fHistoMCK0sWOWeightPt = new TH1F*[fnCuts];
803                                 fHistoMCK0sPtY = new TH2F*[fnCuts];
804                                 fHistoMCSecPi0PtvsSource= new TH2F*[fnCuts];
805                                 fHistoMCSecPi0Source = new TH1F*[fnCuts];
806                                 fHistoMCSecEtaPt = new TH1F*[fnCuts];
807                                 fHistoMCSecEtaSource = new TH1F*[fnCuts];
808                                 fHistoTruePrimaryPi0MCPtResolPt = new TH2F*[fnCuts];
809                                 fHistoTruePrimaryEtaMCPtResolPt = new TH2F*[fnCuts];
810                                 fHistoTrueK0sWithPi0DaughterMCPt = new TH1F*[fnCuts];
811                                 fHistoTrueEtaWithPi0DaughterMCPt = new TH1F*[fnCuts];
812                                 fHistoTrueLambdaWithPi0DaughterMCPt = new TH1F*[fnCuts];
813                                 fHistoTrueBckGGInvMassPt = new TH2F*[fnCuts];
814                                 fHistoTrueBckContInvMassPt = new TH2F*[fnCuts];
815                                 fHistoTruePi0PtY = new TH2F*[fnCuts];
816                                 fHistoTrueEtaPtY = new TH2F*[fnCuts];
817                                 fHistoTruePi0PtAlpha = new TH2F*[fnCuts];
818                                 fHistoTrueEtaPtAlpha = new TH2F*[fnCuts];
819                                 fHistoTruePi0PtOpenAngle = new TH2F*[fnCuts];
820                                 fHistoTrueEtaPtOpenAngle = new TH2F*[fnCuts];
821                         }
822                 }
823     
824     
825     
826                 for(Int_t iCut = 0; iCut<fnCuts;iCut++){
827                         TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
828                         TString cutstringCalo = ((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutNumber();
829                         TString cutstringMeson = "NoMesonCut";
830                         if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
831       
832                         fMCList[iCut] = new TList();
833                         fMCList[iCut]->SetName(Form("%s_%s_%s MC histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
834                         fMCList[iCut]->SetOwner(kTRUE);
835                         fCutFolder[iCut]->Add(fMCList[iCut]);
836                         fHistoMCHeaders[iCut] = new TH1I("MC_Headers","MC_Headers",20,0,20);
837                         fMCList[iCut]->Add(fHistoMCHeaders[iCut]);
838                         fHistoMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
839                         fMCList[iCut]->Add(fHistoMCAllGammaPt[iCut]);
840                         fHistoMCDecayGammaPi0Pt[iCut] = new TH1F("MC_DecayGammaPi0_Pt","MC_DecayGammaPi0_Pt",250,0,25);
841                         fMCList[iCut]->Add(fHistoMCDecayGammaPi0Pt[iCut]);
842                         fHistoMCDecayGammaRhoPt[iCut] = new TH1F("MC_DecayGammaRho_Pt","MC_DecayGammaRho_Pt",250,0,25);
843                         fMCList[iCut]->Add(fHistoMCDecayGammaRhoPt[iCut]);
844                         fHistoMCDecayGammaEtaPt[iCut] = new TH1F("MC_DecayGammaEta_Pt","MC_DecayGammaEta_Pt",250,0,25);
845                         fMCList[iCut]->Add(fHistoMCDecayGammaEtaPt[iCut]);
846                         fHistoMCDecayGammaOmegaPt[iCut] = new TH1F("MC_DecayGammaOmega_Pt","MC_DecayGammaOmmega_Pt",250,0,25);
847                         fMCList[iCut]->Add(fHistoMCDecayGammaOmegaPt[iCut]);
848                         fHistoMCDecayGammaEtapPt[iCut] = new TH1F("MC_DecayGammaEtap_Pt","MC_DecayGammaEtap_Pt",250,0,25);
849                         fMCList[iCut]->Add(fHistoMCDecayGammaEtapPt[iCut]);
850                         fHistoMCDecayGammaPhiPt[iCut] = new TH1F("MC_DecayGammaPhi_Pt","MC_DecayGammaPhi_Pt",250,0,25);
851                         fMCList[iCut]->Add(fHistoMCDecayGammaPhiPt[iCut]);
852                         fHistoMCDecayGammaSigmaPt[iCut] = new TH1F("MC_DecayGammaSigma_Pt","MC_DecayGammaSigma_Pt",250,0,25);
853                         fMCList[iCut]->Add(fHistoMCDecayGammaSigmaPt[iCut]);
854                         fHistoMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
855                         fMCList[iCut]->Add(fHistoMCConvGammaPt[iCut]);
856       
857                         if (fDoPhotonQA > 0){
858                                 fHistoMCConvGammaR[iCut] = new TH1F("MC_ConvGamma_R","MC_ConvGamma_R",800,0,200);
859                                 fMCList[iCut]->Add(fHistoMCConvGammaR[iCut]);
860                                 fHistoMCConvGammaEta[iCut] = new TH1F("MC_ConvGamma_Eta","MC_ConvGamma_Eta",2000,-2,2);
861                                 fMCList[iCut]->Add(fHistoMCConvGammaEta[iCut]);
862                         }
863       
864                         if(fDoMesonAnalysis){
865                                 fHistoMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
866                                 fHistoMCPi0Pt[iCut]->Sumw2();
867                                 fMCList[iCut]->Add(fHistoMCPi0Pt[iCut]);
868                                 fHistoMCPi0WOWeightPt[iCut] = new TH1F("MC_Pi0_WOWeights_Pt","MC_Pi0_WOWeights_Pt",250,0,25);
869                                 fHistoMCPi0WOWeightPt[iCut]->Sumw2();
870                                 fMCList[iCut]->Add(fHistoMCPi0WOWeightPt[iCut]);
871                                 
872                                 fHistoMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
873                                 fHistoMCEtaPt[iCut]->Sumw2();
874                                 fMCList[iCut]->Add(fHistoMCEtaPt[iCut]);
875                                 fHistoMCEtaWOWeightPt[iCut] = new TH1F("MC_Eta_WOWeights_Pt","MC_Eta_WOWeights_Pt",250,0,25);
876                                 fHistoMCEtaWOWeightPt[iCut]->Sumw2();
877                                 fMCList[iCut]->Add(fHistoMCEtaWOWeightPt[iCut]);
878                                 fHistoMCPi0InAccPt[iCut] = new TH1F("MC_Pi0InAcc_Pt","MC_Pi0InAcc_Pt",250,0,25);
879                                 fHistoMCPi0InAccPt[iCut]->Sumw2();
880                                 fMCList[iCut]->Add(fHistoMCPi0InAccPt[iCut]);
881                                 fHistoMCEtaInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
882                                 fHistoMCEtaInAccPt[iCut]->Sumw2();
883                                 fMCList[iCut]->Add(fHistoMCEtaInAccPt[iCut]);
884                                 if (fDoMesonQA > 0){
885                                         fHistoMCPi0PtY[iCut] = new TH2F("MC_Pi0_Pt_Y","MC_Pi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
886                                         fHistoMCPi0PtY[iCut]->Sumw2();
887                                         SetLogBinningXTH2(fHistoMCPi0PtY[iCut]);
888                                         fMCList[iCut]->Add(fHistoMCPi0PtY[iCut]);
889                                         fHistoMCEtaPtY[iCut] = new TH2F("MC_Eta_Pt_Y","MC_Eta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
890                                         fHistoMCEtaPtY[iCut]->Sumw2();
891                                         SetLogBinningXTH2(fHistoMCEtaPtY[iCut]);
892                                         fMCList[iCut]->Add(fHistoMCEtaPtY[iCut]);
893                                         fHistoMCK0sPt[iCut] = new TH1F("MC_K0s_Pt","MC_K0s_Pt",150,0,15);
894                                         fHistoMCK0sPt[iCut]->Sumw2();
895                                         fMCList[iCut]->Add(fHistoMCK0sPt[iCut]);
896                                         fHistoMCK0sWOWeightPt[iCut] = new TH1F("MC_K0s_WOWeights_Pt","MC_K0s_WOWeights_Pt",150,0,15);
897                                         fHistoMCK0sWOWeightPt[iCut]->Sumw2();
898                                         fMCList[iCut]->Add(fHistoMCK0sWOWeightPt[iCut]);
899                                         fHistoMCK0sPtY[iCut] = new TH2F("MC_K0s_Pt_Y","MC_K0s_Pt_Y",150,0.03,15.,150,-1.5,1.5);
900                                         fHistoMCK0sPtY[iCut]->Sumw2();
901                                         SetLogBinningXTH2(fHistoMCK0sPtY[iCut]);
902                                         fMCList[iCut]->Add(fHistoMCK0sPtY[iCut]);
903                                         
904                                         fHistoMCSecPi0Source[iCut] = new TH1F("MC_SecPi0_Source","MC_SecPi0_Source",5000,0.,5000);
905                                         fMCList[iCut]->Add(fHistoMCSecPi0Source[iCut]);
906                                         fHistoMCSecEtaSource[iCut] = new TH1F("MC_SecEta_Source","MC_SecEta_Source",5000,0,5000);
907                                         fMCList[iCut]->Add(fHistoMCSecEtaSource[iCut]);
908                                         fHistoMCSecPi0PtvsSource[iCut] = new TH2F("MC_SecPi0_Pt_Source","MC_SecPi0_Pt_Source",250,0.0,25.,16,-0.5,15.5);
909                                         fHistoMCSecPi0PtvsSource[iCut]->Sumw2();
910                                         fMCList[iCut]->Add(fHistoMCSecPi0PtvsSource[iCut]);
911                                         fHistoMCSecEtaPt[iCut] = new TH1F("MC_SecEta_Pt","MC_SecEta_Pt",250,0,25);
912                                         fHistoMCSecEtaPt[iCut]->Sumw2();
913                                         fMCList[iCut]->Add(fHistoMCSecEtaPt[iCut]);
914                                 }
915         
916                         }
917                         fTrueList[iCut] = new TList();
918                         fTrueList[iCut]->SetName(Form("%s_%s_%s True histograms",cutstring.Data(),cutstringCalo.Data(),cutstringMeson.Data()));
919                         fTrueList[iCut]->SetOwner(kTRUE);
920                         fCutFolder[iCut]->Add(fTrueList[iCut]);
921       
922                         fHistoTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
923                         fTrueList[iCut]->Add(fHistoTrueConvGammaPt[iCut]);
924       
925                         fHistoTrueConvPi0GammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
926                         fTrueList[iCut]->Add(fHistoTrueConvPi0GammaPt[iCut]);
927       
928                         fHistoCombinatorialPt[iCut] = new TH2F("ESD_TrueCombinatorial_Pt","ESD_TrueCombinatorial_Pt",250,0,25,16,-0.5,15.5);
929                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 1,"Elec+Elec");
930                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 2,"Elec+Pion");
931                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 3,"Elec+Kaon");
932                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 4,"Elec+Proton");
933                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 5,"Elec+Muon");
934                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 6,"Pion+Pion");
935                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 7,"Pion+Kaon");
936                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 8,"Pion+Proton");
937                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 9,"Pion+Muon");
938                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(10,"Kaon+Kaon");
939                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(11,"Kaon+Proton");
940                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(12,"Kaon+Muon");
941                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(13,"Proton+Proton");
942                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(14,"Proton+Muon");
943                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(15,"Muon+Muon");
944                         fHistoCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(16,"Rest");
945                         fTrueList[iCut]->Add(fHistoCombinatorialPt[iCut]);
946                         fHistoTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt","ESD_TruePrimaryConvGamma_Pt",250,0,25);
947                         fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaPt[iCut]);
948                         fHistoTrueSecondaryConvGammaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGamma_Pt","ESD_TrueSecondaryConvGamma_Pt",250,0,25);
949                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaPt[iCut]);
950       
951                         fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromK0s_Pt", "ESD_TrueSecondaryConvGammaFromXFromK0s_Pt",250,0,25);
952                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromK0sPt[iCut]);
953                         fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromLambda_Pt", "ESD_TrueSecondaryConvGammaFromXFromLambda_Pt",250,0,25);
954                         fTrueList[iCut]->Add(fHistoTrueSecondaryConvGammaFromXFromLambdaPt[iCut]);
955       
956                         fHistoTrueDalitzPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueDalitzPsiPairDeltaPhi_Pt", "ESD_TrueDalitzPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
957                         fTrueList[iCut]->Add(fHistoTrueDalitzPsiPairDeltaPhi[iCut]);
958                         
959                         fHistoTrueGammaPsiPairDeltaPhi[iCut] = new TH2F("ESD_TrueGammaPsiPairDeltaPhi_Pt", "ESD_TrueGammaPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
960                         fTrueList[iCut]->Add(fHistoTrueGammaPsiPairDeltaPhi[iCut]);
961       
962                         fHistoTruePrimaryConvGammaESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryConvGammaESD_PtMCPt", "ESD_TruePrimaryConvGammaESD_PtMCPt",250,0,25,250,0,25);
963                         fTrueList[iCut]->Add(fHistoTruePrimaryConvGammaESDPtMCPt[iCut]);
964                 
965                         fHistoTrueClusGammaPt[iCut] = new TH1F("TrueClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25);
966                         fTagOutputList[iCut]->Add(fHistoTrueClusGammaPt[iCut]);
967                         fHistoTruePrimaryClusGammaPt[iCut] = new TH1F("TruePrimaryClusGamma_Pt","ESD_TruePrimaryClusGamma_Pt",250,0,25);
968                         fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaPt[iCut]);
969                         fHistoTruePrimaryClusGammaESDPtMCPt[iCut] = new TH2F("TruePrimaryClusGamma_Pt_MCPt","ESD_TruePrimaryClusGamma_MCPt",250,0,25,250,0,25);
970                         fTagOutputList[iCut]->Add(fHistoTruePrimaryClusGammaESDPtMCPt[iCut]);
971
972                         if(fDoMesonAnalysis){
973                                 fHistoTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",800,0,0.8,250,0,25);
974                                 fTrueList[iCut]->Add(fHistoTrueMotherInvMassPt[iCut]);
975                                 fHistoTruePrimaryMotherInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_Pt", "ESD_TruePrimaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
976                                 fHistoTruePrimaryMotherInvMassPt[iCut]->Sumw2();
977                                 fTrueList[iCut]->Add(fHistoTruePrimaryMotherInvMassPt[iCut]);
978                                 fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut] = new TH2F("ESD_TruePrimaryMotherW0Weights_InvMass_Pt", "ESD_TruePrimaryMotherW0Weights_InvMass_Pt", 800,0,0.8,250,0,25);
979                                 fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]->Sumw2();
980                                 fTrueList[iCut]->Add(fHistoTruePrimaryMotherW0WeightingInvMassPt[iCut]);
981                                 fProfileTruePrimaryMotherWeightsInvMassPt[iCut] = new TProfile2D("ESD_TruePrimaryMotherWeights_InvMass_Pt", "ESD_TruePrimaryMotherWeights_InvMass_Pt", 800,0,0.8,250,0,25);
982                                 fProfileTruePrimaryMotherWeightsInvMassPt[iCut]->Sumw2();
983                                 fTrueList[iCut]->Add(fProfileTruePrimaryMotherWeightsInvMassPt[iCut]);
984                                 fHistoTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMother_InvMass_Pt", "ESD_TrueSecondaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
985                                 fHistoTrueSecondaryMotherInvMassPt[iCut]->Sumw2();
986                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherInvMassPt[iCut]);
987                                 fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueSecondaryMotherFromK0s_InvMass_Pt",800,0,0.8,250,0,25);
988                                 fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]->Sumw2();
989                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromK0sInvMassPt[iCut]);
990                                 fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromEta_InvMass_Pt","ESD_TrueSecondaryMotherFromEta_InvMass_Pt",800,0,0.8,250,0,25);
991                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromEtaInvMassPt[iCut]);
992                                 fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromLambda_InvMass_Pt","ESD_TrueSecondaryMotherFromLambda_InvMass_Pt",800,0,0.8,250,0,25);
993                                 fTrueList[iCut]->Add(fHistoTrueSecondaryMotherFromLambdaInvMassPt[iCut]);
994                                 fHistoTrueMotherDalitzInvMassPt[iCut] = new TH2F("ESD_TrueDalitz_InvMass_Pt","ESD_TrueDalitz_InvMass_Pt",800,0,0.8,250,0,25);
995                                 fTrueList[iCut]->Add(fHistoTrueMotherDalitzInvMassPt[iCut]);
996                                 if (fDoMesonQA > 0){
997                                         fHistoTruePrimaryPi0MCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryPi0_MCPt_ResolPt","ESD_TruePrimaryPi0_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
998                                         fHistoTruePrimaryPi0MCPtResolPt[iCut]->Sumw2();
999                                         SetLogBinningXTH2(fHistoTruePrimaryPi0MCPtResolPt[iCut]);
1000                                         fTrueList[iCut]->Add(fHistoTruePrimaryPi0MCPtResolPt[iCut]);
1001                                         fHistoTruePrimaryEtaMCPtResolPt[iCut]  = new TH2F("ESD_TruePrimaryEta_MCPt_ResolPt","ESD_TruePrimaryEta_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
1002                                         fHistoTruePrimaryEtaMCPtResolPt[iCut]->Sumw2();
1003                                         SetLogBinningXTH2(fHistoTruePrimaryEtaMCPtResolPt[iCut]);
1004                                         fTrueList[iCut]->Add(fHistoTruePrimaryEtaMCPtResolPt[iCut]);
1005                                         fHistoTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt","ESD_TrueBckGG_InvMass_Pt",800,0,0.8,250,0,25);
1006                                         fTrueList[iCut]->Add(fHistoTrueBckGGInvMassPt[iCut]);
1007                                         fHistoTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt","ESD_TrueBckCont_InvMass_Pt",800,0,0.8,250,0,25);
1008                                         fTrueList[iCut]->Add(fHistoTrueBckContInvMassPt[iCut]);
1009                                         fHistoTrueK0sWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueK0sWithPi0Daughter_MCPt","ESD_TrueK0sWithPi0Daughter_MCPt",250,0,25);
1010                                         fTrueList[iCut]->Add(fHistoTrueK0sWithPi0DaughterMCPt[iCut]);
1011                                         fHistoTrueEtaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueEtaWithPi0Daughter_MCPt","ESD_TrueEtaWithPi0Daughter_MCPt",250,0,25);
1012                                         fTrueList[iCut]->Add(fHistoTrueEtaWithPi0DaughterMCPt[iCut]);
1013                                         fHistoTrueLambdaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueLambdaWithPi0Daughter_MCPt","ESD_TrueLambdaWithPi0Daughter_MCPt",250,0,25);
1014                                         fTrueList[iCut]->Add(fHistoTrueLambdaWithPi0DaughterMCPt[iCut]);
1015           
1016                                         fHistoTruePi0PtY[iCut] = new TH2F("ESD_TruePi0_Pt_Y","ESD_TruePi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
1017                                         SetLogBinningXTH2(fHistoTruePi0PtY[iCut]);
1018                                         fTrueList[iCut]->Add(fHistoTruePi0PtY[iCut]);
1019                                         fHistoTrueEtaPtY[iCut] = new TH2F("ESD_TrueEta_Pt_Y","ESD_TrueEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
1020                                         SetLogBinningXTH2(fHistoTrueEtaPtY[iCut]);
1021                                         fTrueList[iCut]->Add(fHistoTrueEtaPtY[iCut]);
1022                                         fHistoTruePi0PtAlpha[iCut] = new TH2F("ESD_TruePi0_Pt_Alpha","ESD_TruePi0_Pt_Alpha",150,0.03,15.,100,0,1);
1023                                         SetLogBinningXTH2(fHistoTruePi0PtAlpha[iCut]);
1024                                         fTrueList[iCut]->Add(fHistoTruePi0PtAlpha[iCut]);
1025                                         fHistoTrueEtaPtAlpha[iCut] = new TH2F("ESD_TrueEta_Pt_Alpha","ESD_TrueEta_Pt_Alpha",150,0.03,15.,100,0,1);
1026                                         SetLogBinningXTH2(fHistoTrueEtaPtAlpha[iCut]);
1027                                         fTrueList[iCut]->Add(fHistoTrueEtaPtAlpha[iCut]);
1028                                         
1029                                         fHistoTruePi0PtOpenAngle[iCut] = new TH2F("ESD_TruePi0_Pt_OpenAngle","ESD_TruePi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
1030                                         SetLogBinningXTH2(fHistoTruePi0PtOpenAngle[iCut]);
1031                                         fTrueList[iCut]->Add(fHistoTruePi0PtOpenAngle[iCut]);
1032                                         fHistoTrueEtaPtOpenAngle[iCut] = new TH2F("ESD_TrueEta_Pt_OpenAngle","ESD_TrueEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());
1033                                         SetLogBinningXTH2(fHistoTrueEtaPtOpenAngle[iCut]);
1034                                         fTrueList[iCut]->Add(fHistoTrueEtaPtOpenAngle[iCut]);
1035           
1036                                         fHistoTrueConvGammaEta[iCut] = new TH1F("ESD_TrueConvGamma_Eta","ESD_TrueConvGamma_Eta",2000,-2,2);
1037                                         fTrueList[iCut]->Add(fHistoTrueConvGammaEta[iCut]);
1038           
1039                                 }
1040                         }
1041                 }
1042         }  
1043     
1044         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
1045         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
1046   
1047         if(fV0Reader)
1048                 if((AliConversionCuts*)fV0Reader->GetConversionCuts())
1049                         if(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
1050                                 fOutputContainer->Add(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
1051   
1052         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1053                 if(!((AliConversionCuts*)fCutArray->At(iCut))) continue;
1054                 if(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms()){
1055                         fCutFolder[iCut]->Add(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms());
1056                 }
1057                 if(!((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))) continue;
1058                 if(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms()){
1059                         fCutFolder[iCut]->Add(((AliCaloPhotonCuts*)fClusterCutArray->At(iCut))->GetCutHistograms());
1060                 }
1061                 if(fDoMesonAnalysis){
1062                         if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue;
1063                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
1064                                 fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
1065                         }
1066                 }
1067         }
1068         PostData(1, fOutputContainer);
1069 }
1070 //_____________________________________________________________________________
1071 Bool_t AliAnalysisTaskGammaConvCalo::Notify()
1072 {
1073         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1074                 if(!((AliConversionCuts*)fCutArray->At(iCut))->GetDoEtaShift()){
1075                         fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
1076                         continue; // No Eta Shift requested, continue
1077                 }
1078                 if(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
1079                         ((AliConversionCuts*)fCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
1080                         fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
1081                         ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1082                         continue;
1083                 }
1084                 else{
1085                         printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n",
1086                                         (((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber()).Data(),((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift());
1087                         fProfileEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
1088                         ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
1089                 }
1090         }
1091         
1092         return kTRUE;
1093 }
1094 //_____________________________________________________________________________
1095 void AliAnalysisTaskGammaConvCalo::UserExec(Option_t *)
1096 {
1097         //
1098         // Called for each event
1099         //
1100         Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
1101         if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
1102                 for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1103                 fHistoNEvents[iCut]->Fill(eventQuality);
1104                 }
1105                 return;
1106         }
1107         
1108         if(fIsMC) fMCEvent = MCEvent();
1109         if(fMCEvent == NULL) fIsMC = kFALSE;
1110         
1111         fInputEvent = InputEvent();
1112         
1113         if(fIsMC && fInputEvent->IsA()==AliESDEvent::Class()){
1114                 fMCStack = fMCEvent->Stack();
1115                 if(fMCStack == NULL) fIsMC = kFALSE;
1116         }
1117         
1118         fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
1119
1120         // ------------------- BeginEvent ----------------------------
1121         
1122         AliEventplane *EventPlane = fInputEvent->GetEventplane();
1123         if(fIsHeavyIon ==1)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2);
1124         else fEventPlaneAngle=0.0;
1125         
1126         if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
1127                 RelabelAODPhotonCandidates(kTRUE);    // In case of AODMC relabeling MC
1128                 fV0Reader->RelabelAODs(kTRUE);
1129         }
1130         
1131
1132         for(Int_t iCut = 0; iCut<fnCuts; iCut++){
1133                 
1134                 fiCut = iCut;
1135                 Int_t eventNotAccepted = ((AliConversionCuts*)fCutArray->At(iCut))->IsEventAcceptedByConversionCut(fV0Reader->GetConversionCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
1136                 
1137                 if(eventNotAccepted){
1138                 // cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
1139                         fHistoNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
1140                         continue;
1141                 }
1142
1143                 if(eventQuality != 0){// Event Not Accepted
1144                         //cout << "event rejected due to: " <<eventQuality << endl;
1145                         fHistoNEvents[iCut]->Fill(eventQuality);
1146                         continue;
1147                 }
1148
1149                 fHistoNEvents[iCut]->Fill(eventQuality); // Should be 0 here
1150                 fHistoNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks());
1151                 if(((AliConversionCuts*)fCutArray->At(iCut))->IsHeavyIon() == 2)        fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A());
1152                         else fHistoNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C());
1153
1154                 if(fIsMC){
1155                         // Process MC Particle
1156                         if(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection() != 0){
1157                                 if(fInputEvent->IsA()==AliESDEvent::Class()){
1158                                 ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(),
1159                                                                                                                                                                         ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(),
1160                                                                                                                                                                         fMCEvent);
1161                                 }
1162                                 else if(fInputEvent->IsA()==AliAODEvent::Class()){
1163                                 ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(),
1164                                                                                                                                                                         ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(),
1165                                                                                                                                                                         fInputEvent);
1166                                 }
1167
1168                                 if(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader()){
1169                                         for(Int_t i = 0;i<(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){
1170                                                 TString nameBin= fHistoMCHeaders[iCut]->GetXaxis()->GetBinLabel(i+1);
1171                                                 if (nameBin.CompareTo("")== 0){
1172                                                 TString nameHeader = ((TObjString*)((TList*)((AliConversionCuts*)fCutArray->At(iCut))
1173                                                                                                                         ->GetAcceptedHeader())->At(i))->GetString();
1174                                                 fHistoMCHeaders[iCut]->GetXaxis()->SetBinLabel(i+1,nameHeader.Data());
1175                                                 }
1176                                         }
1177                                 }
1178                         }
1179                 }
1180                 if(fIsMC){
1181                 if(fInputEvent->IsA()==AliESDEvent::Class())
1182                         ProcessMCParticles();
1183                 if(fInputEvent->IsA()==AliAODEvent::Class())
1184                         ProcessAODMCParticles();
1185                 }
1186                 
1187                 // it is in the loop to have the same conversion cut string (used also for MC stuff that should be same for V0 and Cluster)
1188                 ProcessClusters();                                      // process calo clusters
1189                 ProcessPhotonCandidates();                              // Process this cuts gammas
1190
1191                 fHistoNGammaCandidates[iCut]->Fill(fGammaCandidates->GetEntries());
1192                 fHistoNGoodESDTracksVsNGammaCanditates[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks(),fGammaCandidates->GetEntries());
1193                 if(fDoMesonAnalysis){ // Meson Analysis
1194                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1195                         fUnsmearedPx = new Double_t[fGammaCandidates->GetEntries()]; // Store unsmeared Momenta
1196                         fUnsmearedPy = new Double_t[fGammaCandidates->GetEntries()];
1197                         fUnsmearedPz = new Double_t[fGammaCandidates->GetEntries()];
1198                         fUnsmearedE =  new Double_t[fGammaCandidates->GetEntries()];
1199                         
1200                         for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1201                         fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Px();
1202                         fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Py();
1203                         fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Pz();
1204                         fUnsmearedE[gamma] =  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->E();
1205                         ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(gamma)));
1206                         }
1207                 }
1208
1209                 PhotonTagging(); // tag PCM photons with calorimeter
1210
1211                 CalculatePi0Candidates(); // Combine Gammas from conversion and from calo
1212
1213                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
1214                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
1215                         CalculateBackground(); // Combinatorial Background
1216                         UpdateEventByEventData(); // Store Event for mixed Events
1217                         }
1218                         else{
1219                         CalculateBackgroundRP(); // Combinatorial Background
1220                         fBGHandlerRP[iCut]->AddEvent(fGammaCandidates,fInputEvent); // Store Event for mixed Events
1221                         fBGClusHandlerRP[iCut]->AddEvent(fClusterCandidates,fInputEvent); // Store Event for mixed Events
1222                         }
1223                 }
1224                 if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1225                         for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1226                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
1227                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPy(fUnsmearedPy[gamma]);
1228                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPz(fUnsmearedPz[gamma]);
1229                         ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetE(fUnsmearedE[gamma]);
1230                         }
1231                         delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
1232                         delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
1233                         delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
1234                         delete[] fUnsmearedE;  fUnsmearedE  = 0x0;
1235                 }
1236                 }
1237
1238                 fGammaCandidates->Clear(); // delete this cuts good gammas
1239                 fClusterCandidates->Clear(); // delete cluster candidates
1240         }
1241         
1242         if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
1243                 RelabelAODPhotonCandidates(kFALSE); // Back to ESDMC Label
1244                 fV0Reader->RelabelAODs(kFALSE);
1245         }
1246         
1247         PostData(1, fOutputContainer);
1248 }
1249
1250 //________________________________________________________________________
1251 void AliAnalysisTaskGammaConvCalo::ProcessClusters()
1252 {
1253         
1254         Int_t nclus = 0;
1255         nclus = fInputEvent->GetNumberOfCaloClusters();
1256         
1257         cout << nclus << endl;
1258         
1259         if(nclus == 0)  return;
1260         
1261         // vertex
1262         Double_t vertex[3] = {0};
1263         InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
1264         
1265         // Loop over EMCal clusters
1266         for(Int_t i = 0; i < nclus; i++){
1267                 
1268                 AliVCluster* clus = fInputEvent->GetCaloCluster(i);             
1269                 if (!clus) continue;
1270                 if(!((AliCaloPhotonCuts*)fClusterCutArray->At(fiCut))->ClusterIsSelected(clus,fInputEvent,fIsMC)) continue;
1271                 // TLorentzvector with cluster
1272                 TLorentzVector clusterVector;
1273                 clus->GetMomentum(clusterVector,vertex);
1274                 
1275                 TLorentzVector* tmpvec = new TLorentzVector();
1276                 tmpvec->SetPxPyPzE(clusterVector.Px(),clusterVector.Py(),clusterVector.Pz(),clusterVector.E());
1277                 
1278                 // convert to AODConversionPhoton
1279                 AliAODConversionPhoton *PhotonCandidate=new AliAODConversionPhoton(tmpvec);
1280                 if(!PhotonCandidate) continue;
1281
1282                 // get MC label
1283                 Int_t mclab[2] = {0,0};         
1284                 if(fIsMC){
1285                         mclab[0] = clus->GetLabel();
1286                         mclab[1] = clus->GetLabel();
1287                         cout << "mclabels: " << mclab[0] << " " << mclab[1] << endl;
1288                 }
1289                 
1290                 PhotonCandidate->SetMCLabel(mclab);
1291                 
1292         //    fIsFromMBHeader = kTRUE;
1293         //    if(fIsMC){
1294         //      Int_t isPosFromMBHeader
1295         //      = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent);
1296         //      if(isPosFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1297         // 
1298         //      Int_t isNegFromMBHeader
1299         //      = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(mclab, fMCStack, fInputEvent);
1300         //      if(isNegFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1301         //      
1302         //      if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1303         //    }
1304                 
1305                 fHistoClusGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1306                 fClusterCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1307                 
1308         //    if(fIsFromMBHeader){
1309         //      fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1310         //      if (fDoPhotonQA > 0){
1311         //        fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1312         //        fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1313         //      }
1314         //    }
1315
1316                 if(fIsMC){
1317                         ProcessTrueClusterCandidates(PhotonCandidate);
1318                 }
1319                 
1320         //    if (fIsFromMBHeader && fDoPhotonQA == 2){
1321         //      if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1322         //        fPtGamma = PhotonCandidate->Pt();
1323         //        fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1324         //        fRConvPhoton = PhotonCandidate->GetConversionRadius();
1325         //        fEtaPhoton = PhotonCandidate->GetPhotonEta();
1326         //        fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1327         //        fTreeConvGammaPtDcazCat[fiCut]->Fill();
1328         //    } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1329         //        fPtGamma = PhotonCandidate->Pt();
1330         //        fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1331         //        fRConvPhoton = PhotonCandidate->GetConversionRadius();
1332         //        fEtaPhoton = PhotonCandidate->GetPhotonEta();
1333         //        fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1334         //        fTreeConvGammaPtDcazCat[fiCut]->Fill();
1335         //      }
1336         //    }
1337                 
1338         }
1339         
1340 }
1341
1342 //________________________________________________________________________
1343 void AliAnalysisTaskGammaConvCalo::ProcessTrueClusterCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1344 {
1345         
1346         TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack);
1347         
1348         if(Photon == NULL){
1349         //    cout << "no photon" << endl;
1350                 return;
1351         }
1352         
1353         if(Photon->GetPdgCode() != 22){
1354         //    cout << "pdg code: " << Photon->GetPdgCode() << endl;
1355                 return; // Particle is no Photon
1356         }
1357         
1358         // True Photon
1359         if(fIsFromMBHeader){
1360                 fHistoTrueClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1361                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1362         }
1363
1364         if(Photon->GetMother(0) <= fMCStack->GetNprimary()){
1365                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1366                 if(fIsFromMBHeader){
1367                 fHistoTruePrimaryClusGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1368                 fHistoTruePrimaryClusGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1369                 }
1370         }
1371
1372         // maybe later
1373         //  else{
1374         //    if(fIsFromMBHeader){
1375         //      fCharPhotonMCInfo = 2;
1376         //      fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1377         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1378         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){
1379         //        fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1380         //        fCharPhotonMCInfo = 5;
1381         //      }
1382         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1383         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
1384         //        fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1385         //        fCharPhotonMCInfo = 4;
1386         //      }
1387         //      if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1388         //         fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){
1389         //        fCharPhotonMCInfo = 3;
1390         //      }
1391         //    }
1392         //  }
1393 }
1394
1395
1396 //________________________________________________________________________
1397 void AliAnalysisTaskGammaConvCalo::ProcessPhotonCandidates()
1398 {
1399         Int_t nV0 = 0;
1400         TList *GammaCandidatesStepOne = new TList();
1401         TList *GammaCandidatesStepTwo = new TList();
1402         // Loop over Photon Candidates allocated by ReaderV1
1403         for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
1404                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
1405                 if(!PhotonCandidate) continue;
1406                 fIsFromMBHeader = kTRUE;
1407                 if(fIsMC && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1408                         Int_t isPosFromMBHeader
1409                         = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1410                         if(isPosFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1411                         Int_t isNegFromMBHeader
1412                         = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1413                         if(isNegFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1414                         
1415                         if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1416                 }
1417                 
1418                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue;
1419                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue;
1420                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1421                 !((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1422                         fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1423                         
1424                         if(fIsFromMBHeader){
1425                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1426                                 if (fDoPhotonQA > 0){
1427                                 fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1428                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1429                                 }
1430                         }
1431                         if(fIsMC){
1432                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1433                                 ProcessTruePhotonCandidates(PhotonCandidate);
1434                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1435                                 ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1436                         }
1437                         if (fIsFromMBHeader && fDoPhotonQA == 2){
1438                                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1439                                 fPtGamma = PhotonCandidate->Pt();
1440                                 fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1441                                 fRConvPhoton = PhotonCandidate->GetConversionRadius();
1442                                 fEtaPhoton = PhotonCandidate->GetPhotonEta();
1443                                 fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1444                                 fTreeConvGammaPtDcazCat[fiCut]->Fill();
1445                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1446                                 fPtGamma = PhotonCandidate->Pt();
1447                                 fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1448                                 fRConvPhoton = PhotonCandidate->GetConversionRadius();
1449                                 fEtaPhoton = PhotonCandidate->GetPhotonEta();
1450                                 fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1451                                 fTreeConvGammaPtDcazCat[fiCut]->Fill();
1452                                 }
1453                         }
1454                 } else if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
1455                         ((AliConversionCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
1456                         nV0++;
1457                         GammaCandidatesStepOne->Add(PhotonCandidate);
1458                 } else if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1459                                 ((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
1460                         GammaCandidatesStepTwo->Add(PhotonCandidate);
1461                 }
1462         }
1463         if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
1464                 for(Int_t i = 0;i<GammaCandidatesStepOne->GetEntries();i++){
1465                         AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i);
1466                         if(!PhotonCandidate) continue;
1467                         fIsFromMBHeader = kTRUE;
1468                         if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1469                                 Int_t isPosFromMBHeader
1470                                 = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1471                                 Int_t isNegFromMBHeader
1472                                 = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1473                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1474                         }
1475                         if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue;
1476                         if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
1477                                 fGammaCandidates->Add(PhotonCandidate);
1478                                 if(fIsFromMBHeader){
1479                                         fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1480                                         if (fDoPhotonQA > 0){
1481                                                 fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1482                                                 fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1483                                         }
1484                                 }
1485                         }
1486                         if(fIsMC){
1487                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1488                                         ProcessTruePhotonCandidates(PhotonCandidate);
1489                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1490                                         ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1491                         } else GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
1492                         
1493                         if (fIsFromMBHeader && fDoPhotonQA == 2){
1494                                 if (fIsHeavyIon ==1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1495                                         fPtGamma = PhotonCandidate->Pt();
1496                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1497                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1498                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1499                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1500                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1501                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1502                                         fPtGamma = PhotonCandidate->Pt();
1503                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1504                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1505                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1506                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1507                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1508                                 }
1509                         }
1510                 }
1511         }
1512         if(((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1513                 for(Int_t i = 0;i<GammaCandidatesStepTwo->GetEntries();i++){
1514                         AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i);
1515                         if(!PhotonCandidate) continue;
1516                         fIsFromMBHeader = kTRUE;
1517                         if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1518                                 Int_t isPosFromMBHeader
1519                                 = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1520                                 Int_t isNegFromMBHeader
1521                                 = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1522                                 if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1523                         }
1524                         if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue;
1525                         fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList
1526                         if(fIsFromMBHeader){
1527                                 fHistoConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1528                                 if (fDoPhotonQA > 0){
1529                                         fHistoConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1530                                         fHistoConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1531                                 }
1532                         }
1533                         if(fIsMC){
1534                                 if(fInputEvent->IsA()==AliESDEvent::Class())
1535                                         ProcessTruePhotonCandidates(PhotonCandidate);
1536                                 if(fInputEvent->IsA()==AliAODEvent::Class())
1537                                         ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1538                         }
1539                         if (fIsFromMBHeader){
1540                                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1541                                         fPtGamma = PhotonCandidate->Pt();
1542                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1543                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1544                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1545                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1546                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1547                                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1548                                         fPtGamma = PhotonCandidate->Pt();
1549                                         fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1550                                         fRConvPhoton = PhotonCandidate->GetConversionRadius();
1551                                         fEtaPhoton = PhotonCandidate->GetPhotonEta();
1552                                         fCharCatPhoton = PhotonCandidate->GetPhotonQuality();
1553                                         fTreeConvGammaPtDcazCat[fiCut]->Fill();
1554                                 }
1555                         }
1556                 }
1557         }
1558         
1559         delete GammaCandidatesStepOne;
1560         GammaCandidatesStepOne = 0x0;
1561         delete GammaCandidatesStepTwo;
1562         GammaCandidatesStepTwo = 0x0;
1563   
1564 }
1565 //________________________________________________________________________
1566 void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate)
1567 {
1568         
1569         Double_t magField = fInputEvent->GetMagneticField();
1570         if( magField  < 0.0 ){
1571                 magField =  1.0;
1572         }
1573         else {
1574                 magField =  -1.0;
1575         }
1576         
1577         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1578         AliAODMCParticle *posDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelPositive());
1579         AliAODMCParticle *negDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelNegative());
1580         fCharPhotonMCInfo = 0;
1581         
1582         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1583         Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1584         
1585         if(posDaughter->GetMother() != negDaughter->GetMother()){
1586                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1587                 fCharPhotonMCInfo = 1;
1588                 return;
1589         }
1590         else if(posDaughter->GetMother() == -1){
1591                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1592                 fCharPhotonMCInfo = 1;
1593                 return;
1594         }
1595         
1596         if(pdgCode[0]!=11 || pdgCode[1]!=11){
1597                 fCharPhotonMCInfo = 1;
1598                 return; //One Particle is not a electron
1599         }
1600         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()){
1601                 fCharPhotonMCInfo = 1;
1602                 return; // Same Charge
1603         }
1604         
1605         AliAODMCParticle *Photon = (AliAODMCParticle*) AODMCTrackArray->At(posDaughter->GetMother());
1606         AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1607         AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1608         Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1609         
1610         if(Photon->GetPdgCode() != 22){
1611                 fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1612                 fCharPhotonMCInfo = 1;
1613                 return; // Mother is no Photon
1614         }
1615         
1616         if(((posDaughter->GetMCProcessCode())) != 5 || ((negDaughter->GetMCProcessCode())) != 5){
1617                 fCharPhotonMCInfo = 1;
1618                 return;// check if the daughters come from a conversion
1619         }
1620         // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
1621         
1622         
1623         
1624         // True Photon
1625         if(fIsFromMBHeader){
1626                 fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1627                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1628         }
1629         fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1630         if(Photon->IsPrimary()){
1631                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1632                 if(fIsFromMBHeader){
1633                         fCharPhotonMCInfo = 6;
1634                         fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1635                         fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1636                 }
1637                 // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1638         } else {
1639                 if(fIsFromMBHeader){
1640                         fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1641                         fCharPhotonMCInfo = 2;
1642                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1643                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 3122){
1644                                 fCharPhotonMCInfo = 5;
1645                                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1646                         }
1647                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1648                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 310){
1649                                 fCharPhotonMCInfo = 4;
1650                                 fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1651                         }
1652                         if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1653                                 ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 221){
1654                                 fCharPhotonMCInfo = 3;
1655                         }
1656                 }
1657         }
1658   
1659 }
1660 //________________________________________________________________________
1661 void AliAnalysisTaskGammaConvCalo::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1662 {
1663         
1664         Double_t magField = fInputEvent->GetMagneticField();
1665         if( magField  < 0.0 ){
1666                 magField =  1.0;
1667         }
1668         else {
1669                 magField =  -1.0;
1670         }
1671         
1672         // Process True Photons
1673         TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(fMCStack);
1674         TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(fMCStack);  
1675         fCharPhotonMCInfo = 0;
1676         
1677         if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1678         Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1679         fCharPhotonMCInfo = 1;
1680         if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){
1681                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1682                 return;
1683         }
1684         else if(posDaughter->GetMother(0) == -1){
1685                 FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1686                 return;
1687         }
1688         
1689         if(pdgCode[0]!=11 || pdgCode[1]!=11) return; //One Particle is not a electron
1690         
1691         if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
1692         
1693         TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack);
1694         AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1695         AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1696         Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1697         
1698         if(Photon->GetPdgCode() != 22){
1699                 fHistoTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1700                 return; // Mother is no Photon
1701         }
1702         
1703         if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
1704         
1705         
1706         
1707         // True Photon
1708         if(fIsFromMBHeader){
1709                 fHistoTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1710                 if (fDoPhotonQA > 0) fHistoTrueConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
1711         }
1712         fHistoTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1713         if(posDaughter->GetMother(0) <= fMCStack->GetNprimary()){
1714                 // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1715                 if(fIsFromMBHeader){
1716                         fCharPhotonMCInfo = 6;
1717                         fHistoTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1718                         fHistoTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1719                 
1720                 }
1721                 // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1722         }
1723         else{
1724                 if(fIsFromMBHeader){
1725                         fCharPhotonMCInfo = 2;
1726                         fHistoTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1727                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1728                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){
1729                                 fHistoTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1730                                 fCharPhotonMCInfo = 5;
1731                         }
1732                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1733                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
1734                                 fHistoTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1735                                 fCharPhotonMCInfo = 4;
1736                         }
1737                         if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1738                                 fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){
1739                                 fCharPhotonMCInfo = 3;
1740                         }
1741                 }
1742         }
1743         
1744         // pi0 photon
1745         //Bool_t bpi0 = 0;
1746         Int_t imother = Photon->GetMother(0);
1747         AliMCParticle *McMother = static_cast<AliMCParticle*>(fMCEvent->GetTrack(imother));
1748         if(McMother->PdgCode() == 111){
1749                 fHistoTrueConvPi0GammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1750         }
1751   
1752 }
1753 //________________________________________________________________________
1754 void AliAnalysisTaskGammaConvCalo::ProcessAODMCParticles()
1755 {
1756         
1757         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1758         
1759         // Loop over all primary MC particle
1760         for(Int_t i = 0; i < AODMCTrackArray->GetEntriesFast(); i++) {
1761                 
1762                 AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(i));
1763                 if (!particle) continue;
1764                 if (!particle->IsPrimary()) continue;
1765                 
1766                 Int_t isMCFromMBHeader = -1;
1767                 if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1768                         isMCFromMBHeader
1769                         = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1770                         if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1771                 }
1772                 
1773                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1774                 if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kFALSE)){
1775                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1776                                 if(particle->GetMother() >-1){ // Meson Decay Gamma
1777                                         switch((static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetMother())))->GetPdgCode()){
1778                                         case 111: // Pi0
1779                                                 fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1780                                                 break;
1781                                         case 113: // Rho0
1782                                                 fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1783                                                 break;
1784                                         case 221: // Eta
1785                                                 fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1786                                                 break;
1787                                         case 223: // Omega
1788                                                 fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1789                                                 break;
1790                                         case 331: // Eta'
1791                                                 fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1792                                                 break;
1793                                         case 333: // Phi
1794                                                 fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1795                                                 break;
1796                                         case 3212: // Sigma
1797                                                 fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1798                                                 break;
1799                                         }
1800                                 }
1801                         }
1802                         if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kTRUE)){
1803                                 Double_t rConv = 0;
1804                                 for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
1805                                         AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(daughterIndex));
1806                                         if(!tmpDaughter) continue;
1807                                         if(abs(tmpDaughter->GetPdgCode()) == 11){
1808                                                 rConv = sqrt( (tmpDaughter->Xv()*tmpDaughter->Xv()) + (tmpDaughter->Yv()*tmpDaughter->Yv()) );
1809                                         }
1810                                 }
1811                                 fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1812                                 if (fDoPhotonQA > 0){
1813                                         fHistoMCConvGammaR[fiCut]->Fill(rConv);
1814                                         fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta());
1815                                 }
1816                         }
1817                         // Converted MC Gamma
1818                         if(fDoMesonAnalysis){
1819                         if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1820                                 Double_t mesonY = 10.;
1821                                 if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1822                                         mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1823                                 } else {
1824                                         mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1825                                 }
1826                                 Float_t weightedK0s= 1;
1827                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1828                                         if (particle->Pt()>0.005){
1829                                                 weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1830                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1831                                         }
1832                                 }
1833                                 fHistoMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1834                                 fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1835                                 fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1836                         }
1837                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1838                                 ->MesonIsSelectedAODMC(particle,AODMCTrackArray,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
1839                                 AliAODMCParticle* daughter0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(0)));
1840                                 AliAODMCParticle* daughter1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(1)));
1841                                 Float_t weighted= 1;
1842                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1843                                         if (particle->Pt()>0.005){
1844                                                 weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1845                                                 //                   if(particle->GetPdgCode() == 221){
1846                                                 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1847                                                 //                   }
1848                                         }
1849                                 }
1850                                 Double_t mesonY = 10.;
1851                                 if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1852                                         mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1853                                 } else{
1854                                         mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1855                                 }
1856                                 
1857                                 if(particle->GetPdgCode() == 111){
1858                                         fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1859                                         fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1860                                         if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1861                                 } else if(particle->GetPdgCode() == 221){
1862                                         fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1863                                         fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1864                                         if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1865                                 }
1866                                 
1867                                 // Check the acceptance for both gammas
1868                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) &&
1869                                 ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE)  &&
1870                                 ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1871                                 ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1872                                         
1873                                         if(particle->GetPdgCode() == 111){
1874                                                 fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
1875                                         } else if(particle->GetPdgCode() == 221){
1876                                                 fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
1877                                         }
1878                                 }
1879                         }
1880                 }
1881         }
1882         
1883 }
1884 //________________________________________________________________________
1885 void AliAnalysisTaskGammaConvCalo::ProcessMCParticles()
1886 {
1887         // Loop over all primary MC particle
1888         for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1889                 TParticle* particle = (TParticle *)fMCStack->Particle(i);
1890                 if (!particle) continue;
1891                 
1892                 Int_t isMCFromMBHeader = -1;
1893                 if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1894                         isMCFromMBHeader
1895                         = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1896                         if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1897                 }
1898                 
1899                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1900                 if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1901                         fHistoMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1902                         if(particle->GetMother(0) >-1){ // Meson Decay Gamma
1903                                 switch(fMCStack->Particle(particle->GetMother(0))->GetPdgCode()){
1904                                 case 111: // Pi0
1905                                         fHistoMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1906                                         break;
1907                                 case 113: // Rho0
1908                                         fHistoMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1909                                         break;
1910                                 case 221: // Eta
1911                                         fHistoMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1912                                         break;
1913                                 case 223: // Omega
1914                                         fHistoMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1915                                         break;
1916                                 case 331: // Eta'
1917                                         fHistoMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1918                                         break;
1919                                 case 333: // Phi
1920                                         fHistoMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1921                                         break;
1922                                 case 3212: // Sigma
1923                                         fHistoMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1924                                         break;
1925                                 }
1926                         }
1927                 }
1928                 if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1929                         fHistoMCConvGammaPt[fiCut]->Fill(particle->Pt());
1930                         if (fDoPhotonQA > 0){
1931                                 fHistoMCConvGammaR[fiCut]->Fill(((TParticle*)fMCStack->Particle(particle->GetFirstDaughter()))->R());
1932                                 fHistoMCConvGammaEta[fiCut]->Fill(particle->Eta());
1933                         }
1934                 } // Converted MC Gamma
1935                 if(fDoMesonAnalysis){
1936                         if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1937                                 Double_t mesonY = 10.;
1938                                 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1939                                         mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1940                                 } else{
1941                                         mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1942                                 }
1943                                 Float_t weightedK0s= 1;
1944                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1945                                         if (particle->Pt()>0.005){
1946                                                 weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1947                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1948                                         }
1949                                 }
1950                                 if (fMCStack->IsPhysicalPrimary(i)){
1951                                         fHistoMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1952                                         fHistoMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1953                                         fHistoMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1954                                 }
1955                         }
1956                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1957                                 ->MesonIsSelectedMC(particle,fMCStack,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
1958                                 TParticle* daughter0 = (TParticle*)fMCStack->Particle(particle->GetFirstDaughter());
1959                                 TParticle* daughter1 = (TParticle*)fMCStack->Particle(particle->GetLastDaughter());
1960                                 
1961                                 Float_t weighted= 1;
1962                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1963                                         if (particle->Pt()>0.005){
1964                                                 weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1965                                                 //                   if(particle->GetPdgCode() == 221){
1966                                                 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1967                                                 //                   }
1968                                         }
1969                                 }
1970                                 Double_t mesonY = 10.;
1971                                 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1972                                         mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1973                                 } else{
1974                                         mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1975                                 }
1976                                 
1977                                 if(particle->GetPdgCode() == 111){
1978                                         fHistoMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1979                                         fHistoMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1980                                         if (fDoMesonQA > 0) fHistoMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1981                                 } else if(particle->GetPdgCode() == 221){
1982                                         fHistoMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1983                                         fHistoMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1984                                         if (fDoMesonQA > 0) fHistoMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1985                                 }
1986                                 
1987                                 // Check the acceptance for both gammas
1988                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCStack,kFALSE) &&
1989                                 ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCStack,kFALSE)  &&
1990                                 ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1991                                 ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1992                                         
1993                                         if(particle->GetPdgCode() == 111){
1994                                                 fHistoMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
1995                                         } else if(particle->GetPdgCode() == 221){
1996                                                 fHistoMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
1997                                         }
1998                                 }
1999                         }
2000                 }
2001         }
2002   
2003         if (fDoMesonQA){
2004                 for(Int_t i = fMCStack->GetNprimary(); i < fMCStack->GetNtrack(); i++) {
2005                         TParticle* particle = (TParticle *)fMCStack->Particle(i);
2006                         if (!particle) continue;
2007       
2008                         Int_t isMCFromMBHeader = -1;
2009                         if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
2010                                 isMCFromMBHeader
2011         = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
2012                                 if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
2013                         }
2014       
2015                         if(fDoMesonAnalysis){
2016                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2017            ->MesonIsSelectedMC(particle,fMCStack,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
2018                                         Float_t weighted= 1;
2019                                         if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
2020                                                 if (particle->Pt()>0.005){
2021                                                         weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
2022               //                   if(particle->GetPdgCode() == 221){
2023               //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2024               //                   }
2025                                                 }
2026                                         }
2027                                         
2028                                         if(particle->GetPdgCode() == 111){
2029                                                 Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode();
2030                                                 Int_t source = GetSourceClassification(111,pdgCode);
2031                                                 fHistoMCSecPi0PtvsSource[fiCut]->Fill(particle->Pt(),source,weighted); // All MC Pi0
2032                                                 fHistoMCSecPi0Source[fiCut]->Fill(pdgCode);
2033                                         } else if(particle->GetPdgCode() == 221){
2034                                                 Int_t pdgCode = ((TParticle*)fMCStack->Particle( particle->GetFirstMother() ))->GetPdgCode();
2035                                                 fHistoMCSecEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
2036                                                 fHistoMCSecEtaSource[fiCut]->Fill(pdgCode);
2037                                         }
2038                                 }
2039                         }
2040                 }
2041         }
2042 }
2043
2044 //________________________________________________________________________
2045 void AliAnalysisTaskGammaConvCalo::PhotonTagging(){
2046         
2047         // Conversion Gammas
2048         if(fGammaCandidates->GetEntries()>0){
2049
2050                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries()-1;firstGammaIndex++){
2051                 
2052                         // get conversion photon
2053                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2054                         if (gamma0==NULL) continue;
2055                         
2056                         TLorentzVector photonVector;
2057                         photonVector.SetPxPyPzE(gamma0->GetPx(),gamma0->GetPy(),gamma0->GetPz(),gamma0->GetPhotonP());
2058                         
2059                         Bool_t btagpi0 = 0;
2060                         Bool_t btageta = 0;
2061                         
2062                         // loop over clusters
2063                         for(Int_t secondGammaIndex = 0; secondGammaIndex<fClusterCandidates->GetEntries(); ++secondGammaIndex) {
2064                                 
2065                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2066                                 if (gamma1==NULL) continue;
2067                                 
2068                                 TLorentzVector clusterVector;
2069                                 clusterVector.SetPxPyPzE(gamma1->GetPx(),gamma1->GetPy(),gamma1->GetPz(),gamma1->GetPhotonP());
2070                                 
2071                                 // do the tagging
2072                                 TLorentzVector pairVector = photonVector+clusterVector;
2073                                 
2074                                 // see if pi0?
2075                                 if((pairVector.M() > 0.11 && pairVector.M() < 0.15)){
2076                                 btagpi0 = 1;
2077                                 }
2078                                 // or eta
2079                                 if((pairVector.M() > 0.50 && pairVector.M() < 0.6)){
2080                                 btageta = 1;
2081                                 }
2082                         }// end loop over clusters
2083                         
2084                         if(btagpi0 && btageta)
2085                                 fHistoConvGammaTagged[fiCut]->Fill(photonVector.Pt());
2086                         else if(btagpi0 && !btageta)
2087                                 fHistoConvGammaPi0Tagged[fiCut]->Fill(photonVector.Pt());
2088                         else if(btageta && !btagpi0)
2089                                 fHistoConvGammaEtaTagged[fiCut]->Fill(photonVector.Pt());
2090                         else
2091                                 fHistoConvGammaUntagged[fiCut]->Fill(photonVector.Pt());
2092
2093                 }// end loop over gammas
2094         }// end if
2095         return;
2096 }
2097
2098 //________________________________________________________________________
2099 void AliAnalysisTaskGammaConvCalo::CalculatePi0Candidates(){
2100         
2101         // Conversion Gammas
2102         if(fGammaCandidates->GetEntries()>0){
2103
2104                 // vertex
2105                 Double_t vertex[3] = {0};
2106                 InputEvent()->GetPrimaryVertex()->GetXYZ(vertex);
2107
2108                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2109                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2110                         if (gamma0==NULL) continue;
2111                         
2112                         for(Int_t secondGammaIndex=0;secondGammaIndex<fClusterCandidates->GetEntries();secondGammaIndex++){
2113
2114                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fClusterCandidates->At(secondGammaIndex));
2115                                 if (gamma1==NULL) continue;
2116                                 
2117                                 AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
2118                                 pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
2119                                 
2120                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2121                                         fHistoMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
2122                                         
2123                                         // fill new histograms
2124                                         fHistoPhotonPairAll[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
2125                                         fHistoPhotonPairAllGam[fiCut]->Fill(pi0cand->M(),gamma0->Pt());
2126                                         
2127                                         if(pi0cand->GetAlpha()<0.1)
2128                                                 fHistoMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
2129                                         
2130                                         if (fDoMesonQA > 0){
2131                                                 if ( pi0cand->M() > 0.05 && pi0cand->M() < 0.17){
2132                                                         fHistoMotherPi0PtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2133                                                         fHistoMotherPi0PtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());
2134                                                         fHistoMotherPi0PtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle());      
2135                                                 }
2136                                                 if ( pi0cand->M() > 0.45 && pi0cand->M() < 0.65){
2137                                                         fHistoMotherEtaPtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2138                                                         fHistoMotherEtaPtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());
2139                                                         fHistoMotherEtaPtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle());
2140                                                 }
2141                                         }
2142                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
2143                                                 Int_t zbin = 0;
2144                                                 Int_t mbin = 0;
2145                                                 
2146                                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
2147                                                         zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2148                                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2149                                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2150                                                         } else {
2151                                                                 mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2152                                                         }
2153                                                 } else{
2154                                                         zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2155                                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2156                                                                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2157                                                         } else {
2158                                                                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2159                                                         }
2160                                                 }
2161                                                 Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),(Double_t)zbin,(Double_t)mbin};
2162                                                 fSparseMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
2163                                         }
2164                                         
2165                                         
2166                                         if(fIsMC){
2167                                                 if(fInputEvent->IsA()==AliESDEvent::Class())
2168                                                         ProcessTrueMesonCandidates(pi0cand,gamma0,gamma1);
2169                                                 if(fInputEvent->IsA()==AliAODEvent::Class())
2170                                                         ProcessTrueMesonCandidatesAOD(pi0cand,gamma0,gamma1);
2171                                         }
2172                                         if (fDoMesonQA == 2){
2173                                                 fInvMass = pi0cand->M();
2174                                                 fPt  = pi0cand->Pt();
2175                                                 if (abs(gamma0->GetDCAzToPrimVtx()) < abs(gamma1->GetDCAzToPrimVtx())){
2176                                                         fDCAzGammaMin = gamma0->GetDCAzToPrimVtx();
2177                                                         fDCAzGammaMax = gamma1->GetDCAzToPrimVtx();
2178                                                 } else {
2179                                                         fDCAzGammaMin = gamma1->GetDCAzToPrimVtx();
2180                                                         fDCAzGammaMax = gamma0->GetDCAzToPrimVtx();
2181                                                 }
2182                                                 fCharFlag = pi0cand->GetMesonQuality();
2183                                                 //                   cout << "gamma 0: " << gamma0->GetV0Index()<< "\t" << gamma0->GetPx() << "\t" << gamma0->GetPy() << "\t" <<  gamma0->GetPz() << "\t" << endl;
2184                                                 //                   cout << "gamma 1: " << gamma1->GetV0Index()<< "\t"<< gamma1->GetPx() << "\t" << gamma1->GetPy() << "\t" <<  gamma1->GetPz() << "\t" << endl;
2185                                                 //                    cout << "pi0: "<<fInvMass << "\t" << fPt <<"\t" << fDCAzGammaMin << "\t" << fDCAzGammaMax << "\t" << (Int_t)fCharFlag << "\t" << (Int_t)fCharMesonMCInfo <<endl;
2186                                                 if (fIsHeavyIon == 1 && fPt > 0.399 && fPt < 20. ) {
2187                                                 if (fInvMass > 0.08 && fInvMass < 0.2) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2188                                                 if ((fInvMass > 0.45 && fInvMass < 0.6) &&  (fPt > 0.999 && fPt < 20.) )fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2189                                                 } else if (fPt > 0.299 && fPt < 20. )  {
2190                                                 if ( (fInvMass > 0.08 && fInvMass < 0.2) || (fInvMass > 0.45 && fInvMass < 0.6)) fTreeMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
2191                                                 }
2192                                         }
2193                                 }
2194                                 delete pi0cand;
2195                                 pi0cand=0x0;
2196                         }
2197                 }
2198         }
2199 }
2200 //______________________________________________________________________
2201 void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
2202 {
2203         // Process True Mesons
2204         AliStack *MCStack = fMCEvent->Stack();
2205         fCharMesonMCInfo = 0;
2206         if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
2207                 Bool_t isTruePi0 = kFALSE;
2208                 Bool_t isTrueEta = kFALSE;
2209                 Bool_t isTruePi0Dalitz = kFALSE;
2210                 Bool_t isTrueEtaDalitz = kFALSE;
2211                 Bool_t gamma0DalitzCand = kFALSE;
2212                 Bool_t gamma1DalitzCand = kFALSE;
2213                 Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
2214                 Int_t gamma0MotherLabel = -1;
2215                 if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2216                 // Daughters Gamma 0
2217                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
2218                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
2219                 TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
2220                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2221                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
2222                                 if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
2223                                         gamma0MotherLabel=gammaMC0->GetFirstMother();
2224                                 }
2225                         }
2226                         if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
2227                                 gamma0DalitzCand = kTRUE;
2228                                 gamma0MotherLabel=-111;
2229                         }
2230                         if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
2231                                 gamma0DalitzCand = kTRUE;
2232                                 gamma0MotherLabel=-221;
2233                         }
2234                 }
2235                 }
2236                 if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
2237                         Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
2238                         Int_t gamma1MotherLabel = -1;
2239                         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2240                                 // Daughters Gamma 1
2241                                 TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
2242                                 TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
2243                                 TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
2244                                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2245                                         if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
2246                                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
2247                                                 gamma1MotherLabel=gammaMC1->GetFirstMother();
2248                                                 }
2249                                         }
2250                                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
2251                                                 gamma1DalitzCand = kTRUE;
2252                                                 gamma1MotherLabel=-111;
2253                                         }
2254                                         if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
2255                                                 gamma1DalitzCand = kTRUE;
2256                                                 gamma1MotherLabel=-221;
2257                                         }
2258                                 }
2259                         }
2260                         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2261                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
2262                                         isTruePi0=kTRUE;
2263                                 }
2264                                 if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 221){
2265                                         isTrueEta=kTRUE;
2266                                 }
2267                         }
2268                         
2269                         //Identify Dalitz candidate
2270                         if (gamma1DalitzCand || gamma0DalitzCand){
2271                                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
2272                                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2273                                         if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2274                                 }
2275                                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
2276                                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2277                                         if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2278                                 }
2279                         }
2280                         
2281                         
2282                         if(isTruePi0 || isTrueEta){// True Pion or Eta
2283                                 fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2284                                 if (fDoMesonQA > 0){
2285                                         if (isTruePi0){
2286                                                 if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
2287                                                         fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2288                                                         fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2289                                                         fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2290                                                 }
2291                                         } else if (isTrueEta){
2292                                                 if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
2293                                                         fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2294                                                         fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2295                                                         fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2296                                                 }
2297                                         }
2298                                 }
2299                                 if(gamma0MotherLabel >= MCStack->GetNprimary()){ // Secondary Meson
2300                                         Int_t secMotherLabel = ((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0);
2301                                         Float_t weightedSec= 1;
2302                                         if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, fMCStack, fInputEvent) && MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
2303                                                 weightedSec= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, fMCStack, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt
2304                                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2305                                         }
2306                                         fHistoTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2307                                         fCharMesonMCInfo = 2;
2308                                         if (secMotherLabel >-1){
2309                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
2310                                                         fCharMesonMCInfo = 4;
2311                                                         fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2312                                                         if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2313                                                 }
2314                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==221){
2315                                                 fCharMesonMCInfo = 3;
2316                                                 fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2317                                                 if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2318                                                 }
2319                                                 if(MCStack->Particle(secMotherLabel)->GetPdgCode()==3122){
2320                                                         fCharMesonMCInfo = 7;
2321                                                         fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2322                                                         if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(MCStack->Particle(secMotherLabel)->Pt());
2323                                                 }
2324                                         }
2325                                 } else { // Only primary pi0 for efficiency calculation
2326                                         fCharMesonMCInfo = 6;
2327                                         Float_t weighted= 1;
2328                                         if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCStack, fInputEvent)){
2329                                                 if (((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt()>0.005){
2330                                                         weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, fMCStack, fInputEvent);
2331                                                         //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
2332                                                 }
2333                                         }
2334                                         fHistoTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2335                                         fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2336                                         fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2337                                         
2338                                         
2339                                         if (fDoMesonQA > 0){
2340                                                 if(isTruePi0){ // Only primary pi0 for resolution
2341                                                         fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
2342                                                 }
2343                                                 if (isTrueEta){ // Only primary eta for resolution
2344                                                         fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
2345                                                 }
2346                                         }
2347                                 }
2348                         } else if(!isTruePi0 && !isTrueEta){ // Background
2349                                 if (fDoMesonQA > 0){
2350                                         if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2351                                                 fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2352                                                 fCharMesonMCInfo = 1;
2353                                         } else { // No photon or without mother
2354                                                 fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2355                                         }
2356                                 }
2357                                 if( isTruePi0Dalitz || isTrueEtaDalitz ){
2358                                         // Dalitz
2359                                         fCharMesonMCInfo = 5;
2360                                         fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2361                                 } else if (gamma0DalitzCand || gamma1DalitzCand){
2362                                         if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2363                                 }
2364                         }
2365                 }
2366         }
2367 }
2368 //______________________________________________________________________
2369 void AliAnalysisTaskGammaConvCalo::ProcessTrueMesonCandidatesAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
2370 {
2371         
2372         // Process True Mesons
2373         TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
2374         Bool_t isTruePi0 = kFALSE;
2375         Bool_t isTrueEta = kFALSE;
2376         Bool_t isTruePi0Dalitz = kFALSE;
2377         Bool_t isTrueEtaDalitz = kFALSE;
2378         Bool_t gamma0DalitzCand = kFALSE;
2379         Bool_t gamma1DalitzCand = kFALSE;
2380         
2381         AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
2382         AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
2383         
2384         fCharMesonMCInfo = 0;
2385         Int_t gamma0MCLabel = -1;
2386         Int_t gamma0MotherLabel = -1;
2387         if(!positiveMC||!negativeMC)
2388                 return;
2389         
2390         if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
2391                 gamma0MCLabel = positiveMC->GetMother();
2392         }
2393         
2394         if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2395                 // Daughters Gamma 0
2396                 AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
2397                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2398                         if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
2399                                 if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
2400                                 gamma0MotherLabel=gammaMC0->GetMother();
2401                                 }
2402                         }
2403                         if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
2404                                 gamma0DalitzCand = kTRUE;
2405                                 gamma0MotherLabel=-111;
2406                         }
2407                         if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
2408                                 gamma0DalitzCand = kTRUE;
2409                                 gamma0MotherLabel=-221;
2410                         }
2411                 }
2412         }
2413         positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
2414         negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
2415         
2416         Int_t gamma1MCLabel = -1;
2417         Int_t gamma1MotherLabel = -1;
2418         if(!positiveMC||!negativeMC)
2419                 return;
2420         
2421         if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
2422                 gamma1MCLabel = positiveMC->GetMother();
2423         }
2424         if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
2425                 // Daughters Gamma 1
2426                 AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
2427                 if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
2428                         if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...
2429                                 if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
2430                                 gamma1MotherLabel=gammaMC1->GetMother();
2431                                 }
2432                         }
2433                         if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
2434                                 gamma1DalitzCand = kTRUE;
2435                                 gamma1MotherLabel=-111;
2436                         }
2437                         if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
2438                                 gamma1DalitzCand = kTRUE;
2439                                 gamma1MotherLabel=-221;
2440                         }
2441                 }
2442         }
2443         if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
2444                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
2445                         isTruePi0=kTRUE;
2446                 }
2447                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 221){
2448                         isTrueEta=kTRUE;
2449                 }
2450         }
2451         
2452         //Identify Dalitz candidate
2453         if (gamma1DalitzCand || gamma0DalitzCand){
2454                 if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
2455                         if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2456                         if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2457                 }
2458                 if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
2459                         if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2460                         if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;
2461                 }
2462         }
2463         
2464         if(isTruePi0 || isTrueEta){// True Pion or Eta
2465                 fHistoTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2466                 if (fDoMesonQA > 0){
2467                         if (isTruePi0){
2468                                 if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
2469                                 fHistoTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2470                                 fHistoTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2471                                 fHistoTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2472                                 }
2473                         } else if (isTrueEta){
2474                                 if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
2475                                 fHistoTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2476                                 fHistoTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha());
2477                                 fHistoTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle());
2478                                 }
2479                         }
2480                 }
2481                 if(!(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel))->IsPrimary())){ // Secondary Meson
2482                         Int_t secMotherLabel = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetMother();
2483                         Float_t weightedSec= 1;
2484                         if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, 0x0, fInputEvent) && static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2485                                 weightedSec= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),secMotherLabel, 0x0, fInputEvent)/2.; //invariant mass is additive thus the weight for the daughters has to be devide by two for the K0s at a certain pt
2486                                 //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2487                         }
2488                         fHistoTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2489                         fCharMesonMCInfo = 2;
2490                         if (secMotherLabel >-1){
2491                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2492                                         fCharMesonMCInfo = 4;
2493                                         fHistoTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2494                                         if (fDoMesonQA > 0)fHistoTrueK0sWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2495                                 }
2496                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==221){
2497                                         fCharMesonMCInfo = 3;
2498                                         fHistoTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2499                                         if (fDoMesonQA > 0)fHistoTrueEtaWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2500                                 }
2501                                 if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==3122){
2502                                         fCharMesonMCInfo = 7;
2503                                         fHistoTrueSecondaryMotherFromLambdaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2504                                         if (fDoMesonQA > 0)fHistoTrueLambdaWithPi0DaughterMCPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2505                                 }
2506                         }
2507                 }else{ // Only primary pi0 for efficiency calculation
2508                         Float_t weighted= 1;
2509                         fCharMesonMCInfo = 6;
2510                         if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
2511                                 if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
2512                                 weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, 0x0, fInputEvent);
2513                                 //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
2514                                 }
2515                         }
2516                         fHistoTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2517                         fHistoTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2518                         fProfileTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2519                         
2520                         if (fDoMesonQA > 0){
2521                                 if(isTruePi0){ // Only primary pi0 for resolution
2522                                         fHistoTruePrimaryPi0MCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2523                                                                                                                         (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2524                                 
2525                                 }
2526                                 if (isTrueEta){ // Only primary eta for resolution
2527                                         fHistoTruePrimaryEtaMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2528                                                                                                                         (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2529                                 }
2530                         }
2531                 }
2532         } else if(!isTruePi0 && !isTrueEta) { // Background
2533                 if (fDoMesonQA > 0){
2534                         if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2535                                 fHistoTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2536                                 fCharMesonMCInfo = 1;
2537                         } else { // No photon or without mother
2538                                 fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2539                         }
2540                 }
2541                 if( isTruePi0Dalitz || isTrueEtaDalitz ){
2542                         // Dalitz
2543                         fCharMesonMCInfo = 5;
2544                         fHistoTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2545                 } else if (gamma0DalitzCand || gamma1DalitzCand){
2546                         if (fDoMesonQA > 0)fHistoTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2547                 }
2548         }
2549 }
2550 //________________________________________________________________________
2551 void AliAnalysisTaskGammaConvCalo::CalculateBackground(){
2552         
2553         Int_t zbin= fBGClusHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2554         Int_t mbin = 0;
2555         
2556         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2557                 mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2558         } else {
2559                 mbin = fBGClusHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2560         }
2561         
2562         AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;    
2563         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2564                 for(Int_t nEventsInBG=0;nEventsInBG<fBGClusHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2565                         AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2566                         if(fMoveParticleAccordingToVertex == kTRUE){
2567                                 bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2568                         }
2569                         
2570                         for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2571                                 AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2572                                 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2573                                         AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2574                                         if(fMoveParticleAccordingToVertex == kTRUE){
2575                                                 MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2576                                         }
2577                                         if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2578                                                 RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2579                                         }
2580                                         
2581                                         AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2582                                         backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2583                                         if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2584                                                 ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2585                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2586                                                 Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2587                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2588                                         }
2589                                         delete backgroundCandidate;
2590                                         backgroundCandidate = 0x0;
2591                                 }
2592                         }
2593                 }
2594         } else {
2595                 for(Int_t nEventsInBG=0;nEventsInBG <fBGClusHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2596                         AliGammaConversionAODVector *previousEventV0s = fBGClusHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2597                         if(previousEventV0s){
2598                                 if(fMoveParticleAccordingToVertex == kTRUE){
2599                                         bgEventVertex = fBGClusHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2600                                 }
2601                                 for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2602                                         AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2603                                         for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2604                                 
2605                                                 AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2606                                         
2607                                                 if(fMoveParticleAccordingToVertex == kTRUE){
2608                                                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2609                                                 }
2610                                                 if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2611                                                         RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2612                                                 }
2613                                         
2614                                         
2615                                                 AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2616                                                 backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2617                                                 if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2618                                                         fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2619                                                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2620                                                         fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2621                                                 }
2622                                                 delete backgroundCandidate;
2623                                                 backgroundCandidate = 0x0;
2624                                         }
2625                                 }
2626                         }
2627                 }
2628         }
2629 }
2630
2631 //________________________________________________________________________
2632 void AliAnalysisTaskGammaConvCalo::CalculateBackgroundRP(){
2633         
2634         Int_t zbin= fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2635         Int_t mbin = 0;
2636         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2637                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2638         } else {
2639                 mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2640         }
2641         
2642         
2643         //Rotation Method
2644         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
2645                 // Correct for the number of rotations
2646                 // BG is for rotation the same, except for factor NRotations
2647                 Double_t weight=1./Double_t(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents());
2648                 
2649                 for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2650                         
2651                         AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2652                         if (gamma0==NULL) continue;
2653                         for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGammaCandidates->GetEntries();secondGammaIndex++){
2654                                 AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(secondGammaIndex));
2655                                 if (gamma1 == NULL) continue;
2656                                 if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(gamma1,fInputEvent))continue;
2657                                 for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){
2658                                 
2659                                         RotateParticle(gamma1);                                 
2660                                         AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2661                                         backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2662                                         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2663                                                 ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
2664                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2665                                                 Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2666                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2667                                         }
2668                                 }
2669                         }
2670                 }
2671         }
2672         else{
2673                 // Do Event Mixing
2674                 for(Int_t nEventsInBG=0;nEventsInBG <fBGHandlerRP[fiCut]->GetNBGEvents(fGammaCandidates,fInputEvent);nEventsInBG++){
2675                 
2676                         AliGammaConversionPhotonVector *previousEventGammas = fBGHandlerRP[fiCut]->GetBGGoodGammas(fGammaCandidates,fInputEvent,nEventsInBG);
2677                         
2678                         if(previousEventGammas){
2679                                 // test weighted background
2680                                 Double_t weight=1.0;
2681                                 // Correct for the number of eventmixing:
2682                                 // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1))  using sum formula sum(i)=N*(N-1)/2  -> N*(N-1)/2
2683                                 // real combinations (since you cannot combine a photon with its own)
2684                                 // but BG leads to N_{a}*N_{b} combinations
2685                                 weight*=0.5*(Double_t(fGammaCandidates->GetEntries()-1))/Double_t(previousEventGammas->size());
2686                                 
2687                                 for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2688                                         AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));                                     
2689                                         for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
2690                                                 
2691                                                 AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
2692                                                 
2693                                                 AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2694                                                 backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2695                                                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2696                                                 ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
2697                                                 fHistoMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2698                                                 Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2699                                                 fSparseMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2700                                                 }
2701                                         }
2702                                 }
2703                         }
2704                 }
2705         }
2706 }
2707 //________________________________________________________________________
2708 void AliAnalysisTaskGammaConvCalo::RotateParticle(AliAODConversionPhoton *gamma){
2709         Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation();
2710         Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
2711         Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2712         gamma->RotateZ(rotationValue);
2713 }
2714
2715 //________________________________________________________________________
2716 void AliAnalysisTaskGammaConvCalo::RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP){
2717         
2718         previousEventEP=previousEventEP+TMath::Pi();
2719         thisEventEP=thisEventEP+TMath::Pi();
2720         Double_t rotationValue= thisEventEP-previousEventEP;
2721         gamma->RotateZ(rotationValue);
2722 }
2723
2724 //________________________________________________________________________
2725 void AliAnalysisTaskGammaConvCalo::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
2726         //see header file for documentation
2727         
2728         Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
2729         Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
2730         Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
2731         
2732         Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
2733         particle->SetConversionPoint(movedPlace);
2734 }
2735 //________________________________________________________________________
2736 void AliAnalysisTaskGammaConvCalo::UpdateEventByEventData(){
2737         //see header file for documentation
2738         if(fGammaCandidates->GetEntries() >0 ){
2739                 if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2740                         fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
2741                         fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
2742                 } else { // means we use #V0s for multiplicity
2743                         fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
2744                         fBGClusHandler[fiCut]->AddEvent(fClusterCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
2745                 }
2746         }
2747 }
2748
2749
2750 //________________________________________________________________________
2751 void AliAnalysisTaskGammaConvCalo::FillPhotonCombinatorialBackgroundHist(AliAODConversionPhoton *TruePhotonCandidate, Int_t pdgCode[])
2752 {
2753         // Combinatorial Bck = 0 ee, 1 ep,i 2 ek, 3 ep, 4 emu, 5 pipi, 6 pik, 7 pip, 8 pimu, 9 kk, 10 kp, 11 kmu, 12 pp, 13 pmu, 14 mumu, 15 Rest
2754         if(pdgCode[0]==11   && pdgCode[1]==11){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),0);}
2755         else if( (pdgCode[0]==11   && pdgCode[1]==211) || (pdgCode[0]==211  && pdgCode[1]==11) )
2756         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),1);}
2757         else if( (pdgCode[0]==11   && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==11) )
2758         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),2);}
2759         else if( (pdgCode[0]==11   && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==11) )
2760         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),3);}
2761         else if( (pdgCode[0]==11   && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==11) )
2762         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),4);}
2763         else if(  pdgCode[0]==211  && pdgCode[1]==211 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),5);}
2764         else if( (pdgCode[0]==211  && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==211) )
2765         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),6);}
2766         else if( (pdgCode[0]==211  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==211) )
2767         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),7);}
2768         else if( (pdgCode[0]==211  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==211) )
2769         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),8);}
2770         else if(  pdgCode[0]==321  && pdgCode[1]==321 ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),9);}
2771         else if( (pdgCode[0]==321  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==321) )
2772         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),10);}
2773         else if( (pdgCode[0]==321  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==321) )
2774         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),11);}
2775         else if(  pdgCode[0]==2212   && pdgCode[1]==2212  ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),12);}
2776         else if( (pdgCode[0]==2212  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==2212) )
2777         {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),13);}
2778         else if(  pdgCode[0]==13   && pdgCode[1]==13  ){if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),14);}
2779         else {if(fIsFromMBHeader)fHistoCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),15);}
2780 }
2781 //________________________________________________________________________
2782 void AliAnalysisTaskGammaConvCalo::RelabelAODPhotonCandidates(Bool_t mode){
2783         
2784         // Relabeling For AOD Event
2785         // ESDiD -> AODiD
2786         // MCLabel -> AODMCLabel
2787         
2788         if(mode){
2789                 fMCStackPos = new Int_t[fReaderGammas->GetEntries()];
2790                 fMCStackNeg = new Int_t[fReaderGammas->GetEntries()];
2791                 fESDArrayPos = new Int_t[fReaderGammas->GetEntries()];
2792                 fESDArrayNeg = new Int_t[fReaderGammas->GetEntries()];
2793         }
2794         
2795         for(Int_t iGamma = 0;iGamma<fReaderGammas->GetEntries();iGamma++){
2796                 AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(iGamma);
2797                 if(!PhotonCandidate) continue;
2798                 if(!mode){// Back to ESD Labels
2799                 PhotonCandidate->SetMCLabelPositive(fMCStackPos[iGamma]);
2800                 PhotonCandidate->SetMCLabelNegative(fMCStackNeg[iGamma]);
2801                 PhotonCandidate->SetLabelPositive(fESDArrayPos[iGamma]);
2802                 PhotonCandidate->SetLabelNegative(fESDArrayNeg[iGamma]);
2803                 continue;
2804                 }
2805                 fMCStackPos[iGamma] =  PhotonCandidate->GetMCLabelPositive();
2806                 fMCStackNeg[iGamma] =  PhotonCandidate->GetMCLabelNegative();
2807                 fESDArrayPos[iGamma] = PhotonCandidate->GetTrackLabelPositive();
2808                 fESDArrayNeg[iGamma] = PhotonCandidate->GetTrackLabelNegative();
2809                 
2810                 Bool_t AODLabelPos = kFALSE;
2811                 Bool_t AODLabelNeg = kFALSE;
2812                 
2813                 for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
2814                 AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
2815                 if(!AODLabelPos){
2816                         if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
2817                         PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
2818                         PhotonCandidate->SetLabelPositive(i);
2819                         AODLabelPos = kTRUE;
2820                         }
2821                 }
2822                 if(!AODLabelNeg){
2823                         if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
2824                         PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
2825                         PhotonCandidate->SetLabelNegative(i);
2826                         AODLabelNeg = kTRUE;
2827                         }
2828                 }
2829                 if(AODLabelNeg && AODLabelPos){
2830                         break;
2831                 }
2832                 }
2833                 if(!AODLabelPos || !AODLabelNeg){
2834                 cout<<"WARNING!!! AOD TRACKS NOT FOUND FOR"<<endl;
2835                 }
2836         }
2837         
2838         
2839         if(!mode){
2840                 delete[] fMCStackPos;
2841                 delete[] fMCStackNeg;
2842                 delete[] fESDArrayPos;
2843                 delete[] fESDArrayNeg;
2844         }
2845 }
2846
2847 void AliAnalysisTaskGammaConvCalo::SetLogBinningXTH2(TH2* histoRebin){
2848         TAxis *axisafter = histoRebin->GetXaxis();
2849         Int_t bins = axisafter->GetNbins();
2850         Double_t from = axisafter->GetXmin();
2851         Double_t to = axisafter->GetXmax();
2852         Double_t *newbins = new Double_t[bins+1];
2853         newbins[0] = from;
2854         Double_t factor = TMath::Power(to/from, 1./bins);
2855         for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1];
2856         axisafter->Set(bins, newbins);
2857         delete [] newbins;
2858 }
2859
2860 //________________________________________________________________________
2861 void AliAnalysisTaskGammaConvCalo::Terminate(const Option_t *)
2862 {
2863   
2864   //fOutputContainer->Print(); // Will crash on GRID
2865 }
2866
2867 //________________________________________________________________________
2868 Int_t AliAnalysisTaskGammaConvCalo::GetSourceClassification(Int_t daughter, Int_t pdgCode){
2869   
2870         if (daughter == 111) {
2871                 if (abs(pdgCode) == 310) return 1; // k0s
2872                 else if (abs(pdgCode) == 3122) return 2; // Lambda
2873                 else if (abs(pdgCode) == 130) return 3; // K0L
2874                 else if (abs(pdgCode) == 2212) return 4; // proton
2875                 else if (abs(pdgCode) == 2112) return 5; // neutron
2876                 else if (abs(pdgCode) == 211) return 6; // pion
2877                 else if (abs(pdgCode) == 321) return 7; // kaon
2878                 else if (abs(pdgCode) == 113 || abs(pdgCode) == 213 ) return 8; // rho 0,+,-
2879                 else if (abs(pdgCode) == 3222 || abs(pdgCode) == 3212 || abs(pdgCode) == 3112  ) return 9; // Sigma
2880                 else if (abs(pdgCode) == 2224 || abs(pdgCode) == 2214 || abs(pdgCode) == 2114 || abs(pdgCode) == 1114  ) return 10; // Delta
2881                 else if (abs(pdgCode) == 313 || abs(pdgCode) == 323   ) return 11; // K*
2882                 else return 15;
2883         }
2884         return 15;
2885   
2886 }
2887
2888 // //________________________________________________________________________
2889 // Double_t AliAnalysisTaskGammaConvCalo::GetMaxCellEnergy(const AliVCluster *cluster, Short_t &id) const
2890 // {
2891 //   // Get maximum energy of attached cell.
2892 //   
2893 //   id = -1;
2894 //   Double_t maxe = 0;
2895 //   Int_t ncells = cluster->GetNCells();
2896 //   if (fEsdCells) {
2897 //     for (Int_t i=0; i<ncells; i++) {
2898 //       Double_t e = fEsdCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2899 //       if (e>maxe) {
2900 //         maxe = e;
2901 //         id   = cluster->GetCellAbsId(i);
2902 //       }
2903 //     }
2904 //   } else {
2905 //     for (Int_t i=0; i<ncells; i++) {
2906 //       Double_t e = fAodCells->GetCellAmplitude(TMath::Abs(cluster->GetCellAbsId(i)));
2907 //       if (e>maxe)
2908 //         maxe = e;
2909 //       id   = cluster->GetCellAbsId(i);
2910 //     }
2911 //   }
2912 //   return maxe;
2913 // }