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