0955548f5c514f996c3daac4298c017f0c47acce
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskGammaConvV1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: Martin Wilde, Daniel Lohner, Friederike Bock                   *
5  * Version 1.0                                                            *
6  *                                                                        *
7  * based on: on older version (see aliroot up to v5-04-42-AN)             *
8  *           AliAnalysisTaskGammaConversion.cxx                           *
9  *           Authors: Kathrin Koch, Kenneth Aamodt, Ana Marin             *
10  *                                                                        *
11  * Permission to use, copy, modify and distribute this software and its   *
12  * documentation strictly for non-commercial purposes is hereby granted   *
13  * without fee, provided that the above copyright notice appears in all   *
14  * copies and that both the copyright notice and this permission notice   *
15  * appear in the supporting documentation. The authors make no claims     *
16  * about the suitability of this software for any purpose. It is          *
17  * provided "as is" without express or implied warranty.                  *
18  **************************************************************************/
19
20 ////////////////////////////////////////////////
21 //---------------------------------------------
22
23 // Class used to do analysis on conversion pairs
24 //---------------------------------------------
25 ///////////////////////////////////////////////
26 #include "TChain.h"
27 #include "TTree.h"
28 #include "TBranch.h"
29 #include "TFile.h"
30 #include "TH1F.h"
31 #include "TH2F.h"
32 #include "TH3F.h"
33 #include "THnSparse.h"
34 #include "TCanvas.h"
35 #include "TNtuple.h"
36 #include "AliAnalysisTask.h"
37 #include "AliAnalysisManager.h"
38 #include "AliESDEvent.h"
39 #include "AliESDInputHandler.h"
40 #include "AliMCEventHandler.h"
41 #include "AliMCEvent.h"
42 #include "AliMCParticle.h"
43 #include "AliCentrality.h"
44 #include "AliESDVZERO.h"
45 #include "AliESDpid.h"
46 #include "AliAnalysisTaskGammaConvV1.h"
47 #include "AliVParticle.h"
48 #include "AliESDtrack.h"
49 #include "AliESDtrackCuts.h"
50 #include "AliKFVertex.h"
51 #include "AliV0ReaderV1.h"
52 #include "AliGenCocktailEventHeader.h"
53 #include "AliConversionAODBGHandlerRP.h"
54 #include "AliAODMCParticle.h"
55 #include "AliAODMCHeader.h"
56 #include "AliEventplane.h"
57
58 ClassImp(AliAnalysisTaskGammaConvV1)
59
60 //________________________________________________________________________
61 AliAnalysisTaskGammaConvV1::AliAnalysisTaskGammaConvV1(): AliAnalysisTaskSE(),
62    fV0Reader(NULL),
63    fBGHandler(NULL),
64    fBGHandlerRP(NULL),
65    fInputEvent(NULL),
66    fMCEvent(NULL),
67    fMCStack(NULL),
68    fCutFolder(NULL),
69    fESDList(NULL),
70    fBackList(NULL),
71    fMotherList(NULL),
72    fPhotonDCAList(NULL),
73    fMesonDCAList(NULL),        
74    fTrueList(NULL),
75    fMCList(NULL),
76    fHeaderNameList(NULL),
77    fOutputContainer(0),
78    fReaderGammas(NULL),
79    fGammaCandidates(NULL),
80    fCutArray(NULL),
81    fConversionCuts(NULL),
82    fMesonCutArray(NULL),
83    fMesonCuts(NULL),
84    hESDConvGammaPt(NULL),
85    hESDConvGammaR(NULL),
86    hESDConvGammaEta(NULL),
87    tESDConvGammaPtDcazCat(NULL),
88    fPtGamma(0),
89    fDCAzPhoton(0),
90    fRConvPhoton(0),
91    fEtaPhoton(0),
92    iCatPhoton(0),
93    iPhotonMCInfo(0),
94    hESDMotherInvMassPt(NULL),
95    sESDMotherInvMassPtZM(NULL),
96    hESDMotherBackInvMassPt(NULL),
97    sESDMotherBackInvMassPtZM(NULL),
98    hESDMotherInvMassEalpha(NULL),
99    hESDMotherPi0PtY(NULL),
100    hESDMotherEtaPtY(NULL),
101    hESDMotherPi0PtAlpha(NULL),
102    hESDMotherEtaPtAlpha(NULL),
103    hESDMotherPi0PtOpenAngle(NULL),
104    hESDMotherEtaPtOpenAngle(NULL),
105    hMCHeaders(NULL),
106    hMCAllGammaPt(NULL),
107    hMCDecayGammaPi0Pt(NULL),
108    hMCDecayGammaRhoPt(NULL),
109    hMCDecayGammaEtaPt(NULL),
110    hMCDecayGammaOmegaPt(NULL),
111    hMCDecayGammaEtapPt(NULL),
112    hMCDecayGammaPhiPt(NULL),
113    hMCDecayGammaSigmaPt(NULL),
114    hMCConvGammaPt(NULL),
115    hMCConvGammaR(NULL),
116    hMCConvGammaEta(NULL),
117    hMCPi0Pt(NULL),
118    hMCPi0WOWeightPt(NULL),
119    hMCEtaPt(NULL),
120    hMCEtaWOWeightPt(NULL),
121    hMCPi0InAccPt(NULL),
122    hMCEtaInAccPt(NULL),
123    hMCPi0PtY(NULL),
124    hMCEtaPtY(NULL),
125    hMCK0sPt(NULL),
126    hMCK0sWOWeightPt(NULL),
127    hMCK0sPtY(NULL),
128    hESDTrueMotherInvMassPt(NULL),
129    hESDTruePrimaryMotherInvMassPt(NULL),
130    hESDTruePrimaryMotherW0WeightingInvMassPt(NULL),
131    pESDTruePrimaryMotherWeightsInvMassPt(NULL),
132    hESDTruePrimaryPi0MCPtResolPt(NULL),
133    hESDTruePrimaryEtaMCPtResolPt(NULL),
134    hESDTrueSecondaryMotherInvMassPt(NULL),
135    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
136    hESDTrueK0sWithPi0DaughterMCPt(NULL),
137    hESDTrueSecondaryMotherFromEtaInvMassPt(NULL),
138    hESDTrueEtaWithPi0DaughterMCPt(NULL),
139    hESDTrueBckGGInvMassPt(NULL),
140    hESDTrueBckContInvMassPt(NULL),
141    hESDTruePi0PtY(NULL),
142    hESDTrueEtaPtY(NULL),
143    hESDTruePi0PtAlpha(NULL),
144    hESDTrueEtaPtAlpha(NULL),
145    hESDTruePi0PtOpenAngle(NULL),
146    hESDTrueEtaPtOpenAngle(NULL),
147    hESDTrueMotherDalitzInvMassPt(NULL),
148    hESDTrueConvGammaPt(NULL),
149    hESDCombinatorialPt(NULL),
150    hESDTruePrimaryConvGammaPt(NULL),
151    hESDTruePrimaryConvGammaESDPtMCPt(NULL),
152    hESDTrueSecondaryConvGammaPt(NULL),
153    hESDTrueSecondaryConvGammaFromXFromK0sPt(NULL),
154    hESDTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
155    hESDTrueDalitzPsiPairDeltaPhi(NULL),
156    hESDTrueGammaPsiPairDeltaPhi(NULL),
157    hNEvents(NULL),
158    hNGoodESDTracks(NULL),
159    hNGammaCandidates(NULL),
160    hNV0Tracks(NULL),
161    hEtaShift(NULL),
162    tESDMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
163    fInvMass(0),
164    fPt(0),
165    fDCAzGammaMin(0),
166    fDCAzGammaMax(0),
167    iFlag(0),
168    iMesonMCInfo(0),
169    fEventPlaneAngle(-100),
170    fRandom(0),
171    fnGammaCandidates(0),
172    fUnsmearedPx(NULL),
173    fUnsmearedPy(NULL),
174    fUnsmearedPz(NULL),
175    fUnsmearedE(NULL),
176    fMCStackPos(NULL),
177    fMCStackNeg(NULL),
178    fESDArrayPos(NULL),
179    fESDArrayNeg(NULL),
180    fnCuts(0),
181    fiCut(0),
182    fMoveParticleAccordingToVertex(kTRUE),
183    fIsHeavyIon(kFALSE),
184    fDoMesonAnalysis(kTRUE),
185    fDoMesonQA(0),
186    fDoPhotonQA(0),
187    fIsFromMBHeader(kTRUE),
188    fIsMC(kFALSE)
189 {
190
191 }
192
193 //________________________________________________________________________
194 AliAnalysisTaskGammaConvV1::AliAnalysisTaskGammaConvV1(const char *name):
195    AliAnalysisTaskSE(name),
196    fV0Reader(NULL),
197    fBGHandler(NULL),
198    fBGHandlerRP(NULL),
199    fInputEvent(NULL),
200    fMCEvent(NULL),
201    fMCStack(NULL),
202    fCutFolder(NULL),
203    fESDList(NULL),
204    fBackList(NULL),
205    fMotherList(NULL),
206    fPhotonDCAList(NULL),
207    fMesonDCAList(NULL),
208    fTrueList(NULL),
209    fMCList(NULL),
210    fHeaderNameList(NULL),
211    fOutputContainer(0),
212    fReaderGammas(NULL),
213    fGammaCandidates(NULL),
214    fCutArray(NULL),
215    fConversionCuts(NULL),
216    fMesonCutArray(NULL),
217    fMesonCuts(NULL),
218    hESDConvGammaPt(NULL),
219    hESDConvGammaR(NULL),
220    hESDConvGammaEta(NULL),
221    tESDConvGammaPtDcazCat(NULL),
222    fPtGamma(0),
223    fDCAzPhoton(0),
224    fRConvPhoton(0),
225    fEtaPhoton(0),
226    iCatPhoton(0),
227    iPhotonMCInfo(0),
228    hESDMotherInvMassPt(NULL),
229    sESDMotherInvMassPtZM(NULL),
230    hESDMotherBackInvMassPt(NULL),
231    sESDMotherBackInvMassPtZM(NULL),
232    hESDMotherInvMassEalpha(NULL),
233    hESDMotherPi0PtY(NULL),
234    hESDMotherEtaPtY(NULL),
235    hESDMotherPi0PtAlpha(NULL),
236    hESDMotherEtaPtAlpha(NULL),
237    hESDMotherPi0PtOpenAngle(NULL),
238    hESDMotherEtaPtOpenAngle(NULL),
239    hMCHeaders(NULL),
240    hMCAllGammaPt(NULL),
241    hMCDecayGammaPi0Pt(NULL),
242    hMCDecayGammaRhoPt(NULL),
243    hMCDecayGammaEtaPt(NULL),
244    hMCDecayGammaOmegaPt(NULL),
245    hMCDecayGammaEtapPt(NULL),
246    hMCDecayGammaPhiPt(NULL),
247    hMCDecayGammaSigmaPt(NULL),
248    hMCConvGammaPt(NULL),
249    hMCConvGammaR(NULL),
250    hMCConvGammaEta(NULL),
251    hMCPi0Pt(NULL),
252    hMCPi0WOWeightPt(NULL),
253    hMCEtaPt(NULL),
254    hMCEtaWOWeightPt(NULL),
255    hMCPi0InAccPt(NULL),
256    hMCEtaInAccPt(NULL),
257    hMCPi0PtY(NULL),
258    hMCEtaPtY(NULL),
259    hMCK0sPt(NULL),
260    hMCK0sWOWeightPt(NULL),
261    hMCK0sPtY(NULL),
262    hESDTrueMotherInvMassPt(NULL),
263    hESDTruePrimaryMotherInvMassPt(NULL),
264    hESDTruePrimaryMotherW0WeightingInvMassPt(NULL),
265    pESDTruePrimaryMotherWeightsInvMassPt(NULL),
266    hESDTruePrimaryPi0MCPtResolPt(NULL),
267    hESDTruePrimaryEtaMCPtResolPt(NULL),
268    hESDTrueSecondaryMotherInvMassPt(NULL),
269    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
270    hESDTrueK0sWithPi0DaughterMCPt(NULL),
271    hESDTrueSecondaryMotherFromEtaInvMassPt(NULL),
272    hESDTrueEtaWithPi0DaughterMCPt(NULL),
273    hESDTrueBckGGInvMassPt(NULL),
274    hESDTrueBckContInvMassPt(NULL),
275    hESDTruePi0PtY(NULL),
276    hESDTrueEtaPtY(NULL),
277    hESDTruePi0PtAlpha(NULL),
278    hESDTrueEtaPtAlpha(NULL),
279    hESDTruePi0PtOpenAngle(NULL),
280    hESDTrueEtaPtOpenAngle(NULL),
281    hESDTrueMotherDalitzInvMassPt(NULL),
282    hESDTrueConvGammaPt(NULL),
283    hESDCombinatorialPt(NULL),
284    hESDTruePrimaryConvGammaPt(NULL),
285    hESDTruePrimaryConvGammaESDPtMCPt(NULL),
286    hESDTrueSecondaryConvGammaPt(NULL),
287    hESDTrueSecondaryConvGammaFromXFromK0sPt(NULL),
288    hESDTrueSecondaryConvGammaFromXFromLambdaPt(NULL),
289    hESDTrueDalitzPsiPairDeltaPhi(NULL),
290    hESDTrueGammaPsiPairDeltaPhi(NULL),
291    hNEvents(NULL),
292    hNGoodESDTracks(NULL),
293    hNGammaCandidates(NULL),
294    hNV0Tracks(NULL),
295    hEtaShift(NULL),
296    tESDMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
297    fInvMass(0),
298    fPt(0),
299    fDCAzGammaMin(0),
300    fDCAzGammaMax(0),
301    iFlag(0),
302    iMesonMCInfo(0),
303    fEventPlaneAngle(-100),
304    fRandom(0),
305    fnGammaCandidates(0),
306    fUnsmearedPx(NULL),
307    fUnsmearedPy(NULL),
308    fUnsmearedPz(NULL),
309    fUnsmearedE(NULL),
310    fMCStackPos(NULL),
311    fMCStackNeg(NULL),
312    fESDArrayPos(NULL),
313    fESDArrayNeg(NULL),
314    fnCuts(0),
315    fiCut(0),
316    fMoveParticleAccordingToVertex(kTRUE),
317    fIsHeavyIon(kFALSE),
318    fDoMesonAnalysis(kTRUE),
319    fDoMesonQA(0),
320    fDoPhotonQA(0),
321    fIsFromMBHeader(kTRUE),
322    fIsMC(kFALSE)
323 {
324    // Define output slots here
325    DefineOutput(1, TList::Class());
326 }
327
328 AliAnalysisTaskGammaConvV1::~AliAnalysisTaskGammaConvV1()
329 {
330    if(fGammaCandidates){
331       delete fGammaCandidates;
332       fGammaCandidates = 0x0;
333    }
334    if(fBGHandler){
335       delete[] fBGHandler;
336       fBGHandler = 0x0;
337    }
338    if(fBGHandlerRP){
339       delete[] fBGHandlerRP;
340       fBGHandlerRP = 0x0;
341    }
342 }
343 //___________________________________________________________
344 void AliAnalysisTaskGammaConvV1::InitBack(){
345
346    const Int_t nDim = 4;
347    Int_t nBins[nDim] = {800,250,7,4};
348    Double_t xMin[nDim] = {0,0, 0,0};
349    Double_t xMax[nDim] = {0.8,25,7,4};
350    
351    sESDMotherInvMassPtZM = new THnSparseF*[fnCuts];
352    sESDMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
353
354    fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
355    fBGHandlerRP = new AliConversionAODBGHandlerRP*[fnCuts];
356    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
357       if (((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
358          TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
359          TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
360          
361          Int_t collisionSystem = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(0,1));
362          Int_t centMin = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(1,1));
363          Int_t centMax = atoi((TString)(((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber())(2,1));
364          
365          if(collisionSystem == 1 || collisionSystem == 2 ||
366             collisionSystem == 5 || collisionSystem == 8 ||
367             collisionSystem == 9){
368             centMin = centMin*10;
369             centMax = centMax*10; 
370             if(centMax ==0 && centMax!=centMin) centMax=100;
371          }
372          else if(collisionSystem == 3 || collisionSystem == 6){
373             centMin = centMin*5;
374             centMax = centMax*5;
375          }
376          else if(collisionSystem == 4 || collisionSystem == 7){
377             centMin = ((centMin*5)+45);
378             centMax = ((centMax*5)+45);
379          }
380          
381          fBackList[iCut] = new TList();
382          fBackList[iCut]->SetName(Form("%s_%s Back histograms",cutstring.Data(),cutstringMeson.Data()));
383          fBackList[iCut]->SetOwner(kTRUE);
384          fCutFolder[iCut]->Add(fBackList[iCut]);
385
386          sESDMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
387          fBackList[iCut]->Add(sESDMotherBackInvMassPtZM[iCut]);
388
389          fMotherList[iCut] = new TList();
390          fMotherList[iCut]->SetName(Form("%s_%s Mother histograms",cutstring.Data(),cutstringMeson.Data()));
391          fMotherList[iCut]->SetOwner(kTRUE);
392          fCutFolder[iCut]->Add(fMotherList[iCut]);
393
394          sESDMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
395          fMotherList[iCut]->Add(sESDMotherInvMassPtZM[iCut]);
396
397          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
398             fBGHandler[iCut] = new AliGammaConversionAODBGHandler(
399                                                                   collisionSystem,centMin,centMax,
400                                                                   ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents(),
401                                                                   ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity());
402             fBGHandlerRP[iCut] = NULL;
403          }
404          else{
405             fBGHandlerRP[iCut] = new AliConversionAODBGHandlerRP(
406                                                                  ((AliConversionCuts*)fCutArray->At(fiCut))->IsHeavyIon(),
407                                                                  ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity(),
408                                                                  ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetNumberOfBGEvents());
409             fBGHandler[iCut] = NULL;
410          }
411       }
412    }
413 }
414 //________________________________________________________________________
415 void AliAnalysisTaskGammaConvV1::UserCreateOutputObjects()
416 {
417
418    // Create histograms
419    if(fOutputContainer != NULL){
420       delete fOutputContainer;
421       fOutputContainer = NULL;
422    }
423    if(fOutputContainer == NULL){
424       fOutputContainer = new TList();
425       fOutputContainer->SetOwner(kTRUE);
426    }
427
428    // Array of current cut's gammas
429    fGammaCandidates = new TList();
430
431    fCutFolder = new TList*[fnCuts];
432    fESDList = new TList*[fnCuts];
433    fBackList = new TList*[fnCuts];
434    fMotherList = new TList*[fnCuts];
435    hNEvents = new TH1I*[fnCuts];
436    hNGoodESDTracks = new TH1I*[fnCuts];
437    hNGammaCandidates = new TH1I*[fnCuts];
438    hNV0Tracks = new TH1I*[fnCuts];
439    hEtaShift = new TProfile*[fnCuts];
440    hESDConvGammaPt = new TH1F*[fnCuts];
441
442    if (fDoPhotonQA == 2){
443       fPhotonDCAList = new TList*[fnCuts];
444       tESDConvGammaPtDcazCat = new TTree*[fnCuts];
445    }
446    if (fDoPhotonQA > 0){
447       hESDConvGammaR = new TH1F*[fnCuts];
448       hESDConvGammaEta = new TH1F*[fnCuts];
449    } 
450    
451    if(fDoMesonAnalysis){
452       hESDMotherInvMassPt = new TH2F*[fnCuts];
453       hESDMotherBackInvMassPt = new TH2F*[fnCuts];
454       hESDMotherInvMassEalpha = new TH2F*[fnCuts];
455       if (fDoMesonQA == 2){
456          fMesonDCAList = new TList*[fnCuts];
457          tESDMesonsInvMassPtDcazMinDcazMaxFlag = new TTree*[fnCuts];
458       }
459       if (fDoMesonQA > 0){
460          hESDMotherPi0PtY =  new TH2F*[fnCuts];
461          hESDMotherEtaPtY =  new TH2F*[fnCuts];
462          hESDMotherPi0PtAlpha =  new TH2F*[fnCuts];
463          hESDMotherEtaPtAlpha =  new TH2F*[fnCuts];
464          hESDMotherPi0PtOpenAngle =  new TH2F*[fnCuts];
465          hESDMotherEtaPtOpenAngle =  new TH2F*[fnCuts];
466       }   
467    }
468
469    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
470
471       TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
472       TString cutstringMeson = "NoMesonCut";
473       if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
474
475       fCutFolder[iCut] = new TList();
476       fCutFolder[iCut]->SetName(Form("Cut Number %s_%s",cutstring.Data(),cutstringMeson.Data()));
477       fCutFolder[iCut]->SetOwner(kTRUE);
478       fOutputContainer->Add(fCutFolder[iCut]);
479       fESDList[iCut] = new TList();
480       fESDList[iCut]->SetName(Form("%s_%s ESD histograms",cutstring.Data(),cutstringMeson.Data()));
481       fESDList[iCut]->SetOwner(kTRUE);
482       fCutFolder[iCut]->Add(fESDList[iCut]);
483
484       hNEvents[iCut] = new TH1I("NEvents","NEvents",9,-0.5,8.5);
485       hNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
486       hNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
487       hNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
488       if (((AliConversionCuts*)fCutArray->At(iCut))->IsSpecialTrigger() == 4 ){
489          TString TriggerNames = "Not Trigger: ";
490          TriggerNames = TriggerNames+ ( (AliConversionCuts*)fCutArray->At(iCut))->GetSpecialTriggerName();
491          hNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
492       } else {
493          hNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
494       }
495       hNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
496       hNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
497       hNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
498       hNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
499       hNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
500       fESDList[iCut]->Add(hNEvents[iCut]);
501       if(fIsHeavyIon) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",4000,0,4000);
502       else hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
503       fESDList[iCut]->Add(hNGoodESDTracks[iCut]);
504       if(fIsHeavyIon) hNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",100,0,100);
505       else hNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
506       fESDList[iCut]->Add(hNGammaCandidates[iCut]);
507       if(fIsHeavyIon) hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",30000,0,30000);
508       else hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2000,0,2000);
509       fESDList[iCut]->Add(hNV0Tracks[iCut]);
510       hEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
511       fESDList[iCut]->Add(hEtaShift[iCut]);
512       hESDConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
513       fESDList[iCut]->Add(hESDConvGammaPt[iCut]);
514
515       if (fDoPhotonQA == 2){
516          fPhotonDCAList[iCut] = new TList();
517          fPhotonDCAList[iCut]->SetName(Form("%s_%s Photon DCA tree",cutstring.Data(),cutstringMeson.Data()));
518          fPhotonDCAList[iCut]->SetOwner(kTRUE);
519          fCutFolder[iCut]->Add(fPhotonDCAList[iCut]);
520             
521         tESDConvGammaPtDcazCat[iCut] = new TTree("ESD_ConvGamma_Pt_Dcaz_R_Eta","ESD_ConvGamma_Pt_Dcaz_R_Eta_Cat");   
522          tESDConvGammaPtDcazCat[iCut]->Branch("Pt",&fPtGamma,"fPtGamma/F");
523          tESDConvGammaPtDcazCat[iCut]->Branch("DcaZPhoton",&fDCAzPhoton,"fDCAzPhoton/F");
524 //          tESDConvGammaPtDcazCat[iCut]->Branch("R",&fRConvPhoton,"fRConvPhoton/F");
525 //          tESDConvGammaPtDcazCat[iCut]->Branch("Eta",&fEtaPhoton,"fEtaPhoton/F");
526          
527          tESDConvGammaPtDcazCat[iCut]->Branch("cat",&iCatPhoton,"iCatPhoton/b");
528          if(fIsMC){
529             tESDConvGammaPtDcazCat[iCut]->Branch("photonMCInfo",&iPhotonMCInfo,"iPhotonMCInfo/b");
530          }
531          fPhotonDCAList[iCut]->Add(tESDConvGammaPtDcazCat[iCut]);
532       }
533
534       if (fDoPhotonQA > 0){
535          hESDConvGammaR[iCut] = new TH1F("ESD_ConvGamma_R","ESD_ConvGamma_R",800,0,200);
536          fESDList[iCut]->Add(hESDConvGammaR[iCut]);
537          hESDConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",2000,-2,2);
538          fESDList[iCut]->Add(hESDConvGammaEta[iCut]);
539       }
540       
541       if(fDoMesonAnalysis){
542          hESDMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",800,0,0.8,250,0,25);
543          fESDList[iCut]->Add(hESDMotherInvMassPt[iCut]);
544          hESDMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",800,0,0.8,250,0,25);
545          fESDList[iCut]->Add(hESDMotherBackInvMassPt[iCut]);
546          hESDMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",800,0,0.8,250,0,25);
547          fESDList[iCut]->Add(hESDMotherInvMassEalpha[iCut]);
548          if (fDoMesonQA == 2){    
549             fMesonDCAList[iCut] = new TList();
550             fMesonDCAList[iCut]->SetName(Form("%s_%s Meson DCA tree",cutstring.Data(),cutstringMeson.Data()));
551             fMesonDCAList[iCut]->SetOwner(kTRUE);
552             fCutFolder[iCut]->Add(fMesonDCAList[iCut]);
553             
554             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut] = new TTree("ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag","ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag");   
555             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("InvMass",&fInvMass,"fInvMass/F");
556             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("Pt",&fPt,"fPt/F");
557             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMin",&fDCAzGammaMin,"fDCAzGammaMin/F");
558             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMax",&fDCAzGammaMax,"fDCAzGammaMax/F");
559             tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("kind",&iFlag,"iFlag/b");
560             if(fIsMC){
561                tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("mesonMCInfo",&iMesonMCInfo,"iMesonMCInfo/b");
562             }
563             fMesonDCAList[iCut]->Add(tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]);
564               
565          }
566          if (fDoMesonQA > 0 ){
567             hESDMotherPi0PtY[iCut] = new TH2F("ESD_MotherPi0_Pt_Y","ESD_MotherPi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);            
568             SetLogBinningXTH2(hESDMotherPi0PtY[iCut]);
569             fESDList[iCut]->Add(hESDMotherPi0PtY[iCut]);
570             hESDMotherEtaPtY[iCut] = new TH2F("ESD_MotherEta_Pt_Y","ESD_MotherEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
571             SetLogBinningXTH2(hESDMotherEtaPtY[iCut]);
572             fESDList[iCut]->Add(hESDMotherEtaPtY[iCut]);
573             hESDMotherPi0PtAlpha[iCut] = new TH2F("ESD_MotherPi0_Pt_Alpha","ESD_MotherPi0_Pt_Alpha",150,0.03,15.,100,0,1);            
574             SetLogBinningXTH2(hESDMotherPi0PtAlpha[iCut]);
575             fESDList[iCut]->Add(hESDMotherPi0PtAlpha[iCut]);
576             hESDMotherEtaPtAlpha[iCut] = new TH2F("ESD_MotherEta_Pt_Alpha","ESD_MotherEta_Pt_Alpha",150,0.03,15.,100,0,1);
577             SetLogBinningXTH2(hESDMotherEtaPtAlpha[iCut]);
578             fESDList[iCut]->Add(hESDMotherEtaPtAlpha[iCut]);
579             hESDMotherPi0PtOpenAngle[iCut] = new TH2F("ESD_MotherPi0_Pt_OpenAngle","ESD_MotherPi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());            
580             SetLogBinningXTH2(hESDMotherPi0PtOpenAngle[iCut]);
581             fESDList[iCut]->Add(hESDMotherPi0PtOpenAngle[iCut]);
582             hESDMotherEtaPtOpenAngle[iCut] = new TH2F("ESD_MotherEta_Pt_OpenAngle","ESD_MotherEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());            
583             SetLogBinningXTH2(hESDMotherEtaPtOpenAngle[iCut]);
584             fESDList[iCut]->Add(hESDMotherEtaPtOpenAngle[iCut]);
585          }
586          
587             
588       }
589
590
591    }
592    if(fDoMesonAnalysis){
593       InitBack(); // Init Background Handler
594    }
595
596    if(fIsMC){
597       // MC Histogramms
598       fMCList = new TList*[fnCuts];
599       // True Histogramms
600       fTrueList = new TList*[fnCuts];
601       // Selected Header List
602       fHeaderNameList = new TList*[fnCuts];
603       hMCHeaders = new TH1I*[fnCuts];
604       hMCAllGammaPt = new TH1F*[fnCuts];
605       hMCDecayGammaPi0Pt = new TH1F*[fnCuts];
606       hMCDecayGammaRhoPt = new TH1F*[fnCuts];
607       hMCDecayGammaEtaPt = new TH1F*[fnCuts];
608       hMCDecayGammaOmegaPt = new TH1F*[fnCuts];
609       hMCDecayGammaEtapPt = new TH1F*[fnCuts];
610       hMCDecayGammaPhiPt = new TH1F*[fnCuts];
611       hMCDecayGammaSigmaPt = new TH1F*[fnCuts];
612       hMCConvGammaPt = new TH1F*[fnCuts];
613       hESDTrueConvGammaPt = new TH1F*[fnCuts];
614
615       hESDCombinatorialPt = new TH2F*[fnCuts];
616       hESDTruePrimaryConvGammaPt = new TH1F*[fnCuts];
617       hESDTruePrimaryConvGammaESDPtMCPt = new TH2F*[fnCuts];
618       hESDTrueSecondaryConvGammaPt = new TH1F*[fnCuts];
619       hESDTrueSecondaryConvGammaFromXFromK0sPt = new TH1F*[fnCuts];
620       hESDTrueSecondaryConvGammaFromXFromLambdaPt = new TH1F*[fnCuts];
621
622       hESDTrueDalitzPsiPairDeltaPhi= new TH2F*[fnCuts];
623       hESDTrueGammaPsiPairDeltaPhi= new TH2F*[fnCuts];
624
625       if (fDoPhotonQA > 0){
626          hMCConvGammaR = new TH1F*[fnCuts];
627          hMCConvGammaEta = new TH1F*[fnCuts];
628       }
629
630       if(fDoMesonAnalysis){
631          hMCPi0Pt = new TH1F*[fnCuts];
632          hMCPi0WOWeightPt = new TH1F*[fnCuts];
633          hMCEtaPt = new TH1F*[fnCuts];
634          hMCEtaWOWeightPt = new TH1F*[fnCuts];
635          hMCPi0InAccPt = new TH1F*[fnCuts];
636          hMCEtaInAccPt = new TH1F*[fnCuts];
637
638          hESDTrueMotherInvMassPt = new TH2F*[fnCuts];
639          hESDTruePrimaryMotherInvMassPt = new TH2F*[fnCuts];
640          hESDTruePrimaryMotherW0WeightingInvMassPt = new TH2F*[fnCuts];
641          pESDTruePrimaryMotherWeightsInvMassPt = new TProfile2D*[fnCuts];
642          hESDTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
643          hESDTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
644          hESDTrueSecondaryMotherFromEtaInvMassPt = new TH2F*[fnCuts];
645          hESDTrueMotherDalitzInvMassPt = new TH2F*[fnCuts];
646          if (fDoMesonQA > 0){
647             hMCPi0PtY = new TH2F*[fnCuts];
648             hMCEtaPtY = new TH2F*[fnCuts];
649             hMCK0sPt = new TH1F*[fnCuts];
650             hMCK0sWOWeightPt = new TH1F*[fnCuts];
651             hMCK0sPtY = new TH2F*[fnCuts];
652             hESDTruePrimaryPi0MCPtResolPt = new TH2F*[fnCuts];
653             hESDTruePrimaryEtaMCPtResolPt = new TH2F*[fnCuts];
654             hESDTrueK0sWithPi0DaughterMCPt = new TH1F*[fnCuts];
655             hESDTrueEtaWithPi0DaughterMCPt = new TH1F*[fnCuts];
656             hESDTrueBckGGInvMassPt = new TH2F*[fnCuts];
657             hESDTrueBckContInvMassPt = new TH2F*[fnCuts];
658             hESDTruePi0PtY = new TH2F*[fnCuts];
659             hESDTrueEtaPtY = new TH2F*[fnCuts];
660             hESDTruePi0PtAlpha = new TH2F*[fnCuts];
661             hESDTrueEtaPtAlpha = new TH2F*[fnCuts];
662             hESDTruePi0PtOpenAngle = new TH2F*[fnCuts];
663             hESDTrueEtaPtOpenAngle = new TH2F*[fnCuts];
664          }
665       }
666
667       for(Int_t iCut = 0; iCut<fnCuts;iCut++){
668          TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
669          TString cutstringMeson = "NoMesonCut";
670          if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
671
672          fMCList[iCut] = new TList();
673          fMCList[iCut]->SetName(Form("%s_%s MC histograms",cutstring.Data(),cutstringMeson.Data()));
674          fMCList[iCut]->SetOwner(kTRUE);
675          fCutFolder[iCut]->Add(fMCList[iCut]);
676          hMCHeaders[iCut] = new TH1I("MC_Headers","MC_Headers",20,0,20);
677          fMCList[iCut]->Add(hMCHeaders[iCut]);
678          hMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
679          fMCList[iCut]->Add(hMCAllGammaPt[iCut]);
680          hMCDecayGammaPi0Pt[iCut] = new TH1F("MC_DecayGammaPi0_Pt","MC_DecayGammaPi0_Pt",250,0,25);
681          fMCList[iCut]->Add(hMCDecayGammaPi0Pt[iCut]);
682          hMCDecayGammaRhoPt[iCut] = new TH1F("MC_DecayGammaRho_Pt","MC_DecayGammaRho_Pt",250,0,25);
683          fMCList[iCut]->Add(hMCDecayGammaRhoPt[iCut]);
684          hMCDecayGammaEtaPt[iCut] = new TH1F("MC_DecayGammaEta_Pt","MC_DecayGammaEta_Pt",250,0,25);
685          fMCList[iCut]->Add(hMCDecayGammaEtaPt[iCut]);
686          hMCDecayGammaOmegaPt[iCut] = new TH1F("MC_DecayGammaOmega_Pt","MC_DecayGammaOmmega_Pt",250,0,25);
687          fMCList[iCut]->Add(hMCDecayGammaOmegaPt[iCut]);
688          hMCDecayGammaEtapPt[iCut] = new TH1F("MC_DecayGammaEtap_Pt","MC_DecayGammaEtap_Pt",250,0,25);
689          fMCList[iCut]->Add(hMCDecayGammaEtapPt[iCut]);
690          hMCDecayGammaPhiPt[iCut] = new TH1F("MC_DecayGammaPhi_Pt","MC_DecayGammaPhi_Pt",250,0,25);
691          fMCList[iCut]->Add(hMCDecayGammaPhiPt[iCut]);
692          hMCDecayGammaSigmaPt[iCut] = new TH1F("MC_DecayGammaSigma_Pt","MC_DecayGammaSigma_Pt",250,0,25);
693          fMCList[iCut]->Add(hMCDecayGammaSigmaPt[iCut]);
694          hMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
695          fMCList[iCut]->Add(hMCConvGammaPt[iCut]);
696       
697          if (fDoPhotonQA > 0){
698             hMCConvGammaR[iCut] = new TH1F("MC_ConvGamma_R","MC_ConvGamma_R",800,0,200);
699             fMCList[iCut]->Add(hMCConvGammaR[iCut]);
700             hMCConvGammaEta[iCut] = new TH1F("MC_ConvGamma_Eta","MC_ConvGamma_Eta",100,-4,4);
701             fMCList[iCut]->Add(hMCConvGammaEta[iCut]);
702          }
703
704          if(fDoMesonAnalysis){
705             hMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
706             hMCPi0Pt[iCut]->Sumw2();
707             fMCList[iCut]->Add(hMCPi0Pt[iCut]);
708             hMCPi0WOWeightPt[iCut] = new TH1F("MC_Pi0_WOWeights_Pt","MC_Pi0_WOWeights_Pt",250,0,25);
709             hMCPi0WOWeightPt[iCut]->Sumw2();
710             fMCList[iCut]->Add(hMCPi0WOWeightPt[iCut]);
711             
712             hMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
713             hMCEtaPt[iCut]->Sumw2();
714             fMCList[iCut]->Add(hMCEtaPt[iCut]);
715             hMCEtaWOWeightPt[iCut] = new TH1F("MC_Eta_WOWeights_Pt","MC_Eta_WOWeights_Pt",250,0,25);
716             hMCEtaWOWeightPt[iCut]->Sumw2();
717             fMCList[iCut]->Add(hMCEtaWOWeightPt[iCut]);
718             hMCPi0InAccPt[iCut] = new TH1F("MC_Pi0InAcc_Pt","MC_Pi0InAcc_Pt",250,0,25);
719             hMCPi0InAccPt[iCut]->Sumw2();
720             fMCList[iCut]->Add(hMCPi0InAccPt[iCut]);
721             hMCEtaInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
722             hMCEtaInAccPt[iCut]->Sumw2();
723             fMCList[iCut]->Add(hMCEtaInAccPt[iCut]);
724             if (fDoMesonQA > 0){
725                hMCPi0PtY[iCut] = new TH2F("MC_Pi0_Pt_Y","MC_Pi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
726                hMCPi0PtY[iCut]->Sumw2();
727                SetLogBinningXTH2(hMCPi0PtY[iCut]);
728                fMCList[iCut]->Add(hMCPi0PtY[iCut]);
729                hMCEtaPtY[iCut] = new TH2F("MC_Eta_Pt_Y","MC_Eta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
730                hMCEtaPtY[iCut]->Sumw2();
731                SetLogBinningXTH2(hMCEtaPtY[iCut]);
732                fMCList[iCut]->Add(hMCEtaPtY[iCut]);
733                hMCK0sPt[iCut] = new TH1F("MC_K0s_Pt","MC_K0s_Pt",150,0,15);
734                hMCK0sPt[iCut]->Sumw2();
735                fMCList[iCut]->Add(hMCK0sPt[iCut]);
736                hMCK0sWOWeightPt[iCut] = new TH1F("MC_K0s_WOWeights_Pt","MC_K0s_WOWeights_Pt",150,0,15);
737                hMCK0sWOWeightPt[iCut]->Sumw2();
738                fMCList[iCut]->Add(hMCK0sWOWeightPt[iCut]);
739                hMCK0sPtY[iCut] = new TH2F("MC_K0s_Pt_Y","MC_K0s_Pt_Y",150,0.03,15.,150,-1.5,1.5);
740                hMCK0sPtY[iCut]->Sumw2();
741                SetLogBinningXTH2(hMCK0sPtY[iCut]);
742                fMCList[iCut]->Add(hMCK0sPtY[iCut]);
743                
744             }
745
746          }
747          fTrueList[iCut] = new TList();
748          fTrueList[iCut]->SetName(Form("%s_%s True histograms",cutstring.Data(),cutstringMeson.Data()));
749          fTrueList[iCut]->SetOwner(kTRUE);
750          fCutFolder[iCut]->Add(fTrueList[iCut]);
751
752          hESDTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
753          fTrueList[iCut]->Add(hESDTrueConvGammaPt[iCut]);
754
755          hESDCombinatorialPt[iCut] = new TH2F("ESD_TrueCombinatorial_Pt","ESD_TrueCombinatorial_Pt",250,0,25,16,-0.5,15.5);
756          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 1,"Elec+Elec");
757          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 2,"Elec+Pion");
758          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 3,"Elec+Kaon");
759          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 4,"Elec+Proton");
760          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 5,"Elec+Muon");
761          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 6,"Pion+Pion");
762          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 7,"Pion+Kaon");
763          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 8,"Pion+Proton");
764          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel( 9,"Pion+Muon");
765          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(10,"Kaon+Kaon");
766          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(11,"Kaon+Proton");
767          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(12,"Kaon+Muon");
768          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(13,"Proton+Proton");
769          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(14,"Proton+Muon");
770          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(15,"Muon+Muon");
771          hESDCombinatorialPt[iCut]->GetYaxis()->SetBinLabel(16,"Rest");
772          fTrueList[iCut]->Add(hESDCombinatorialPt[iCut]);
773          hESDTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt","ESD_TruePrimaryConvGamma_Pt",250,0,25);
774          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaPt[iCut]);
775          hESDTrueSecondaryConvGammaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGamma_Pt","ESD_TrueSecondaryConvGamma_Pt",250,0,25);
776          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaPt[iCut]);
777
778          hESDTrueSecondaryConvGammaFromXFromK0sPt[iCut]
779             = new TH1F("ESD_TrueSecondaryConvGammaFromXFromK0s_Pt", "ESD_TrueSecondaryConvGammaFromXFromK0s_Pt",250,0,25);
780          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaFromXFromK0sPt[iCut]);
781          hESDTrueSecondaryConvGammaFromXFromLambdaPt[iCut]
782             = new TH1F("ESD_TrueSecondaryConvGammaFromXFromLambda_Pt", "ESD_TrueSecondaryConvGammaFromXFromLambda_Pt",250,0,25);
783          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaFromXFromLambdaPt[iCut]);
784
785          hESDTrueDalitzPsiPairDeltaPhi[iCut]
786             = new TH2F("ESD_TrueDalitzPsiPairDeltaPhi_Pt", "ESD_TrueDalitzPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
787          fTrueList[iCut]->Add(hESDTrueDalitzPsiPairDeltaPhi[iCut]);
788          
789          hESDTrueGammaPsiPairDeltaPhi[iCut]
790             = new TH2F("ESD_TrueGammaPsiPairDeltaPhi_Pt", "ESD_TrueGammaPsiPairDeltaPhi_Pt",100,-0.5,2,100,-0.5,0.5);
791          fTrueList[iCut]->Add(hESDTrueGammaPsiPairDeltaPhi[iCut]);
792
793          
794          hESDTruePrimaryConvGammaESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryConvGammaESD_PtMCPt", "ESD_TruePrimaryConvGammaESD_PtMCPt",250,0,25,250,0,25);
795          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaESDPtMCPt[iCut]);
796          
797          if(fDoMesonAnalysis){
798             hESDTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",800,0,0.8,250,0,25);
799             fTrueList[iCut]->Add(hESDTrueMotherInvMassPt[iCut]);
800             hESDTruePrimaryMotherInvMassPt[iCut]
801                = new TH2F("ESD_TruePrimaryMother_InvMass_Pt", "ESD_TruePrimaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
802             hESDTruePrimaryMotherInvMassPt[iCut]->Sumw2();
803             fTrueList[iCut]->Add(hESDTruePrimaryMotherInvMassPt[iCut]);
804             hESDTruePrimaryMotherW0WeightingInvMassPt[iCut]
805                = new TH2F("ESD_TruePrimaryMotherW0Weights_InvMass_Pt", "ESD_TruePrimaryMotherW0Weights_InvMass_Pt", 800,0,0.8,250,0,25);
806             hESDTruePrimaryMotherW0WeightingInvMassPt[iCut]->Sumw2();
807             fTrueList[iCut]->Add(hESDTruePrimaryMotherW0WeightingInvMassPt[iCut]);
808             pESDTruePrimaryMotherWeightsInvMassPt[iCut]
809                = new TProfile2D("ESD_TruePrimaryMotherWeights_InvMass_Pt", "ESD_TruePrimaryMotherWeights_InvMass_Pt", 800,0,0.8,250,0,25);
810             pESDTruePrimaryMotherWeightsInvMassPt[iCut]->Sumw2();
811             fTrueList[iCut]->Add(pESDTruePrimaryMotherWeightsInvMassPt[iCut]);
812             hESDTrueSecondaryMotherInvMassPt[iCut]
813                = new TH2F("ESD_TrueSecondaryMother_InvMass_Pt", "ESD_TrueSecondaryMother_InvMass_Pt", 800,0,0.8,250,0,25);
814             hESDTrueSecondaryMotherInvMassPt[iCut]->Sumw2();
815             fTrueList[iCut]->Add(hESDTrueSecondaryMotherInvMassPt[iCut]);
816             hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]
817                = new TH2F("ESD_TrueSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueSecondaryMotherFromK0s_InvMass_Pt",800,0,0.8,250,0,25);
818             hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]->Sumw2();
819             fTrueList[iCut]->Add(hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]);
820             hESDTrueSecondaryMotherFromEtaInvMassPt[iCut]
821                = new TH2F("ESD_TrueSecondaryMotherFromEta_InvMass_Pt","ESD_TrueSecondaryMotherFromEta_InvMass_Pt",800,0,0.8,250,0,25);
822             fTrueList[iCut]->Add(hESDTrueSecondaryMotherFromEtaInvMassPt[iCut]);
823             hESDTrueMotherDalitzInvMassPt[iCut] = new TH2F("ESD_TrueDalitz_InvMass_Pt","ESD_TrueDalitz_InvMass_Pt",800,0,0.8,250,0,25);
824             fTrueList[iCut]->Add(hESDTrueMotherDalitzInvMassPt[iCut]);         
825             if (fDoMesonQA > 0){
826                hESDTruePrimaryPi0MCPtResolPt[iCut] = new TH2F("ESD_TruePrimaryPi0_MCPt_ResolPt","ESD_TruePrimaryPi0_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
827                hESDTruePrimaryPi0MCPtResolPt[iCut]->Sumw2();
828                SetLogBinningXTH2(hESDTruePrimaryPi0MCPtResolPt[iCut]);
829                fTrueList[iCut]->Add(hESDTruePrimaryPi0MCPtResolPt[iCut]);
830                hESDTruePrimaryEtaMCPtResolPt[iCut]  = new TH2F("ESD_TruePrimaryEta_MCPt_ResolPt","ESD_TruePrimaryEta_ResolPt_MCPt",500,0.03,25,1000,-1.,1.);
831                hESDTruePrimaryEtaMCPtResolPt[iCut]->Sumw2();
832                SetLogBinningXTH2(hESDTruePrimaryEtaMCPtResolPt[iCut]);
833                fTrueList[iCut]->Add(hESDTruePrimaryEtaMCPtResolPt[iCut]);
834                hESDTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt","ESD_TrueBckGG_InvMass_Pt",800,0,0.8,250,0,25);
835                fTrueList[iCut]->Add(hESDTrueBckGGInvMassPt[iCut]);
836                hESDTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt","ESD_TrueBckCont_InvMass_Pt",800,0,0.8,250,0,25);
837                fTrueList[iCut]->Add(hESDTrueBckContInvMassPt[iCut]);
838                hESDTrueK0sWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueK0sWithPi0Daughter_MCPt","ESD_TrueK0sWithPi0Daughter_MCPt",250,0,25);
839                fTrueList[iCut]->Add(hESDTrueK0sWithPi0DaughterMCPt[iCut]);
840                hESDTrueEtaWithPi0DaughterMCPt[iCut] = new TH1F("ESD_TrueEtaWithPi0Daughter_MCPt","ESD_TrueEtaWithPi0Daughter_MCPt",250,0,25);
841                fTrueList[iCut]->Add(hESDTrueEtaWithPi0DaughterMCPt[iCut]);
842                hESDTruePi0PtY[iCut] = new TH2F("ESD_TruePi0_Pt_Y","ESD_TruePi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
843                SetLogBinningXTH2(hESDTruePi0PtY[iCut]);
844                fTrueList[iCut]->Add(hESDTruePi0PtY[iCut]);
845                hESDTrueEtaPtY[iCut] = new TH2F("ESD_TrueEta_Pt_Y","ESD_TrueEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
846                SetLogBinningXTH2(hESDTrueEtaPtY[iCut]);
847                fTrueList[iCut]->Add(hESDTrueEtaPtY[iCut]);
848                hESDTruePi0PtAlpha[iCut] = new TH2F("ESD_TruePi0_Pt_Alpha","ESD_TruePi0_Pt_Alpha",150,0.03,15.,100,0,1);
849                SetLogBinningXTH2(hESDTruePi0PtAlpha[iCut]);
850                fTrueList[iCut]->Add(hESDTruePi0PtAlpha[iCut]);
851                hESDTrueEtaPtAlpha[iCut] = new TH2F("ESD_TrueEta_Pt_Alpha","ESD_TrueEta_Pt_Alpha",150,0.03,15.,100,0,1);
852                SetLogBinningXTH2(hESDTrueEtaPtAlpha[iCut]);
853                fTrueList[iCut]->Add(hESDTrueEtaPtAlpha[iCut]);
854                
855                hESDTruePi0PtOpenAngle[iCut] = new TH2F("ESD_TruePi0_Pt_OpenAngle","ESD_TruePi0_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());            
856                SetLogBinningXTH2(hESDTruePi0PtOpenAngle[iCut]);
857                fTrueList[iCut]->Add(hESDTruePi0PtOpenAngle[iCut]);
858                hESDTrueEtaPtOpenAngle[iCut] = new TH2F("ESD_TrueEta_Pt_OpenAngle","ESD_TrueEta_Pt_OpenAngle",150,0.03,15.,200,0,2*TMath::Pi());            
859                SetLogBinningXTH2(hESDTrueEtaPtOpenAngle[iCut]);
860                fTrueList[iCut]->Add(hESDTrueEtaPtOpenAngle[iCut]);
861                
862             }
863          }
864       }
865    }
866
867    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
868    if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
869
870    if(fV0Reader)
871       if((AliConversionCuts*)fV0Reader->GetConversionCuts())
872          if(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
873             fOutputContainer->Add(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
874
875    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
876       if(!((AliConversionCuts*)fCutArray->At(iCut))) continue;
877       if(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms()){
878          fCutFolder[iCut]->Add(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms());
879       }
880       if(fDoMesonAnalysis){
881          if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue;
882          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
883             fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
884          }
885       }
886    }
887    PostData(1, fOutputContainer);
888 }
889 //_____________________________________________________________________________
890 Bool_t AliAnalysisTaskGammaConvV1::Notify()
891 {
892    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
893       if(!((AliConversionCuts*)fCutArray->At(iCut))->GetDoEtaShift()){
894          hEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
895          continue; // No Eta Shift requested, continue
896       }
897       if(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
898          ((AliConversionCuts*)fCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
899           hEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
900          ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once   
901          continue;
902       }
903       else{
904          printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n",
905                 (((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber()).Data(),((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift());
906           hEtaShift[iCut]->Fill(0.,(((AliConversionCuts*)fCutArray->At(iCut))->GetEtaShift()));
907          ((AliConversionCuts*)fCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once   
908       }
909    }
910    
911    return kTRUE;
912 }
913 //_____________________________________________________________________________
914 void AliAnalysisTaskGammaConvV1::UserExec(Option_t *)
915 {
916    //
917    // Called for each event
918    //
919    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
920    if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
921       for(Int_t iCut = 0; iCut<fnCuts; iCut++){
922          hNEvents[iCut]->Fill(eventQuality);
923       }
924       return;
925    }
926
927    if(fIsMC) fMCEvent = MCEvent();
928    if(fMCEvent == NULL) fIsMC = kFALSE;
929    
930    fInputEvent = InputEvent();
931
932    if(fIsMC && fInputEvent->IsA()==AliESDEvent::Class()){
933       fMCStack = fMCEvent->Stack();
934       if(fMCStack == NULL) fIsMC = kFALSE;
935    }
936
937    fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
938    
939    // ------------------- BeginEvent ----------------------------
940
941    AliEventplane *EventPlane = fInputEvent->GetEventplane();
942    if(fIsHeavyIon)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2);
943    else fEventPlaneAngle=0.0;
944    
945    if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
946       RelabelAODPhotonCandidates(kTRUE);    // In case of AODMC relabeling MC
947       fV0Reader->RelabelAODs(kTRUE);
948    }
949    for(Int_t iCut = 0; iCut<fnCuts; iCut++){
950       fiCut = iCut;
951       Int_t eventNotAccepted =
952          ((AliConversionCuts*)fCutArray->At(iCut))
953          ->IsEventAcceptedByConversionCut(fV0Reader->GetConversionCuts(),fInputEvent,fMCEvent,fIsHeavyIon);
954       if(eventNotAccepted){
955          // cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
956          hNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
957          continue;
958       }
959
960       if(eventQuality != 0){// Event Not Accepted
961          // cout << "event rejected due to: " <<eventQuality << endl;
962          hNEvents[iCut]->Fill(eventQuality);
963          continue;
964       }
965
966       hNEvents[iCut]->Fill(eventQuality); // Should be 0 here
967       hNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks());
968       if(((AliConversionCuts*)fCutArray->At(iCut))->IsHeavyIon() == 2)  hNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A());
969       else hNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C());
970
971       if(fIsMC){
972          // Process MC Particle
973          if(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection() != 0){
974             if(fInputEvent->IsA()==AliESDEvent::Class()){
975                ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(),
976                                                                                   ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(),
977                                                                                   fMCEvent);
978             }
979             else if(fInputEvent->IsA()==AliAODEvent::Class()){
980                ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(),
981                                                                                   ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(),
982                                                                                   fInputEvent);
983             }
984
985             if(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader()){
986                for(Int_t i = 0;i<(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){
987                   TString nameBin= hMCHeaders[iCut]->GetXaxis()->GetBinLabel(i+1);
988                   if (nameBin.CompareTo("")== 0){
989                      TString nameHeader = ((TObjString*)((TList*)((AliConversionCuts*)fCutArray->At(iCut))
990                                                          ->GetAcceptedHeader())->At(i))->GetString();
991                      hMCHeaders[iCut]->GetXaxis()->SetBinLabel(i+1,nameHeader.Data());
992                   }
993                }
994             }
995          }
996       }
997       if(fIsMC){
998          if(fInputEvent->IsA()==AliESDEvent::Class())
999             ProcessMCParticles();
1000          if(fInputEvent->IsA()==AliAODEvent::Class())
1001             ProcessAODMCParticles();
1002       }
1003
1004       ProcessPhotonCandidates(); // Process this cuts gammas
1005
1006       hNGammaCandidates[iCut]->Fill(fGammaCandidates->GetEntries());
1007       if(fDoMesonAnalysis){ // Meson Analysis
1008          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1009             fUnsmearedPx = new Double_t[fGammaCandidates->GetEntries()]; // Store unsmeared Momenta
1010             fUnsmearedPy = new Double_t[fGammaCandidates->GetEntries()];
1011             fUnsmearedPz = new Double_t[fGammaCandidates->GetEntries()];
1012             fUnsmearedE =  new Double_t[fGammaCandidates->GetEntries()];
1013
1014             for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1015                fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Px();
1016                fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Py();
1017                fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->Pz();
1018                fUnsmearedE[gamma] =  ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->E();
1019                ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(gamma)));
1020             }
1021          }
1022
1023          CalculatePi0Candidates(); // Combine Gammas
1024          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->DoBGCalculation()){
1025             if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->BackgroundHandlerType() == 0){
1026                CalculateBackground(); // Combinatorial Background
1027                UpdateEventByEventData(); // Store Event for mixed Events
1028             }
1029             else{
1030                CalculateBackgroundRP(); // Combinatorial Background
1031                fBGHandlerRP[iCut]->AddEvent(fGammaCandidates,fInputEvent); // Store Event for mixed Events
1032             }
1033          }
1034          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fIsMC){
1035             for(Int_t gamma=0;gamma<fGammaCandidates->GetEntries();gamma++){ // Smear the AODPhotons in MC
1036                ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
1037                ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPy(fUnsmearedPy[gamma]);
1038                ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetPz(fUnsmearedPz[gamma]);
1039                ((AliAODConversionPhoton*)fGammaCandidates->At(gamma))->SetE(fUnsmearedE[gamma]);
1040             }
1041             delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
1042             delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
1043             delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
1044             delete[] fUnsmearedE;  fUnsmearedE  = 0x0;
1045          }
1046       }
1047       fGammaCandidates->Clear(); // delete this cuts good gammas
1048    }
1049
1050    if(fIsMC && fInputEvent->IsA()==AliAODEvent::Class() && !(fV0Reader->AreAODsRelabeled())){
1051       RelabelAODPhotonCandidates(kFALSE); // Back to ESDMC Label
1052       fV0Reader->RelabelAODs(kFALSE);
1053    }
1054
1055    PostData(1, fOutputContainer);
1056 }
1057 //________________________________________________________________________
1058 void AliAnalysisTaskGammaConvV1::ProcessPhotonCandidates()
1059 {
1060    Int_t nV0 = 0;
1061    TList *GammaCandidatesStepOne = new TList();
1062    TList *GammaCandidatesStepTwo = new TList();
1063    // Loop over Photon Candidates allocated by ReaderV1
1064    for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
1065       AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
1066       if(!PhotonCandidate) continue;
1067       fIsFromMBHeader = kTRUE;
1068       if(fIsMC && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1069          Int_t isPosFromMBHeader
1070             = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1071          if(isPosFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1072          Int_t isNegFromMBHeader
1073             = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1074          if(isNegFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1075
1076          if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1077       }
1078
1079       if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue;
1080       if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue;
1081       if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1082          !((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1083          fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
1084
1085          if(fIsFromMBHeader){
1086             hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1087             if (fDoPhotonQA > 0){
1088                hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1089                hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1090             }   
1091          }
1092          if(fIsMC){
1093             if(fInputEvent->IsA()==AliESDEvent::Class())
1094                ProcessTruePhotonCandidates(PhotonCandidate);
1095             if(fInputEvent->IsA()==AliAODEvent::Class())
1096                ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1097          }
1098          if (fIsFromMBHeader && fDoPhotonQA == 2){
1099             if (fIsHeavyIon && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1100                fPtGamma = PhotonCandidate->Pt();
1101                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1102                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1103                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1104                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1105                tESDConvGammaPtDcazCat[fiCut]->Fill();
1106             } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1107                fPtGamma = PhotonCandidate->Pt();
1108                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1109                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1110                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1111                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1112                tESDConvGammaPtDcazCat[fiCut]->Fill();
1113             }   
1114          }   
1115       } else if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
1116          ((AliConversionCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
1117          nV0++;
1118          GammaCandidatesStepOne->Add(PhotonCandidate);
1119       } else if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
1120               ((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
1121          GammaCandidatesStepTwo->Add(PhotonCandidate);
1122       }
1123    }
1124    if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
1125       for(Int_t i = 0;i<GammaCandidatesStepOne->GetEntries();i++){
1126          AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i);
1127          if(!PhotonCandidate) continue;
1128          fIsFromMBHeader = kTRUE;
1129          if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1130             Int_t isPosFromMBHeader
1131                = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1132             Int_t isNegFromMBHeader
1133                = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1134             if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1135          }
1136          if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue;
1137          if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
1138             fGammaCandidates->Add(PhotonCandidate);
1139             if(fIsFromMBHeader){
1140                hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1141                if (fDoPhotonQA > 0){
1142                   hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1143                   hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1144                }
1145             }
1146          }
1147          if(fIsMC){
1148             if(fInputEvent->IsA()==AliESDEvent::Class())
1149                ProcessTruePhotonCandidates(PhotonCandidate);
1150             if(fInputEvent->IsA()==AliAODEvent::Class())
1151                ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1152          } else GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
1153          
1154            if (fIsFromMBHeader && fDoPhotonQA == 2){
1155             if (fIsHeavyIon && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1156                fPtGamma = PhotonCandidate->Pt();
1157                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1158                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1159                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1160                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1161                tESDConvGammaPtDcazCat[fiCut]->Fill();
1162             } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1163                fPtGamma = PhotonCandidate->Pt();
1164                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1165                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1166                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1167                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1168                tESDConvGammaPtDcazCat[fiCut]->Fill();
1169             }
1170          }
1171       }
1172    }
1173    if(((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
1174       for(Int_t i = 0;i<GammaCandidatesStepTwo->GetEntries();i++){
1175          AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i);
1176          if(!PhotonCandidate) continue;
1177          fIsFromMBHeader = kTRUE;
1178          if(fMCStack && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1179             Int_t isPosFromMBHeader
1180                = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack, fInputEvent);
1181             Int_t isNegFromMBHeader
1182                = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack, fInputEvent);
1183             if( (isNegFromMBHeader+isPosFromMBHeader) != 4) fIsFromMBHeader = kFALSE;
1184          }
1185          if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue;
1186          fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList
1187          if(fIsFromMBHeader){
1188             hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
1189             if (fDoPhotonQA > 0){
1190                hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
1191                hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
1192             }
1193          }
1194          if(fIsMC){
1195             if(fInputEvent->IsA()==AliESDEvent::Class())
1196                ProcessTruePhotonCandidates(PhotonCandidate);
1197             if(fInputEvent->IsA()==AliAODEvent::Class())
1198                ProcessTruePhotonCandidatesAOD(PhotonCandidate);
1199          }
1200          if (fIsFromMBHeader){
1201            if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
1202                fPtGamma = PhotonCandidate->Pt();
1203                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1204                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1205                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1206                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1207                tESDConvGammaPtDcazCat[fiCut]->Fill();
1208             } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
1209                fPtGamma = PhotonCandidate->Pt();
1210                fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
1211                fRConvPhoton = PhotonCandidate->GetConversionRadius();
1212                fEtaPhoton = PhotonCandidate->GetPhotonEta();
1213                iCatPhoton = PhotonCandidate->GetPhotonQuality();
1214                tESDConvGammaPtDcazCat[fiCut]->Fill();
1215             }
1216          }
1217       }
1218    }
1219
1220    delete GammaCandidatesStepOne;
1221    GammaCandidatesStepOne = 0x0;
1222    delete GammaCandidatesStepTwo;
1223    GammaCandidatesStepTwo = 0x0;
1224
1225 }
1226 //________________________________________________________________________
1227 void AliAnalysisTaskGammaConvV1::ProcessTruePhotonCandidatesAOD(AliAODConversionPhoton *TruePhotonCandidate)
1228 {
1229
1230    Double_t magField = fInputEvent->GetMagneticField();
1231    if( magField  < 0.0 ){
1232       magField =  1.0;
1233    }
1234    else {
1235       magField =  -1.0;
1236    }
1237
1238    TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1239    AliAODMCParticle *posDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelPositive());
1240    AliAODMCParticle *negDaughter = (AliAODMCParticle*) AODMCTrackArray->At(TruePhotonCandidate->GetMCLabelNegative());
1241    iPhotonMCInfo = 0;
1242    
1243    if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1244    Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1245
1246    if(posDaughter->GetMother() != negDaughter->GetMother()){
1247       FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1248       iPhotonMCInfo = 1;
1249       return;
1250    }
1251    else if(posDaughter->GetMother() == -1){
1252       FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1253       iPhotonMCInfo = 1;
1254       return;
1255    }
1256
1257    if(pdgCode[0]!=11 || pdgCode[1]!=11){
1258       iPhotonMCInfo = 1;
1259       return; //One Particle is not a electron
1260    }
1261    if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()){
1262       iPhotonMCInfo = 1; 
1263       return; // Same Charge
1264    }
1265    
1266    AliAODMCParticle *Photon = (AliAODMCParticle*) AODMCTrackArray->At(posDaughter->GetMother());
1267    AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1268    AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1269    Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1270
1271    if(Photon->GetPdgCode() != 22){
1272       hESDTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1273       iPhotonMCInfo = 1;
1274       return; // Mother is no Photon
1275    }
1276    
1277    if(((posDaughter->GetMCProcessCode())) != 5 || ((negDaughter->GetMCProcessCode())) != 5){
1278       iPhotonMCInfo = 1;
1279       return;// check if the daughters come from a conversion 
1280    }   
1281    // STILL A BUG IN ALIROOT >>8 HAS TPO BE REMOVED AFTER FIX
1282    
1283    
1284    
1285    // True Photon
1286    if(fIsFromMBHeader)hESDTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1287    hESDTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());  
1288    if(Photon->IsPrimary()){ 
1289       // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1290       if(fIsFromMBHeader){
1291          iPhotonMCInfo = 6;
1292          hESDTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1293          hESDTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1294       }
1295       // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1296    }
1297    else{
1298       if(fIsFromMBHeader){
1299          hESDTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1300          iPhotonMCInfo = 2;
1301          if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1302             ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 3122){
1303             iPhotonMCInfo = 5;
1304             hESDTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1305          }
1306          if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1307             ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 310){
1308             iPhotonMCInfo = 4;
1309             hESDTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1310          }
1311          if(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother() > -1 &&
1312             ((AliAODMCParticle*)AODMCTrackArray->At(((AliAODMCParticle*)AODMCTrackArray->At(Photon->GetMother()))->GetMother()))->GetPdgCode() == 221){
1313             iPhotonMCInfo = 3;
1314          }
1315       }
1316    }
1317    
1318 }
1319 //________________________________________________________________________
1320 void AliAnalysisTaskGammaConvV1::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
1321 {
1322    
1323   Double_t magField = fInputEvent->GetMagneticField();
1324    if( magField  < 0.0 ){
1325       magField =  1.0;
1326    }
1327    else {
1328       magField =  -1.0;
1329    }
1330
1331    // Process True Photons
1332    TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(fMCStack);
1333    TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(fMCStack);
1334
1335    iPhotonMCInfo = 0;
1336    
1337    if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
1338    Int_t pdgCode[2] = {abs(posDaughter->GetPdgCode()),abs(negDaughter->GetPdgCode())};
1339    iPhotonMCInfo = 1;
1340    if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){
1341       FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1342       return;
1343    }
1344    else if(posDaughter->GetMother(0) == -1){
1345       FillPhotonCombinatorialBackgroundHist(TruePhotonCandidate, pdgCode);
1346       return;
1347    }
1348
1349    if(pdgCode[0]!=11 || pdgCode[1]!=11) return; //One Particle is not a electron
1350    
1351    if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
1352
1353    TParticle *Photon = TruePhotonCandidate->GetMCParticle(fMCStack);
1354    AliVTrack * electronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelNegative() );
1355    AliVTrack * positronCandidate = ((AliConversionCuts*)fCutArray->At(fiCut))->GetTrack(fInputEvent,TruePhotonCandidate->GetTrackLabelPositive() );
1356    Double_t deltaPhi = magField * TVector2::Phi_mpi_pi( electronCandidate->Phi()-positronCandidate->Phi());
1357    
1358    if(Photon->GetPdgCode() != 22){
1359       hESDTrueDalitzPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());
1360       return; // Mother is no Photon
1361    }
1362    
1363    if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
1364
1365    
1366    
1367    // True Photon
1368    if(fIsFromMBHeader){
1369       hESDTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1370    }
1371    hESDTrueGammaPsiPairDeltaPhi[fiCut]->Fill(deltaPhi,TruePhotonCandidate->GetPsiPair());  
1372    if(posDaughter->GetMother(0) <= fMCStack->GetNprimary()){
1373       // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
1374       if(fIsFromMBHeader){
1375          iPhotonMCInfo = 6;
1376          hESDTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());  
1377          hESDTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt()); // Allways Filled
1378
1379       }
1380       // (Not Filled for i6, Extra Signal Gamma (parambox) are secondary)
1381    }
1382    else{
1383       if(fIsFromMBHeader){
1384          iPhotonMCInfo = 2;
1385          hESDTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1386          if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1387             fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 3122){
1388             hESDTrueSecondaryConvGammaFromXFromLambdaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1389             iPhotonMCInfo = 5;
1390          }
1391          if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1392             fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
1393             hESDTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
1394             iPhotonMCInfo = 4;
1395          }
1396          if(fMCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
1397             fMCStack->Particle(fMCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 221){
1398             iPhotonMCInfo = 3;
1399          }
1400       }
1401    }
1402 }
1403 //________________________________________________________________________
1404 void AliAnalysisTaskGammaConvV1::ProcessAODMCParticles()
1405 {
1406
1407    TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1408
1409    // Loop over all primary MC particle
1410    for(Int_t i = 0; i < AODMCTrackArray->GetEntriesFast(); i++) {
1411
1412       AliAODMCParticle* particle = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(i));
1413       if (!particle) continue;
1414       if (!particle->IsPrimary()) continue;
1415
1416       Int_t isMCFromMBHeader = -1;
1417       if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1418          isMCFromMBHeader
1419             = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1420          if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1421       }
1422       
1423       if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1424       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kFALSE)){
1425          hMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1426          if(particle->GetMother() >-1){ // Meson Decay Gamma
1427             switch((static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetMother())))->GetPdgCode()){
1428             case 111: // Pi0
1429                hMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1430                break;
1431             case 113: // Rho0
1432                hMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1433                break;
1434             case 221: // Eta
1435                hMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1436                break;
1437             case 223: // Omega
1438                hMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1439                break;
1440             case 331: // Eta'
1441                hMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1442                break;
1443             case 333: // Phi
1444                hMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1445                break;
1446             case 3212: // Sigma
1447                hMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1448                break;
1449             }
1450          }
1451       }
1452       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(particle,AODMCTrackArray,kTRUE)){
1453          Double_t rConv = 0;
1454          for(Int_t daughterIndex=particle->GetDaughter(0);daughterIndex<=particle->GetDaughter(1);daughterIndex++){
1455             AliAODMCParticle *tmpDaughter = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(daughterIndex));
1456             if(!tmpDaughter) continue;
1457             if(abs(tmpDaughter->GetPdgCode()) == 11){
1458                rConv = sqrt( (tmpDaughter->Xv()*tmpDaughter->Xv()) + (tmpDaughter->Yv()*tmpDaughter->Yv()) );
1459             }
1460          }
1461          hMCConvGammaPt[fiCut]->Fill(particle->Pt());
1462          if (fDoPhotonQA > 0){
1463             hMCConvGammaR[fiCut]->Fill(rConv);
1464             hMCConvGammaEta[fiCut]->Fill(particle->Eta());
1465          }
1466       }
1467       // Converted MC Gamma
1468       if(fDoMesonAnalysis){
1469          if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1470             Double_t mesonY = 10.;
1471             if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1472                mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1473             } else{
1474                mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))
1475                   -((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1476             }
1477             Float_t weightedK0s= 1;
1478             if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1479                if (particle->Pt()>0.005){
1480                   weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1481                   //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1482                }
1483             }
1484             hMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1485             hMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1486             hMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1487          }         
1488          if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1489             ->MesonIsSelectedAODMC(particle,AODMCTrackArray,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
1490             AliAODMCParticle* daughter0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(0)));
1491             AliAODMCParticle* daughter1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(particle->GetDaughter(1)));
1492             Float_t weighted= 1;
1493             if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1494                if (particle->Pt()>0.005){
1495                   weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, 0x0, fInputEvent);
1496 //                   if(particle->GetPdgCode() == 221){
1497 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1498 //                   }
1499                }
1500             }
1501             Double_t mesonY = 10.;
1502             if(particle->E() - particle->Pz() == 0 || particle->E() + particle->Pz() == 0){
1503                mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1504             } else{
1505                mesonY = 0.5*(TMath::Log((particle->E()+particle->Pz()) / (particle->E()-particle->Pz())))
1506                   -((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1507             }
1508
1509             if(particle->GetPdgCode() == 111){
1510                hMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1511                hMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1512                if (fDoMesonQA > 0) hMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1513             } else if(particle->GetPdgCode() == 221){
1514                hMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1515                hMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1516                if (fDoMesonQA > 0) hMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1517             }
1518             
1519             // Check the acceptance for both gammas
1520             if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter0,AODMCTrackArray,kFALSE) &&
1521                ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedAODMC(daughter1,AODMCTrackArray,kFALSE)  &&
1522                ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1523                ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1524
1525                if(particle->GetPdgCode() == 111){
1526                   hMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
1527                } else if(particle->GetPdgCode() == 221){
1528                   hMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
1529                }
1530             }
1531          }
1532       }
1533    }
1534 }
1535 //________________________________________________________________________
1536 void AliAnalysisTaskGammaConvV1::ProcessMCParticles()
1537 {
1538    // Loop over all primary MC particle
1539    for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
1540       TParticle* particle = (TParticle *)fMCStack->Particle(i);
1541       if (!particle) continue;
1542
1543       Int_t isMCFromMBHeader = -1;
1544       if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
1545          isMCFromMBHeader
1546             = ((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent);
1547          if(isMCFromMBHeader == 0 && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 3) continue;
1548       }
1549
1550       if(!((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(particle->Phi(),fEventPlaneAngle,kFALSE)) continue;
1551       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
1552          hMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
1553          if(particle->GetMother(0) >-1){ // Meson Decay Gamma
1554             switch(fMCStack->Particle(particle->GetMother(0))->GetPdgCode()){
1555             case 111: // Pi0
1556                hMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
1557                break;
1558             case 113: // Rho0
1559                hMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
1560                break;
1561             case 221: // Eta
1562                hMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
1563                break;
1564             case 223: // Omega
1565                hMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
1566                break;
1567             case 331: // Eta'
1568                hMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
1569                break;
1570             case 333: // Phi
1571                hMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
1572                break;
1573             case 3212: // Sigma
1574                hMCDecayGammaSigmaPt[fiCut]->Fill(particle->Pt());
1575                break;
1576             }
1577          }
1578       }
1579       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
1580          hMCConvGammaPt[fiCut]->Fill(particle->Pt());
1581          if (fDoPhotonQA > 0){
1582             hMCConvGammaR[fiCut]->Fill(((TParticle*)fMCStack->Particle(particle->GetFirstDaughter()))->R());
1583             hMCConvGammaEta[fiCut]->Fill(particle->Eta());
1584          }
1585       } // Converted MC Gamma
1586       if(fDoMesonAnalysis){
1587          if(particle->GetPdgCode() == 310 && fDoMesonQA > 0){
1588             Double_t mesonY = 10.;
1589             if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1590                mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1591             } else{
1592                mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))
1593                   -((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1594             }
1595             Float_t weightedK0s= 1;
1596             if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1597                if (particle->Pt()>0.005){
1598                   weightedK0s= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1599                   //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1600                }
1601             }
1602             hMCK0sPt[fiCut]->Fill(particle->Pt(),weightedK0s);
1603             hMCK0sWOWeightPt[fiCut]->Fill(particle->Pt());
1604             hMCK0sPtY[fiCut]->Fill(particle->Pt(),mesonY,weightedK0s);
1605          }
1606          if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1607             ->MesonIsSelectedMC(particle,fMCStack,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
1608             TParticle* daughter0 = (TParticle*)fMCStack->Particle(particle->GetFirstDaughter());
1609             TParticle* daughter1 = (TParticle*)fMCStack->Particle(particle->GetLastDaughter());
1610
1611             Float_t weighted= 1;
1612             if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack, fInputEvent)){
1613                if (particle->Pt()>0.005){
1614                   weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),i, fMCStack, fInputEvent);
1615 //                   if(particle->GetPdgCode() == 221){
1616 //                      cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1617 //                   }
1618                }
1619             }
1620             Double_t mesonY = 10.;
1621             if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
1622                mesonY=10.-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1623             } else{
1624                mesonY = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())))
1625                   -((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift();
1626             }
1627
1628             if(particle->GetPdgCode() == 111){
1629                hMCPi0Pt[fiCut]->Fill(particle->Pt(),weighted); // All MC Pi0
1630                hMCPi0WOWeightPt[fiCut]->Fill(particle->Pt());
1631                if (fDoMesonQA > 0) hMCPi0PtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1632             } else if(particle->GetPdgCode() == 221){
1633                hMCEtaPt[fiCut]->Fill(particle->Pt(),weighted); // All MC Eta
1634                hMCEtaWOWeightPt[fiCut]->Fill(particle->Pt());
1635                if (fDoMesonQA > 0) hMCEtaPtY[fiCut]->Fill(particle->Pt(),mesonY,weighted); // All MC Pi0
1636             } 
1637
1638             // Check the acceptance for both gammas
1639             if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCStack,kFALSE) &&
1640                ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCStack,kFALSE)  &&
1641                ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter0->Phi(),fEventPlaneAngle,kFALSE) &&
1642                ((AliConversionCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(daughter1->Phi(),fEventPlaneAngle,kFALSE)){
1643
1644                if(particle->GetPdgCode() == 111){
1645                   hMCPi0InAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Pi0 with gamma in acc
1646                } else if(particle->GetPdgCode() == 221){
1647                   hMCEtaInAccPt[fiCut]->Fill(particle->Pt(),weighted); // MC Eta with gamma in acc
1648                }
1649             }
1650          }
1651       }
1652    }
1653 }
1654 //________________________________________________________________________
1655 void AliAnalysisTaskGammaConvV1::CalculatePi0Candidates(){
1656
1657    // Conversion Gammas
1658    if(fGammaCandidates->GetEntries()>1){
1659       for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries()-1;firstGammaIndex++){
1660          AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
1661          if (gamma0==NULL) continue;
1662          for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGammaCandidates->GetEntries();secondGammaIndex++){
1663             AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(secondGammaIndex));
1664             //Check for same Electron ID
1665             if (gamma1==NULL) continue;
1666             if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
1667                gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
1668                gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
1669                gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
1670
1671             AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
1672             pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
1673
1674             pi0cand->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
1675             if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
1676                 ->MesonIsSelected(pi0cand,kTRUE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
1677                hESDMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
1678                if(pi0cand->GetAlpha()<0.1)
1679                   hESDMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
1680                
1681                if (fDoMesonQA > 0){
1682                   if ( pi0cand->M() > 0.05 && pi0cand->M() < 0.17){
1683                      hESDMotherPi0PtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());  
1684                      hESDMotherPi0PtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());  
1685                      hESDMotherPi0PtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle()); 
1686                      
1687                   } 
1688                   if ( pi0cand->M() > 0.45 && pi0cand->M() < 0.65){
1689                      hESDMotherEtaPtY[fiCut]->Fill(pi0cand->Pt(),pi0cand->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());  
1690                      hESDMotherEtaPtAlpha[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetAlpha());  
1691                      hESDMotherEtaPtOpenAngle[fiCut]->Fill(pi0cand->Pt(),pi0cand->GetOpeningAngle());       
1692                   }
1693                }   
1694                if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGCalculation()){
1695                   Int_t zbin = 0;
1696                   Int_t mbin = 0;
1697
1698                   if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->BackgroundHandlerType() == 0){
1699                      zbin = fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
1700                      if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1701                         mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
1702                      } else {
1703                         mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
1704                      }
1705                   }
1706                   else{
1707                      zbin = fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
1708                      if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1709                         mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
1710                      } else {
1711                         mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
1712                      }
1713                   }
1714                   Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),(Double_t)zbin,(Double_t)mbin};
1715                   sESDMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
1716                }
1717                
1718
1719                if(fIsMC){
1720                   if(fInputEvent->IsA()==AliESDEvent::Class())
1721                      ProcessTrueMesonCandidates(pi0cand,gamma0,gamma1);
1722                   if(fInputEvent->IsA()==AliAODEvent::Class())
1723                      ProcessTrueMesonCandidatesAOD(pi0cand,gamma0,gamma1);
1724                }
1725                if (fDoMesonQA == 2){
1726                   fInvMass = pi0cand->M();
1727                   fPt  = pi0cand->Pt();
1728                   if (abs(gamma0->GetDCAzToPrimVtx()) < abs(gamma1->GetDCAzToPrimVtx())){
1729                      fDCAzGammaMin = gamma0->GetDCAzToPrimVtx();
1730                      fDCAzGammaMax = gamma1->GetDCAzToPrimVtx();
1731                   } else {
1732                      fDCAzGammaMin = gamma1->GetDCAzToPrimVtx();
1733                      fDCAzGammaMax = gamma0->GetDCAzToPrimVtx();
1734                   }
1735                   iFlag = pi0cand->GetMesonQuality();
1736 //                   cout << "gamma 0: " << gamma0->GetV0Index()<< "\t" << gamma0->GetPx() << "\t" << gamma0->GetPy() << "\t" <<  gamma0->GetPz() << "\t" << endl; 
1737 //                   cout << "gamma 1: " << gamma1->GetV0Index()<< "\t"<< gamma1->GetPx() << "\t" << gamma1->GetPy() << "\t" <<  gamma1->GetPz() << "\t" << endl; 
1738 //                    cout << "pi0: "<<fInvMass << "\t" << fPt <<"\t" << fDCAzGammaMin << "\t" << fDCAzGammaMax << "\t" << (Int_t)iFlag << "\t" << (Int_t)iMesonMCInfo <<endl;
1739                   if (fIsHeavyIon == 1 && fPt > 0.399 && fPt < 20. ) {
1740                      if (fInvMass > 0.08 && fInvMass < 0.2) tESDMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
1741                      if ((fInvMass > 0.45 && fInvMass < 0.6) &&  (fPt > 0.999 && fPt < 20.) )tESDMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
1742                   } else if (fPt > 0.299 && fPt < 20. )  {
1743                      if ( (fInvMass > 0.08 && fInvMass < 0.2) || (fInvMass > 0.45 && fInvMass < 0.6)) tESDMesonsInvMassPtDcazMinDcazMaxFlag[fiCut]->Fill();
1744                   }   
1745                }
1746             }
1747             delete pi0cand;
1748             pi0cand=0x0;
1749          }
1750       }
1751    }
1752 }
1753 //______________________________________________________________________
1754 void AliAnalysisTaskGammaConvV1::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1755 {
1756    // Process True Mesons
1757    AliStack *MCStack = fMCEvent->Stack();
1758    iMesonMCInfo = 0;
1759    if(TrueGammaCandidate0->GetV0Index()<fInputEvent->GetNumberOfV0s()){
1760       Bool_t isTruePi0 = kFALSE;
1761       Bool_t isTrueEta = kFALSE;
1762       Bool_t isTruePi0Dalitz = kFALSE;
1763       Bool_t isTrueEtaDalitz = kFALSE;
1764       Bool_t gamma0DalitzCand = kFALSE;
1765       Bool_t gamma1DalitzCand = kFALSE;
1766       Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
1767       Int_t gamma0MotherLabel = -1;
1768       if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1769          // Daughters Gamma 0
1770          TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
1771          TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
1772          TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
1773          if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1774             if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1775                if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1776                   gamma0MotherLabel=gammaMC0->GetFirstMother();
1777                }
1778             }
1779             if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1780                gamma0DalitzCand = kTRUE;
1781                gamma0MotherLabel=-111;
1782             }
1783             if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
1784                gamma0DalitzCand = kTRUE;
1785                gamma0MotherLabel=-221;
1786             }
1787          }
1788       }
1789       if(TrueGammaCandidate1->GetV0Index()<fInputEvent->GetNumberOfV0s()){
1790          Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
1791          Int_t gamma1MotherLabel = -1;
1792          if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1793             // Daughters Gamma 1
1794             TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
1795             TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
1796             TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
1797             if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1798                if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
1799                   if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1800                      gamma1MotherLabel=gammaMC1->GetFirstMother();
1801                   }
1802                }
1803                if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1804                  gamma1DalitzCand = kTRUE;
1805                  gamma1MotherLabel=-111;
1806                }
1807                if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
1808                  gamma1DalitzCand = kTRUE;
1809                  gamma1MotherLabel=-221;
1810                }
1811             }
1812          }
1813          if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1814             if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
1815                isTruePi0=kTRUE;
1816             }
1817             if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 221){
1818                isTrueEta=kTRUE;
1819             }
1820          }
1821          
1822          //Identify Dalitz candidate
1823          if (gamma1DalitzCand || gamma0DalitzCand){
1824             if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
1825                if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1826                if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;   
1827             }   
1828             if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
1829                if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
1830                if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;   
1831             }
1832          }
1833          
1834          
1835          if(isTruePi0 || isTrueEta){// True Pion or Eta
1836             hESDTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt()); 
1837             if (fDoMesonQA > 0){
1838                if (isTruePi0){
1839                   if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
1840                      hESDTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); 
1841                      hESDTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); 
1842                      hESDTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); 
1843                   }
1844                } else if (isTrueEta){   
1845                   if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
1846                      hESDTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); 
1847                      hESDTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); 
1848                      hESDTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); 
1849                   }
1850                }
1851             }   
1852             if(gamma0MotherLabel >= MCStack->GetNprimary()){ // Secondary Meson
1853                Int_t secMotherLabel = ((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0);
1854                Float_t weightedSec= 1;
1855                if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, fMCStack, fInputEvent) && MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
1856                   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
1857                   //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
1858                }
1859                hESDTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
1860                iMesonMCInfo = 2;
1861                if (secMotherLabel >-1){
1862                   if(MCStack->Particle(secMotherLabel)->GetPdgCode()==310){
1863                      iMesonMCInfo = 4;
1864                      hESDTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
1865                      if (fDoMesonQA > 0)hESDTrueK0sWithPi0DaughterMCPt[fiCut]
1866                                        ->Fill(MCStack->Particle(secMotherLabel)->Pt());
1867                   }
1868                   if(MCStack->Particle(secMotherLabel)->GetPdgCode()==221){
1869                      iMesonMCInfo = 3;
1870                      hESDTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
1871                      if (fDoMesonQA > 0)hESDTrueEtaWithPi0DaughterMCPt[fiCut]
1872                                        ->Fill(MCStack->Particle(secMotherLabel)->Pt());
1873                   }
1874                }
1875             }else{ // Only primary pi0 for efficiency calculation
1876                iMesonMCInfo = 6;
1877                Float_t weighted= 1;
1878                if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, fMCStack, fInputEvent)){
1879                   if (((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt()>0.005){
1880                      weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, fMCStack, fInputEvent);
1881 //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
1882                   }
1883                }
1884                hESDTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
1885                hESDTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1886                pESDTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
1887               
1888                
1889                if (fDoMesonQA > 0){
1890                   if(isTruePi0){ // Only primary pi0 for resolution
1891                      hESDTruePrimaryPi0MCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
1892                   }
1893                   if (isTrueEta){ // Only primary eta for resolution
1894                      hESDTruePrimaryEtaMCPtResolPt[fiCut]->Fill(((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),(Pi0Candidate->Pt()-((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt())/((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt(),weighted);
1895                   }
1896                }
1897             }
1898          } else if(!isTruePi0 && !isTrueEta){ // Background
1899             if (fDoMesonQA > 0){
1900                if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
1901                   hESDTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1902                   iMesonMCInfo = 1;
1903                } else { // No photon or without mother
1904                   hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1905                }
1906             }
1907             if( isTruePi0Dalitz || isTrueEtaDalitz ){
1908                // Dalitz
1909                iMesonMCInfo = 5;
1910                hESDTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1911             } else if (gamma0DalitzCand || gamma1DalitzCand){
1912                if (fDoMesonQA > 0)hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
1913             }   
1914          }
1915       }
1916    }
1917 }
1918 //______________________________________________________________________
1919 void AliAnalysisTaskGammaConvV1::ProcessTrueMesonCandidatesAOD(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
1920 {
1921
1922    // Process True Mesons
1923    TClonesArray *AODMCTrackArray = dynamic_cast<TClonesArray*>(fInputEvent->FindListObject(AliAODMCParticle::StdBranchName()));
1924    Bool_t isTruePi0 = kFALSE;
1925    Bool_t isTrueEta = kFALSE;
1926    Bool_t isTruePi0Dalitz = kFALSE;
1927    Bool_t isTrueEtaDalitz = kFALSE;
1928    Bool_t gamma0DalitzCand = kFALSE;
1929    Bool_t gamma1DalitzCand = kFALSE;
1930    
1931    AliAODMCParticle *positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelPositive()));
1932    AliAODMCParticle *negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate0->GetMCLabelNegative()));
1933
1934    iMesonMCInfo = 0;
1935    Int_t gamma0MCLabel = -1;
1936    Int_t gamma0MotherLabel = -1;
1937    if(!positiveMC||!negativeMC)
1938       return;
1939    
1940    if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1941       gamma0MCLabel = positiveMC->GetMother();
1942    }
1943
1944    if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1945       // Daughters Gamma 0
1946       AliAODMCParticle * gammaMC0 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MCLabel));
1947       if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1948          if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...     
1949             if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
1950                gamma0MotherLabel=gammaMC0->GetMother();
1951             }
1952          }
1953          if(gammaMC0->GetPdgCode() ==111){ // Dalitz candidate
1954             gamma0DalitzCand = kTRUE;
1955             gamma0MotherLabel=-111;
1956          }
1957          if(gammaMC0->GetPdgCode() ==221){ // Dalitz candidate
1958             gamma0DalitzCand = kTRUE;
1959             gamma0MotherLabel=-221;
1960          }
1961       }
1962    }
1963    positiveMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelPositive()));
1964    negativeMC = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(TrueGammaCandidate1->GetMCLabelNegative()));
1965    
1966    Int_t gamma1MCLabel = -1;
1967    Int_t gamma1MotherLabel = -1;
1968    if(!positiveMC||!negativeMC)
1969       return;
1970    
1971    if(positiveMC->GetMother()>-1&&(negativeMC->GetMother() == positiveMC->GetMother())){
1972       gamma1MCLabel = positiveMC->GetMother();
1973    }
1974    if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
1975       // Daughters Gamma 1
1976       AliAODMCParticle * gammaMC1 = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MCLabel));
1977       if(abs(negativeMC->GetPdgCode())==11 && abs(positiveMC->GetPdgCode())==11){  // Electrons ...
1978          if(((positiveMC->GetMCProcessCode())) == 5 && ((negativeMC->GetMCProcessCode())) == 5){ // ... From Conversion ...     
1979             if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
1980                gamma1MotherLabel=gammaMC1->GetMother();
1981             }
1982          }
1983          if(gammaMC1->GetPdgCode() ==111 ){ // Dalitz candidate
1984                  gamma1DalitzCand = kTRUE;
1985                  gamma1MotherLabel=-111;
1986          }
1987          if(gammaMC1->GetPdgCode() ==221){ // Dalitz candidate
1988             gamma1DalitzCand = kTRUE;
1989             gamma1MotherLabel=-221;
1990          }
1991       }
1992    }
1993    if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
1994       if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 111){
1995          isTruePi0=kTRUE;
1996       }
1997       if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetPdgCode() == 221){
1998          isTrueEta=kTRUE;
1999       }
2000    }
2001    
2002    //Identify Dalitz candidate
2003    if (gamma1DalitzCand || gamma0DalitzCand){
2004       if (gamma0DalitzCand && gamma0MCLabel >=0 && gamma0MCLabel==gamma1MotherLabel){
2005          if (gamma0MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2006          if (gamma0MotherLabel == -221) isTrueEtaDalitz = kTRUE;   
2007       }   
2008       if (gamma1DalitzCand && gamma1MCLabel >=0 && gamma1MCLabel==gamma0MotherLabel){
2009          if (gamma1MotherLabel == -111) isTruePi0Dalitz = kTRUE;
2010          if (gamma1MotherLabel == -221) isTrueEtaDalitz = kTRUE;   
2011       }
2012    }
2013             
2014    if(isTruePi0 || isTrueEta){// True Pion or Eta
2015       hESDTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2016       if (fDoMesonQA > 0){
2017          if (isTruePi0){
2018             if ( Pi0Candidate->M() > 0.05 && Pi0Candidate->M() < 0.17){
2019                hESDTruePi0PtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift());
2020                hESDTruePi0PtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); 
2021                hESDTruePi0PtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); 
2022             }
2023          } else if (isTrueEta){   
2024             if ( Pi0Candidate->M() > 0.45 && Pi0Candidate->M() < 0.65){
2025                hESDTrueEtaPtY[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->Rapidity()-((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()); 
2026                hESDTrueEtaPtAlpha[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetAlpha()); 
2027                hESDTrueEtaPtOpenAngle[fiCut]->Fill(Pi0Candidate->Pt(),Pi0Candidate->GetOpeningAngle()); 
2028             }
2029          }
2030       }
2031       if(!(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma0MotherLabel))->IsPrimary())){ // Secondary Meson
2032          Int_t secMotherLabel = static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->GetMother();
2033          Float_t weightedSec= 1;
2034          if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(secMotherLabel, 0x0, fInputEvent) && static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2035             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
2036             //cout << "MC input \t"<<i << "\t" <<  particle->Pt()<<"\t"<<weighted << endl;
2037          }
2038          hESDTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2039          iMesonMCInfo = 2;   
2040          if (secMotherLabel >-1){
2041             if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==310){
2042                iMesonMCInfo = 4;
2043                hESDTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2044                if (fDoMesonQA > 0)hESDTrueK0sWithPi0DaughterMCPt[fiCut]
2045                                  ->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2046             }
2047             if(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->GetPdgCode()==221){
2048                iMesonMCInfo = 3;
2049                hESDTrueSecondaryMotherFromEtaInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weightedSec);
2050                if (fDoMesonQA > 0)hESDTrueEtaWithPi0DaughterMCPt[fiCut]
2051                                  ->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(secMotherLabel))->Pt());
2052             }
2053          }
2054       }else{ // Only primary pi0 for efficiency calculation
2055          Float_t weighted= 1;
2056          iMesonMCInfo = 6;
2057          if(((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(gamma1MotherLabel, 0x0, fInputEvent)){
2058             if (static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt()>0.005){
2059                weighted= ((AliConversionCuts*)fCutArray->At(fiCut))->GetWeightForMeson(fV0Reader->GetPeriodName(),gamma1MotherLabel, 0x0, fInputEvent);
2060                //                      cout << "rec \t " <<gamma1MotherLabel << "\t" <<  weighted << endl;
2061             }
2062          }
2063          hESDTruePrimaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2064          hESDTruePrimaryMotherW0WeightingInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2065          pESDTruePrimaryMotherWeightsInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt(),weighted);
2066               
2067          if (fDoMesonQA > 0){
2068             if(isTruePi0){ // Only primary pi0 for resolution
2069                hESDTruePrimaryPi0MCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2070                                                           (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2071                
2072             }
2073             if (isTrueEta){ // Only primary eta for resolution
2074                hESDTruePrimaryEtaMCPtResolPt[fiCut]->Fill(static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),
2075                                                           (Pi0Candidate->Pt()-static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt())/static_cast<AliAODMCParticle*>(AODMCTrackArray->At(gamma1MotherLabel))->Pt(),weighted);
2076             }
2077          }
2078       }
2079    } else if(!isTruePi0 && !isTrueEta) { // Background
2080       if (fDoMesonQA > 0){
2081          if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
2082             hESDTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2083             iMesonMCInfo = 1;
2084          } else { // No photon or without mother
2085             hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2086          }
2087       }
2088       if( isTruePi0Dalitz || isTrueEtaDalitz ){
2089          // Dalitz
2090          iMesonMCInfo = 5;
2091          hESDTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2092       } else if (gamma0DalitzCand || gamma1DalitzCand){
2093          if (fDoMesonQA > 0)hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
2094       }
2095    }
2096 }
2097 //________________________________________________________________________
2098 void AliAnalysisTaskGammaConvV1::CalculateBackground(){
2099
2100    Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2101    Int_t mbin = 0;
2102
2103    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2104       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2105    } else {
2106       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2107    }
2108
2109    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
2110
2111       for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2112          AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2113          for(Int_t iCurrent2=iCurrent+1;iCurrent2<fGammaCandidates->GetEntries();iCurrent2++){
2114             for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){
2115                AliAODConversionPhoton currentEventGoodV02 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent2));
2116
2117                if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGProbability()){
2118                   AliAODConversionMother *backgroundCandidateProb = new AliAODConversionMother(&currentEventGoodV0,&currentEventGoodV02);
2119                   Double_t massBGprob = backgroundCandidateProb->M();
2120                   if(massBGprob>0.1 && massBGprob<0.14){
2121                      if(fRandom.Rndm()>fBGHandler[fiCut]->GetBGProb(zbin,mbin)){
2122                         delete backgroundCandidateProb;
2123                         continue;
2124                      }
2125                   }
2126                   delete backgroundCandidateProb;
2127                   backgroundCandidateProb = 0x0;
2128                }
2129
2130                RotateParticle(&currentEventGoodV02);
2131                AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&currentEventGoodV02);
2132                backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2133                if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2134                    ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2135                   hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2136                   Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2137                   sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2138                }
2139                delete backgroundCandidate;
2140                backgroundCandidate = 0x0;
2141             }
2142          }
2143       }
2144    }else{
2145       AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
2146
2147       if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2148          for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2149             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2150             if(fMoveParticleAccordingToVertex == kTRUE){
2151                bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2152             }
2153
2154             for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2155                AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2156                for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2157                   AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2158                   if(fMoveParticleAccordingToVertex == kTRUE){
2159                      MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2160                   }
2161                   if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2162                      RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2163                   }
2164
2165                   AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2166                   backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2167                   if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2168                       ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2169                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2170                      Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2171                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2172                   }
2173                   delete backgroundCandidate;
2174                   backgroundCandidate = 0x0;
2175                }
2176             }
2177          }
2178       }
2179       else{
2180          for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
2181             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
2182             if(previousEventV0s){
2183                if(fMoveParticleAccordingToVertex == kTRUE){
2184                   bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
2185                }
2186                for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2187                   AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2188                   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2189
2190                      AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
2191
2192                      if(fMoveParticleAccordingToVertex == kTRUE){
2193                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
2194                      }
2195                      if(((AliConversionCuts*)fCutArray->At(fiCut))->GetInPlaneOutOfPlaneCut() != 0){
2196                         RotateParticleAccordingToEP(&previousGoodV0,bgEventVertex->fEP,fEventPlaneAngle);
2197                      }
2198
2199
2200                      AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
2201                      backgroundCandidate->CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2202                      if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2203                          ->MesonIsSelected(backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift()))){
2204                         hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
2205                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),(Double_t)zbin,(Double_t)mbin};
2206                         sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
2207                      }
2208                      delete backgroundCandidate;
2209                      backgroundCandidate = 0x0;
2210                   }
2211                }
2212             }
2213          }
2214       }
2215    }
2216 }
2217 //________________________________________________________________________
2218 void AliAnalysisTaskGammaConvV1::CalculateBackgroundRP(){
2219
2220    Int_t zbin= fBGHandlerRP[fiCut]->GetZBinIndex(fInputEvent->GetPrimaryVertex()->GetZ());
2221    Int_t mbin = 0;
2222    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2223       mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fV0Reader->GetNumberOfPrimaryTracks());
2224    } else {
2225       mbin = fBGHandlerRP[fiCut]->GetMultiplicityBinIndex(fGammaCandidates->GetEntries());
2226    }
2227
2228
2229    //Rotation Method
2230    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
2231       // Correct for the number of rotations
2232       // BG is for rotation the same, except for factor NRotations
2233       Double_t weight=1./Double_t(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents());
2234
2235       for(Int_t firstGammaIndex=0;firstGammaIndex<fGammaCandidates->GetEntries();firstGammaIndex++){
2236
2237          AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(firstGammaIndex));
2238          if (gamma0==NULL) continue;
2239          for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGammaCandidates->GetEntries();secondGammaIndex++){
2240             AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(secondGammaIndex));
2241             if (gamma1 == NULL) continue;
2242             if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(gamma1,fInputEvent))continue;
2243             for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->GetNumberOfBGEvents();nRandom++){
2244
2245                RotateParticle(gamma1);
2246
2247                AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2248                backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2249                if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2250                   ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
2251                   hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2252                   Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2253                   sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2254                }
2255             }
2256          }
2257       }
2258    }
2259    else{
2260       // Do Event Mixing
2261       for(Int_t nEventsInBG=0;nEventsInBG <fBGHandlerRP[fiCut]->GetNBGEvents(fGammaCandidates,fInputEvent);nEventsInBG++){
2262
2263          AliGammaConversionPhotonVector *previousEventGammas = fBGHandlerRP[fiCut]->GetBGGoodGammas(fGammaCandidates,fInputEvent,nEventsInBG);
2264
2265          if(previousEventGammas){
2266             // test weighted background
2267             Double_t weight=1.0;
2268             // Correct for the number of eventmixing:
2269             // N gammas -> (N-1) + (N-2) +(N-3) ...+ (N-(N-1))  using sum formula sum(i)=N*(N-1)/2  -> N*(N-1)/2
2270             // real combinations (since you cannot combine a photon with its own)
2271             // but BG leads to N_{a}*N_{b} combinations
2272             weight*=0.5*(Double_t(fGammaCandidates->GetEntries()-1))/Double_t(previousEventGammas->size());
2273
2274             for(Int_t iCurrent=0;iCurrent<fGammaCandidates->GetEntries();iCurrent++){
2275
2276                AliAODConversionPhoton *gamma0 = (AliAODConversionPhoton*)(fGammaCandidates->At(iCurrent));
2277
2278                for(UInt_t iPrevious=0;iPrevious<previousEventGammas->size();iPrevious++){
2279
2280                   AliAODConversionPhoton *gamma1 = (AliAODConversionPhoton*)(previousEventGammas->at(iPrevious));
2281
2282                   AliAODConversionMother backgroundCandidate(gamma0,gamma1);
2283                   backgroundCandidate.CalculateDistanceOfClossetApproachToPrimVtx(fInputEvent->GetPrimaryVertex());
2284                   if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))
2285                      ->MesonIsSelected(&backgroundCandidate,kFALSE,((AliConversionCuts*)fCutArray->At(fiCut))->GetEtaShift())){
2286                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate.M(),backgroundCandidate.Pt());
2287                      Double_t sparesFill[4] = {backgroundCandidate.M(),backgroundCandidate.Pt(),(Double_t)zbin,(Double_t)mbin};
2288                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,weight);
2289                   }
2290                }
2291             }
2292          }
2293       }
2294    }
2295 }
2296 //________________________________________________________________________
2297 void AliAnalysisTaskGammaConvV1::RotateParticle(AliAODConversionPhoton *gamma){
2298    Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation();
2299    Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
2300    Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
2301    gamma->RotateZ(rotationValue);
2302 }
2303
2304 //________________________________________________________________________
2305 void AliAnalysisTaskGammaConvV1::RotateParticleAccordingToEP(AliAODConversionPhoton *gamma, Double_t previousEventEP, Double_t thisEventEP){
2306
2307    previousEventEP=previousEventEP+TMath::Pi();
2308    thisEventEP=thisEventEP+TMath::Pi();
2309    Double_t rotationValue= thisEventEP-previousEventEP;
2310    gamma->RotateZ(rotationValue);
2311 }
2312
2313 //________________________________________________________________________
2314 void AliAnalysisTaskGammaConvV1::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
2315    //see header file for documentation
2316
2317    Double_t dx = vertex->fX - fInputEvent->GetPrimaryVertex()->GetX();
2318    Double_t dy = vertex->fY - fInputEvent->GetPrimaryVertex()->GetY();
2319    Double_t dz = vertex->fZ - fInputEvent->GetPrimaryVertex()->GetZ();
2320
2321    Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
2322    particle->SetConversionPoint(movedPlace);
2323 }
2324 //________________________________________________________________________
2325 void AliAnalysisTaskGammaConvV1::UpdateEventByEventData(){
2326    //see header file for documentation
2327    if(fGammaCandidates->GetEntries() >0 ){
2328       if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
2329          fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
2330       }
2331       else{ // means we use #V0s for multiplicity
2332          fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
2333       }
2334    }
2335 }
2336
2337
2338 //________________________________________________________________________
2339 void AliAnalysisTaskGammaConvV1::FillPhotonCombinatorialBackgroundHist(AliAODConversionPhoton *TruePhotonCandidate, Int_t pdgCode[])
2340 {
2341    // 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
2342    if(pdgCode[0]==11   && pdgCode[1]==11){if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),0);}
2343    else if( (pdgCode[0]==11   && pdgCode[1]==211) || (pdgCode[0]==211  && pdgCode[1]==11) )
2344       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),1);}
2345    else if( (pdgCode[0]==11   && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==11) )
2346       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),2);}
2347    else if( (pdgCode[0]==11   && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==11) )
2348       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),3);}
2349    else if( (pdgCode[0]==11   && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==11) )
2350       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),4);}
2351    else if(  pdgCode[0]==211  && pdgCode[1]==211 ){if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),5);}
2352    else if( (pdgCode[0]==211  && pdgCode[1]==321) || (pdgCode[0]==321  && pdgCode[1]==211) )
2353       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),6);}
2354    else if( (pdgCode[0]==211  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==211) )
2355       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),7);}
2356    else if( (pdgCode[0]==211  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==211) )
2357       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),8);}
2358    else if(  pdgCode[0]==321  && pdgCode[1]==321 ){if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),9);}
2359    else if( (pdgCode[0]==321  && pdgCode[1]==2212) || (pdgCode[0]==2212 && pdgCode[1]==321) )
2360       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),10);}
2361    else if( (pdgCode[0]==321  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==321) )
2362       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),11);}
2363    else if(  pdgCode[0]==2212   && pdgCode[1]==2212  ){if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),12);}
2364    else if( (pdgCode[0]==2212  && pdgCode[1]==13) || (pdgCode[0]==13   && pdgCode[1]==2212) )
2365       {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),13);}
2366    else if(  pdgCode[0]==13   && pdgCode[1]==13  ){if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),14);}
2367    else {if(fIsFromMBHeader)hESDCombinatorialPt[fiCut]->Fill(TruePhotonCandidate->Pt(),15);}
2368 }
2369 //________________________________________________________________________
2370 void AliAnalysisTaskGammaConvV1::RelabelAODPhotonCandidates(Bool_t mode){
2371
2372    // Relabeling For AOD Event
2373    // ESDiD -> AODiD
2374    // MCLabel -> AODMCLabel
2375    
2376    if(mode){
2377       fMCStackPos = new Int_t[fReaderGammas->GetEntries()];
2378       fMCStackNeg = new Int_t[fReaderGammas->GetEntries()];
2379       fESDArrayPos = new Int_t[fReaderGammas->GetEntries()];
2380       fESDArrayNeg = new Int_t[fReaderGammas->GetEntries()];
2381    }
2382    
2383    for(Int_t iGamma = 0;iGamma<fReaderGammas->GetEntries();iGamma++){
2384       AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(iGamma);
2385       if(!PhotonCandidate) continue;
2386       if(!mode){// Back to ESD Labels
2387          PhotonCandidate->SetMCLabelPositive(fMCStackPos[iGamma]);
2388          PhotonCandidate->SetMCLabelNegative(fMCStackNeg[iGamma]);
2389          PhotonCandidate->SetLabelPositive(fESDArrayPos[iGamma]);
2390          PhotonCandidate->SetLabelNegative(fESDArrayNeg[iGamma]);
2391          continue;
2392       }
2393       fMCStackPos[iGamma] =  PhotonCandidate->GetMCLabelPositive();
2394       fMCStackNeg[iGamma] =  PhotonCandidate->GetMCLabelNegative();
2395       fESDArrayPos[iGamma] = PhotonCandidate->GetTrackLabelPositive();
2396       fESDArrayNeg[iGamma] = PhotonCandidate->GetTrackLabelNegative();
2397
2398       Bool_t AODLabelPos = kFALSE;
2399       Bool_t AODLabelNeg = kFALSE;
2400
2401       for(Int_t i = 0; i<fInputEvent->GetNumberOfTracks();i++){
2402          AliAODTrack *tempDaughter = static_cast<AliAODTrack*>(fInputEvent->GetTrack(i));
2403          if(!AODLabelPos){
2404             if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelPositive() ){
2405                PhotonCandidate->SetMCLabelPositive(abs(tempDaughter->GetLabel()));
2406                PhotonCandidate->SetLabelPositive(i);
2407                AODLabelPos = kTRUE;
2408             }
2409          }
2410          if(!AODLabelNeg){
2411             if( tempDaughter->GetID() == PhotonCandidate->GetTrackLabelNegative()){
2412                PhotonCandidate->SetMCLabelNegative(abs(tempDaughter->GetLabel()));
2413                PhotonCandidate->SetLabelNegative(i);
2414                AODLabelNeg = kTRUE;
2415             }
2416          }
2417          if(AODLabelNeg && AODLabelPos){
2418             break;
2419          }
2420       }
2421       if(!AODLabelPos || !AODLabelNeg){
2422          cout<<"WARNING!!! AOD TRACKS NOT FOUND FOR"<<endl;
2423       }
2424    }
2425    
2426    
2427    if(!mode){
2428       delete[] fMCStackPos;
2429       delete[] fMCStackNeg;
2430       delete[] fESDArrayPos;
2431       delete[] fESDArrayNeg;
2432    }
2433 }
2434
2435 void AliAnalysisTaskGammaConvV1::SetLogBinningXTH2(TH2* histoRebin){
2436    TAxis *axisafter = histoRebin->GetXaxis(); 
2437    Int_t bins = axisafter->GetNbins();
2438    Double_t from = axisafter->GetXmin();
2439    Double_t to = axisafter->GetXmax();
2440    Double_t *newbins = new Double_t[bins+1];
2441    newbins[0] = from;
2442    Double_t factor = TMath::Power(to/from, 1./bins);
2443    for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1];
2444    axisafter->Set(bins, newbins);
2445    delete [] newbins;
2446
2447 }
2448
2449 //________________________________________________________________________
2450 void AliAnalysisTaskGammaConvV1::Terminate(const Option_t *)
2451 {
2452
2453    //fOutputContainer->Print(); // Will crash on GRID
2454 }