]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGGA/GammaConv/AliAnalysisTaskGammaConvV1.cxx
added material task
[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                                    *
5  * Version 1.0                                                            *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class used to do analysis on conversion pairs
19 //---------------------------------------------
20 ///////////////////////////////////////////////
21 #include "TChain.h"
22 #include "TTree.h"
23 #include "TH1F.h"
24 #include "TH2F.h"
25 #include "TH3F.h"
26 #include "THnSparse.h"
27 #include "TCanvas.h"
28 #include "TNtuple.h"
29 #include "AliAnalysisTask.h"
30 #include "AliAnalysisManager.h"
31 #include "AliESDEvent.h"
32 #include "AliESDInputHandler.h"
33 #include "AliMCEventHandler.h"
34 #include "AliMCEvent.h"
35 #include "AliMCParticle.h"
36 #include "AliCentrality.h"
37 #include "AliESDVZERO.h"
38 #include "AliESDpid.h"
39 #include "AliAnalysisTaskGammaConvV1.h"
40 #include "AliVParticle.h"
41 #include "AliESDtrackCuts.h"
42 #include "AliKFVertex.h"
43 #include "AliV0ReaderV1.h"
44 #include "AliGammaConversionAODBGHandler.h"
45 #include "AliGenCocktailEventHeader.h"
46
47 ClassImp(AliAnalysisTaskGammaConvV1)
48
49 //________________________________________________________________________
50 AliAnalysisTaskGammaConvV1::AliAnalysisTaskGammaConvV1(): AliAnalysisTaskSE(),
51    fV0Reader(NULL),
52    fBGHandler(NULL),
53    fESDEvent(NULL),
54    fMCEvent(NULL),
55    fMCStack(NULL),
56    fCutFolder(NULL),
57    fESDList(NULL),
58    fBackList(NULL),
59    fMotherList(NULL),
60    fTrueList(NULL),
61    fMCList(NULL),
62    fHeaderNameList(NULL),
63    fOutputContainer(0),
64    fReaderGammas(NULL),
65    fGoodGammas(NULL),
66    fCutArray(NULL),
67    fConversionCuts(NULL),
68    fMesonCutArray(NULL),
69    fMesonCuts(NULL),
70    hESDConvGammaPt(NULL),
71    hESDMotherInvMassPt(NULL),
72    sESDMotherInvMassPtZM(NULL),
73    hESDMotherBackInvMassPt(NULL),
74    sESDMotherBackInvMassPtZM(NULL),
75    hESDMotherInvMassEalpha(NULL),
76    hMCAllGammaPt(NULL),
77    hMCDecayGammaPi0Pt(NULL),
78    hMCDecayGammaRhoPt(NULL),
79    hMCDecayGammaEtaPt(NULL),
80    hMCDecayGammaOmegaPt(NULL),
81    hMCDecayGammaEtapPt(NULL),
82    hMCDecayGammaPhiPt(NULL),
83    hMCConvGammaPt(NULL),
84    hMCConvGammaR(NULL),
85    hMCConvGammaEta(NULL),
86    hMCPi0Pt(NULL),
87    hMCEtaPt(NULL),
88    hMCPi0InAccPt(NULL),
89    hMCEtaInAccPt(NULL),
90    hESDTrueMotherInvMassPt(NULL),
91    hESDTruePrimaryMotherInvMassMCPt(NULL),
92    hESDTruePrimaryPi0ESDPtMCPt(NULL),
93    hESDTrueSecondaryMotherInvMassPt(NULL),
94    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
95    hESDTrueBckGGInvMassPt(NULL),
96    hESDTrueBckContInvMassPt(NULL),
97    hESDTrueMotherDalitzInvMassPt(NULL),
98    hESDTrueConvGammaPt(NULL),
99    hESDTrueTwoElecCombPt(NULL),
100    hESDTrueTwoPionCombPt(NULL),
101    hESDTrueElecPionCombPt(NULL),
102    hESDTrueCombPt(NULL),
103    hESDTruePrimaryConvGammaPt(NULL),
104    hESDTruePrimaryConvGammaR(NULL),
105    hESDTruePrimaryConvGammaEta(NULL),
106    hESDTruePrimaryConvGammaESDPtMCPt(NULL),
107    hESDTrueSecondaryConvGammaPt(NULL),
108    hESDTrueSecondaryConvGammaFromK0sPt(NULL),
109    hESDTrueSecondaryConvGammaFromXFromK0sPt(NULL),
110    hNEvents(NULL),
111    hNGoodESDTracks(NULL),
112    hNV0Tracks(NULL),
113    fRandom(0),
114    fUnsmearedPx(NULL),
115    fUnsmearedPy(NULL),
116    fUnsmearedPz(NULL),
117    fUnsmearedE(NULL),
118    fIsHeavyIon(kFALSE),
119    fDoMesonAnalysis(kTRUE)
120
121 {
122    // default Constructor
123    DefineInput(0, TChain::Class());
124    DefineOutput(1, TList::Class());
125 }
126
127 //________________________________________________________________________
128 AliAnalysisTaskGammaConvV1::AliAnalysisTaskGammaConvV1(const char *name):
129    AliAnalysisTaskSE(name),
130    fV0Reader(NULL),
131    fBGHandler(NULL),
132    fESDEvent(NULL),
133    fMCEvent(NULL),
134    fMCStack(NULL),
135    fCutFolder(NULL),
136    fESDList(NULL),
137    fBackList(NULL),
138    fMotherList(NULL),
139    fTrueList(NULL),
140    fMCList(NULL),
141    fHeaderNameList(NULL),
142    fOutputContainer(0),
143    fReaderGammas(NULL),
144    fGoodGammas(NULL),
145    fCutArray(NULL),
146    fConversionCuts(NULL),
147    fMesonCutArray(NULL),
148    fMesonCuts(NULL),
149    hESDConvGammaPt(NULL),
150    hESDMotherInvMassPt(NULL),
151    sESDMotherInvMassPtZM(NULL),
152    hESDMotherBackInvMassPt(NULL),
153    sESDMotherBackInvMassPtZM(NULL),
154    hESDMotherInvMassEalpha(NULL),
155    hMCAllGammaPt(NULL),
156    hMCDecayGammaPi0Pt(NULL),
157    hMCDecayGammaRhoPt(NULL),
158    hMCDecayGammaEtaPt(NULL),
159    hMCDecayGammaOmegaPt(NULL),
160    hMCDecayGammaEtapPt(NULL),
161    hMCDecayGammaPhiPt(NULL),
162    hMCConvGammaPt(NULL),
163    hMCConvGammaR(NULL),
164    hMCConvGammaEta(NULL),
165    hMCPi0Pt(NULL),
166    hMCEtaPt(NULL),
167    hMCPi0InAccPt(NULL),
168    hMCEtaInAccPt(NULL),
169    hESDTrueMotherInvMassPt(NULL),
170    hESDTruePrimaryMotherInvMassMCPt(NULL),
171    hESDTruePrimaryPi0ESDPtMCPt(NULL),
172    hESDTrueSecondaryMotherInvMassPt(NULL),
173    hESDTrueSecondaryMotherFromK0sInvMassPt(NULL),
174    hESDTrueBckGGInvMassPt(NULL),
175    hESDTrueBckContInvMassPt(NULL),
176    hESDTrueMotherDalitzInvMassPt(NULL),
177    hESDTrueConvGammaPt(NULL),
178    hESDTrueTwoElecCombPt(NULL),
179    hESDTrueTwoPionCombPt(NULL),
180    hESDTrueElecPionCombPt(NULL),
181    hESDTrueCombPt(NULL),
182    hESDTruePrimaryConvGammaPt(NULL),
183    hESDTruePrimaryConvGammaR(NULL),
184    hESDTruePrimaryConvGammaEta(NULL),
185    hESDTruePrimaryConvGammaESDPtMCPt(NULL),
186    hESDTrueSecondaryConvGammaPt(NULL),
187    hESDTrueSecondaryConvGammaFromK0sPt(NULL),
188    hESDTrueSecondaryConvGammaFromXFromK0sPt(NULL),
189    hNEvents(NULL),
190    hNGoodESDTracks(NULL),
191    hNV0Tracks(NULL),
192    fRandom(0),
193    fUnsmearedPx(NULL),
194    fUnsmearedPy(NULL),
195    fUnsmearedPz(NULL),
196    fUnsmearedE(NULL),
197    fIsHeavyIon(kFALSE),
198    fDoMesonAnalysis(kTRUE)
199
200 {
201    // Define input and output slots here
202    DefineInput(0, TChain::Class());
203    DefineOutput(1, TList::Class());
204 }
205
206 AliAnalysisTaskGammaConvV1::~AliAnalysisTaskGammaConvV1()
207 {
208    if(fGoodGammas){
209       delete fGoodGammas;
210       fGoodGammas = 0x0;
211    }
212    if(fBGHandler){
213       delete[] fBGHandler;
214       fBGHandler = 0x0;
215    }
216 }
217 //___________________________________________________________
218 void AliAnalysisTaskGammaConvV1::InitBack(){
219
220    Double_t *zBinLimitsArray = new Double_t[9];
221    zBinLimitsArray[0] = -50.00;
222    zBinLimitsArray[1] = -3.375;
223    zBinLimitsArray[2] = -1.605;
224    zBinLimitsArray[3] = -0.225;
225    zBinLimitsArray[4] = 1.065;
226    zBinLimitsArray[5] = 2.445;
227    zBinLimitsArray[6] = 4.245;
228    zBinLimitsArray[7] = 50.00;
229    zBinLimitsArray[8] = 1000.00;
230
231    Double_t *multiplicityBinLimitsArrayTracks = new Double_t[6];
232    multiplicityBinLimitsArrayTracks[0] = 0;
233    multiplicityBinLimitsArrayTracks[1] = 8.5;
234    multiplicityBinLimitsArrayTracks[2] = 16.5;
235    multiplicityBinLimitsArrayTracks[3] = 27.5;
236    multiplicityBinLimitsArrayTracks[4] = 41.5;
237    multiplicityBinLimitsArrayTracks[5] = 200.;
238    if(fIsHeavyIon){
239       multiplicityBinLimitsArrayTracks[0] = 0;
240       multiplicityBinLimitsArrayTracks[1] = 200.;
241       multiplicityBinLimitsArrayTracks[2] = 500.;
242       multiplicityBinLimitsArrayTracks[3] = 1000.;
243       multiplicityBinLimitsArrayTracks[4] = 1500.;
244       multiplicityBinLimitsArrayTracks[5] = 5000.;
245    }
246
247    Double_t *multiplicityBinLimitsArrayV0s = new Double_t[5];
248    multiplicityBinLimitsArrayV0s[0] = 2;
249    multiplicityBinLimitsArrayV0s[1] = 3;
250    multiplicityBinLimitsArrayV0s[2] = 4;
251    multiplicityBinLimitsArrayV0s[3] = 5;
252    multiplicityBinLimitsArrayV0s[4] = 9999;
253    if(fIsHeavyIon){
254       multiplicityBinLimitsArrayV0s[0] = 2;
255       multiplicityBinLimitsArrayV0s[1] = 10;
256       multiplicityBinLimitsArrayV0s[2] = 30;
257       multiplicityBinLimitsArrayV0s[3] = 50;
258       multiplicityBinLimitsArrayV0s[4] = 9999;
259    }
260
261    const Int_t nDim = 4;
262    Int_t nBins[nDim] = {1000,250,8,5};
263    Double_t xMin[nDim] = {0,0, 0,0};
264    Double_t xMax[nDim] = {1,25,8,5};
265
266    sESDMotherInvMassPtZM = new THnSparseF*[fnCuts];
267    sESDMotherBackInvMassPtZM = new THnSparseF*[fnCuts];
268
269    fBGHandler = new AliGammaConversionAODBGHandler*[fnCuts];
270    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
271       TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
272                 TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
273       fBackList[iCut] = new TList();
274       fBackList[iCut]->SetName(Form("%s_%s Back histograms",cutstring.Data(),cutstringMeson.Data()));
275       fBackList[iCut]->SetOwner(kTRUE);
276       fCutFolder[iCut]->Add(fBackList[iCut]);
277
278       sESDMotherBackInvMassPtZM[iCut] = new THnSparseF("Back_Back_InvMass_Pt_z_m","Back_Back_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
279       sESDMotherBackInvMassPtZM[iCut]->Sumw2();
280       fBackList[iCut]->Add(sESDMotherBackInvMassPtZM[iCut]);
281       
282       fMotherList[iCut] = new TList();
283       fMotherList[iCut]->SetName(Form("%s_%s Mother histograms",cutstring.Data(),cutstringMeson.Data()));
284       fMotherList[iCut]->SetOwner(kTRUE);
285       fCutFolder[iCut]->Add(fMotherList[iCut]);
286
287       sESDMotherInvMassPtZM[iCut] = new THnSparseF("Back_Mother_InvMass_Pt_z_m","Back_Mother_InvMass_Pt_z_m",nDim,nBins,xMin,xMax);
288       sESDMotherInvMassPtZM[iCut]->Sumw2();
289       fMotherList[iCut]->Add(sESDMotherInvMassPtZM[iCut]);
290
291       if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseTrackMultiplicity()){
292          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,6,((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->NumberOfBGEvents());
293          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayTracks);
294       }
295       else{
296          fBGHandler[iCut] = new AliGammaConversionAODBGHandler(9,5,((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->NumberOfBGEvents());
297          fBGHandler[iCut]->Initialize(zBinLimitsArray, multiplicityBinLimitsArrayV0s);
298       }
299    }
300 }
301 //________________________________________________________________________
302 void AliAnalysisTaskGammaConvV1::UserCreateOutputObjects()
303 {
304
305    // Create histograms
306    if(fOutputContainer != NULL){
307       delete fOutputContainer;
308       fOutputContainer = NULL;
309    }
310    if(fOutputContainer == NULL){
311       fOutputContainer = new TList();
312       fOutputContainer->SetOwner(kTRUE);
313    }
314
315    // Array of current cut's gammas
316    fGoodGammas = new TList();
317
318    fCutFolder = new TList*[fnCuts];
319    fESDList = new TList*[fnCuts];
320    fBackList = new TList*[fnCuts];
321    fMotherList = new TList*[fnCuts];
322    hESDConvGammaPt = new TH1F*[fnCuts];
323    hNEvents = new TH1I*[fnCuts];
324    hNGoodESDTracks = new TH1I*[fnCuts];
325    hNV0Tracks = new TH1I*[fnCuts];
326
327    if(fDoMesonAnalysis){
328       hESDMotherInvMassPt = new TH2F*[fnCuts];
329       hESDMotherBackInvMassPt = new TH2F*[fnCuts];
330       hESDMotherInvMassEalpha = new TH2F*[fnCuts];
331    }
332    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
333
334       TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
335                 TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
336                 
337       fCutFolder[iCut] = new TList();
338       fCutFolder[iCut]->SetName(Form("Cut Number %s_%s",cutstring.Data(),cutstringMeson.Data()));
339       fCutFolder[iCut]->SetOwner(kTRUE);
340       fOutputContainer->Add(fCutFolder[iCut]);
341       fESDList[iCut] = new TList();
342       fESDList[iCut]->SetName(Form("%s_%s ESD histograms",cutstring.Data(),cutstringMeson.Data()));
343       fESDList[iCut]->SetOwner(kTRUE);
344
345       hNEvents[iCut] = new TH1I("NEvents","NEvents",7,-0.5,6.5);
346       fESDList[iCut]->Add(hNEvents[iCut]);
347       if(fIsHeavyIon) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",3000,0,3000);
348       else hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
349       fESDList[iCut]->Add(hNGoodESDTracks[iCut]);
350       if(fIsHeavyIon) hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",25000,0,25000);
351       else hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2000,0,2000);
352       fESDList[iCut]->Add(hNV0Tracks[iCut]);
353
354       hESDConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
355       fESDList[iCut]->Add(hESDConvGammaPt[iCut]);
356
357       if(fDoMesonAnalysis){
358          hESDMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",1000,0,1,250,0,25);
359          fESDList[iCut]->Add(hESDMotherInvMassPt[iCut]);
360          hESDMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",1000,0,1,250,0,25);
361          fESDList[iCut]->Add(hESDMotherBackInvMassPt[iCut]);
362          hESDMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",1000,0,1,250,0,25);
363          fESDList[iCut]->Add(hESDMotherInvMassEalpha[iCut]);
364       }
365
366       fCutFolder[iCut]->Add(fESDList[iCut]);
367    }
368
369    if(fDoMesonAnalysis){
370       InitBack(); // Init Background Handler
371    }
372
373    if(MCEvent()){
374       // MC Histogramms
375       fMCList = new TList*[fnCuts];
376       // True Histogramms
377       fTrueList = new TList*[fnCuts];
378       // Selected Header List
379       fHeaderNameList = new TList*[fnCuts];
380
381       hMCAllGammaPt = new TH1F*[fnCuts];
382       hMCDecayGammaPi0Pt = new TH1F*[fnCuts];
383       hMCDecayGammaRhoPt = new TH1F*[fnCuts];
384       hMCDecayGammaEtaPt = new TH1F*[fnCuts];
385       hMCDecayGammaOmegaPt = new TH1F*[fnCuts];
386       hMCDecayGammaEtapPt = new TH1F*[fnCuts];
387       hMCDecayGammaPhiPt = new TH1F*[fnCuts];
388       hMCConvGammaPt = new TH1F*[fnCuts];
389       hMCConvGammaR = new TH1F*[fnCuts];
390       hMCConvGammaEta = new TH1F*[fnCuts];
391       hESDTrueConvGammaPt = new TH1F*[fnCuts];
392       hESDTrueTwoElecCombPt = new TH1F*[fnCuts];
393       hESDTrueTwoPionCombPt = new TH1F*[fnCuts];
394       hESDTrueElecPionCombPt = new TH1F*[fnCuts];
395       hESDTrueCombPt = new TH1F*[fnCuts];
396       hESDTruePrimaryConvGammaPt = new TH1F*[fnCuts];
397       hESDTruePrimaryConvGammaR = new TH1F*[fnCuts];
398       hESDTruePrimaryConvGammaEta = new TH1F*[fnCuts];
399       hESDTruePrimaryConvGammaESDPtMCPt = new TH2F*[fnCuts];
400       hESDTrueSecondaryConvGammaPt = new TH1F*[fnCuts];
401       hESDTrueSecondaryConvGammaFromK0sPt = new TH1F*[fnCuts];
402       hESDTrueSecondaryConvGammaFromXFromK0sPt = new TH1F*[fnCuts];
403
404       if(fDoMesonAnalysis){
405          hMCPi0Pt = new TH1F*[fnCuts];
406          hMCEtaPt = new TH1F*[fnCuts];
407          hMCPi0InAccPt = new TH1F*[fnCuts];
408          hMCEtaInAccPt = new TH1F*[fnCuts];
409
410          hESDTrueMotherInvMassPt = new TH2F*[fnCuts];
411          hESDTruePrimaryPi0ESDPtMCPt = new TH2F*[fnCuts];
412          hESDTruePrimaryMotherInvMassMCPt = new TH2F*[fnCuts];
413          hESDTrueSecondaryMotherInvMassPt = new TH2F*[fnCuts];
414          hESDTrueSecondaryMotherFromK0sInvMassPt = new TH2F*[fnCuts];
415          hESDTrueBckGGInvMassPt = new TH2F*[fnCuts];
416          hESDTrueBckContInvMassPt = new TH2F*[fnCuts];
417          hESDTrueMotherDalitzInvMassPt = new TH2F*[fnCuts];
418       }
419
420       for(Int_t iCut = 0; iCut<fnCuts;iCut++){
421          TString cutstring = ((AliConversionCuts*)fCutArray->At(iCut))->GetCutNumber();
422                         TString cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
423          
424                         fMCList[iCut] = new TList();
425          fMCList[iCut]->SetName(Form("%s_%s MC histograms",cutstring.Data(),cutstringMeson.Data()));
426          fMCList[iCut]->SetOwner(kTRUE);
427          fCutFolder[iCut]->Add(fMCList[iCut]);
428
429          hMCAllGammaPt[iCut] = new TH1F("MC_AllGamma_Pt","MC_AllGamma_Pt",250,0,25);
430          fMCList[iCut]->Add(hMCAllGammaPt[iCut]);
431          hMCDecayGammaPi0Pt[iCut] = new TH1F("MC_DecayGammaPi0_Pt","MC_DecayGammaPi0_Pt",250,0,25);
432          fMCList[iCut]->Add(hMCDecayGammaPi0Pt[iCut]);
433          hMCDecayGammaRhoPt[iCut] = new TH1F("MC_DecayGammaRho_Pt","MC_DecayGammaRho_Pt",250,0,25);
434          fMCList[iCut]->Add(hMCDecayGammaRhoPt[iCut]);
435          hMCDecayGammaEtaPt[iCut] = new TH1F("MC_DecayGammaEta_Pt","MC_DecayGammaEta_Pt",250,0,25);
436          fMCList[iCut]->Add(hMCDecayGammaEtaPt[iCut]);
437          hMCDecayGammaOmegaPt[iCut] = new TH1F("MC_DecayGammaOmega_Pt","MC_DecayGammaOmmega_Pt",250,0,25);
438          fMCList[iCut]->Add(hMCDecayGammaOmegaPt[iCut]);
439          hMCDecayGammaEtapPt[iCut] = new TH1F("MC_DecayGammaEtap_Pt","MC_DecayGammaEtap_Pt",250,0,25);
440          fMCList[iCut]->Add(hMCDecayGammaEtapPt[iCut]);
441          hMCDecayGammaPhiPt[iCut] = new TH1F("MC_DecayGammaPhi_Pt","MC_DecayGammaPhi_Pt",250,0,25);
442          fMCList[iCut]->Add(hMCDecayGammaPhiPt[iCut]);
443          hMCConvGammaPt[iCut] = new TH1F("MC_ConvGamma_Pt","MC_ConvGamma_Pt",250,0,25);
444          fMCList[iCut]->Add(hMCConvGammaPt[iCut]);
445          hMCConvGammaR[iCut] = new TH1F("MC_ConvGamma_R","MC_ConvGamma_R",1000,0,250);
446          fMCList[iCut]->Add(hMCConvGammaR[iCut]);
447          hMCConvGammaEta[iCut] = new TH1F("MC_ConvGamma_Eta","MC_ConvGamma_Eta",100,-4,4);
448          fMCList[iCut]->Add(hMCConvGammaEta[iCut]);
449          if(fDoMesonAnalysis){
450             hMCPi0Pt[iCut] = new TH1F("MC_Pi0_Pt","MC_Pi0_Pt",250,0,25);
451             fMCList[iCut]->Add(hMCPi0Pt[iCut]);
452             hMCEtaPt[iCut] = new TH1F("MC_Eta_Pt","MC_Eta_Pt",250,0,25);
453             fMCList[iCut]->Add(hMCEtaPt[iCut]);
454             hMCPi0InAccPt[iCut] = new TH1F("MC_Pi0InAcc_Pt","MC_Pi0InAcc_Pt",250,0,25);
455             fMCList[iCut]->Add(hMCPi0InAccPt[iCut]);
456             hMCEtaInAccPt[iCut] = new TH1F("MC_EtaInAcc_Pt","MC_EtaInAcc_Pt",250,0,25);
457             fMCList[iCut]->Add(hMCEtaInAccPt[iCut]);
458          }
459          fTrueList[iCut] = new TList();
460          fTrueList[iCut]->SetName(Form("%s_%s True histograms",cutstring.Data(),cutstringMeson.Data()));
461          fTrueList[iCut]->SetOwner(kTRUE);
462          fCutFolder[iCut]->Add(fTrueList[iCut]);
463
464          hESDTrueConvGammaPt[iCut] = new TH1F("ESD_TrueConvGamma_Pt","ESD_TrueConvGamma_Pt",250,0,25);
465          fTrueList[iCut]->Add(hESDTrueConvGammaPt[iCut]);
466          hESDTrueTwoElecCombPt[iCut] = new TH1F("ESD_TrueTwoElecComb_Pt","ESD_TrueTwoElecComb_Pt",250,0,25);
467          fTrueList[iCut]->Add(hESDTrueTwoElecCombPt[iCut]);
468          hESDTrueTwoPionCombPt[iCut] = new TH1F("ESD_TrueTwoPionComb_Pt","ESD_TrueTwoPionComb_Pt",250,0,25);
469          fTrueList[iCut]->Add(hESDTrueTwoPionCombPt[iCut]);
470          hESDTrueElecPionCombPt[iCut] = new TH1F("ESD_TrueElecPionComb_Pt","ESD_TrueElecPionComb_Pt",250,0,25);
471          fTrueList[iCut]->Add(hESDTrueElecPionCombPt[iCut]);
472          hESDTrueCombPt[iCut] = new TH1F("ESD_TrueComb_Pt","ESD_TrueComb_Pt",250,0,25);
473          fTrueList[iCut]->Add(hESDTrueCombPt[iCut]);
474          hESDTruePrimaryConvGammaPt[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Pt","ESD_TruePrimaryConvGamma_Pt",250,0,25);
475          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaPt[iCut]);
476          hESDTruePrimaryConvGammaR[iCut] = new TH1F("ESD_TruePrimaryConvGamma_R","ESD_TruePrimaryConvGamma_R",1000,0,250);
477          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaR[iCut]);
478          hESDTruePrimaryConvGammaEta[iCut] = new TH1F("ESD_TruePrimaryConvGamma_Eta","ESD_TruePrimaryConvGamma_Eta",100,-4,4);
479          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaEta[iCut]);
480          hESDTrueSecondaryConvGammaPt[iCut] = new TH1F("ESD_TrueSecondaryConvGamma_Pt","ESD_TrueSecondaryConvGamma_Pt",250,0,25);
481          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaPt[iCut]);
482          hESDTrueSecondaryConvGammaFromK0sPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromK0s_Pt","ESD_TrueSecondaryConvGammaFromK0s_Pt",250,0,25);
483          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaFromK0sPt[iCut]);
484          hESDTrueSecondaryConvGammaFromXFromK0sPt[iCut] = new TH1F("ESD_TrueSecondaryConvGammaFromXFromK0s_Pt", "ESD_TrueSecondaryConvGammaFromXFromK0s_Pt",250,0,25);
485          fTrueList[iCut]->Add(hESDTrueSecondaryConvGammaFromXFromK0sPt[iCut]);
486          hESDTruePrimaryConvGammaESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryConvGammaESD_PtMCPt", "ESD_TruePrimaryConvGammaESD_PtMCPt",250,0,25,250,0,25);
487          fTrueList[iCut]->Add(hESDTruePrimaryConvGammaESDPtMCPt[iCut]);
488
489          if(fDoMesonAnalysis){
490             hESDTrueMotherInvMassPt[iCut] = new TH2F("ESD_TrueMother_InvMass_Pt","ESD_TrueMother_InvMass_Pt",1000,0,1,250,0,25);
491             fTrueList[iCut]->Add(hESDTrueMotherInvMassPt[iCut]);
492             hESDTruePrimaryPi0ESDPtMCPt[iCut] = new TH2F("ESD_TruePrimaryPi0_ESDPt_MCPt","ESD_TruePrimaryPi0_ESDPt_MCPt",250,0,25,250,0,25);
493             fTrueList[iCut]->Add(hESDTruePrimaryPi0ESDPtMCPt[iCut]);
494             hESDTruePrimaryMotherInvMassMCPt[iCut] = new TH2F("ESD_TruePrimaryMother_InvMass_MCPt", "ESD_TruePrimaryMother_InvMass_MCPt", 1000,0,1,250,0,25);
495             fTrueList[iCut]->Add(hESDTruePrimaryMotherInvMassMCPt[iCut]);
496             hESDTrueSecondaryMotherInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMother_InvMass_Pt", "ESD_TrueSecondaryMother_InvMass_Pt", 1000,0,1,250,0,25);
497             fTrueList[iCut]->Add(hESDTrueSecondaryMotherInvMassPt[iCut]);
498             hESDTrueSecondaryMotherFromK0sInvMassPt[iCut] = new TH2F("ESD_TrueSecondaryMotherFromK0s_InvMass_Pt","ESD_TrueSecondaryMotherFromK0s_InvMass_Pt",1000,0,1,250,0,25);
499             fTrueList[iCut]->Add(hESDTrueSecondaryMotherFromK0sInvMassPt[iCut]);
500             hESDTrueBckGGInvMassPt[iCut] = new TH2F("ESD_TrueBckGG_InvMass_Pt","ESD_TrueBckGG_InvMass_Pt",1000,0,1,250,0,25);
501             fTrueList[iCut]->Add(hESDTrueBckGGInvMassPt[iCut]);
502             hESDTrueBckContInvMassPt[iCut] = new TH2F("ESD_TrueBckCont_InvMass_Pt","ESD_TrueBckCont_InvMass_Pt",1000,0,1,250,0,25);
503             fTrueList[iCut]->Add(hESDTrueBckContInvMassPt[iCut]);
504             hESDTrueMotherDalitzInvMassPt[iCut] = new TH2F("ESD_TrueDalitz_InvMass_Pt","ESD_TrueDalitz_InvMass_Pt",1000,0,1,250,0,25);
505             fTrueList[iCut]->Add(hESDTrueMotherDalitzInvMassPt[iCut]);
506          }
507          
508          if(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection() == 2){
509             fHeaderNameList[iCut] = new TList();
510             TString HeaderNames = "Header";
511             for(Int_t i = 0;i<(((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader())->GetEntries();i++){
512                HeaderNames = HeaderNames+"_"+ ((TObjString*)((TList*) ( (AliConversionCuts*)fCutArray->At(iCut)) ->GetAcceptedHeader())->At(i))->GetString();
513             }
514             fHeaderNameList[iCut]->SetName(HeaderNames);
515             fHeaderNameList[iCut]->SetOwner(kTRUE);
516             fCutFolder[iCut]->Add(fHeaderNameList[iCut]);
517          }
518       }
519    }
520
521    PostData(1, fOutputContainer);
522 }
523
524 //_____________________________________________________________________________
525 void AliAnalysisTaskGammaConvV1::UserExec(Option_t *)
526 {
527    //
528    // Called for each event
529    //
530    fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
531    if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
532
533    Int_t eventQuality = ((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetEventQuality();
534    if(eventQuality != 0){// Event Not Accepted
535       for(Int_t iCut = 0; iCut<fnCuts; iCut++){
536          hNEvents[iCut]->Fill(eventQuality);
537       }
538       return;
539    }
540
541    fMCEvent = MCEvent();
542    fESDEvent = (AliESDEvent*) InputEvent();
543    fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
544    CountESDTracks(); // Estimate Event Multiplicity
545
546    for(Int_t iCut = 0; iCut<fnCuts; iCut++){
547       fiCut = iCut;
548       if(fIsHeavyIon && !((AliConversionCuts*)fCutArray->At(iCut))->IsCentralitySelected(fESDEvent)){
549          hNEvents[iCut]->Fill(1); // Check Centrality --> Not Accepted => eventQuality = 1
550          continue;
551       }
552       hNEvents[iCut]->Fill(eventQuality);
553
554       hNGoodESDTracks[iCut]->Fill(fNumberOfESDTracks);
555       hNV0Tracks[iCut]->Fill(fESDEvent->GetVZEROData()->GetMTotV0A()+fESDEvent->GetVZEROData()->GetMTotV0C());
556       if(fMCEvent){ // Process MC Particle
557          fMCStack = fMCEvent->Stack();
558          if(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection() != 0){
559             ((AliConversionCuts*)fCutArray->At(iCut))->GetNotRejectedParticles(((AliConversionCuts*)fCutArray->At(iCut))->GetSignalRejection(),
560                                                                                ((AliConversionCuts*)fCutArray->At(iCut))->GetAcceptedHeader(),
561                                                                                fMCEvent);
562          }
563          ProcessMCParticles();
564       }
565       
566       ProcessPhotonCandidates(); // Process this cuts gammas
567
568       if(fDoMesonAnalysis){ // Meson Analysis
569          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fMCEvent){
570             fUnsmearedPx = new Double_t[fGoodGammas->GetEntries()]; // Store unsmeared Momenta
571             fUnsmearedPy = new Double_t[fGoodGammas->GetEntries()];
572             fUnsmearedPz = new Double_t[fGoodGammas->GetEntries()];
573             fUnsmearedE =  new Double_t[fGoodGammas->GetEntries()];
574
575             for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
576                fUnsmearedPx[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Px();
577                fUnsmearedPy[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Py();
578                fUnsmearedPz[gamma] = ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->Pz();
579                fUnsmearedE[gamma] =  ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->E();
580                ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->SmearParticle(dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(gamma)));
581             }
582          }
583
584          CalculatePi0Candidates(); // Combine Gammas
585          CalculateBackground(); // Combinatorial Background
586          UpdateEventByEventData(); // Store Event for mixed Events
587
588          if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->UseMCPSmearing() && fMCEvent){
589             for(Int_t gamma=0;gamma<fGoodGammas->GetEntries();gamma++){ // Smear the AODPhotons in MC
590                ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPx(fUnsmearedPx[gamma]); // Reset Unsmeared Momenta
591                ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPy(fUnsmearedPy[gamma]);
592                ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetPz(fUnsmearedPz[gamma]);
593                ((AliAODConversionPhoton*)fGoodGammas->At(gamma))->SetE(fUnsmearedE[gamma]);
594             }
595             delete[] fUnsmearedPx; fUnsmearedPx = 0x0;
596             delete[] fUnsmearedPy; fUnsmearedPy = 0x0;
597             delete[] fUnsmearedPz; fUnsmearedPz = 0x0;
598             delete[] fUnsmearedE;  fUnsmearedE  = 0x0;
599          }
600       }
601       fGoodGammas->Clear(); // delete this cuts good gammas
602
603    }
604    
605    PostData(1, fOutputContainer);
606 }
607 //________________________________________________________________________
608 void AliAnalysisTaskGammaConvV1::ProcessPhotonCandidates()
609 {
610    Int_t nV0 = 0;
611    TList *GoodGammasStepOne = new TList();
612    TList *GoodGammasStepTwo = new TList();
613    // Loop over Photon Candidates allocated by ReaderV1
614    for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
615       AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
616       if(!PhotonCandidate) continue;
617
618       if(fMCEvent && ((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
619          if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelPositive(), fMCStack)) continue;
620          if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(PhotonCandidate->GetMCLabelNegative(), fMCStack)) continue;
621       }
622
623       if(!((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fESDEvent)) continue;
624
625       if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
626          !((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // if no post reader loop is required add to events good gammas
627          fGoodGammas->Add(PhotonCandidate);
628          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
629          if(fMCEvent){
630             ProcessTruePhotonCandidates(PhotonCandidate);
631          }
632       }
633       else if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
634          ((AliConversionCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
635          nV0++;
636          GoodGammasStepOne->Add(PhotonCandidate);
637       }
638       else if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
639               ((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
640          GoodGammasStepTwo->Add(PhotonCandidate);
641       }
642    }
643    if(((AliConversionCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
644       for(Int_t i = 0;i<GoodGammasStepOne->GetEntries();i++){
645          AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GoodGammasStepOne->At(i);
646          if(!PhotonCandidate) continue;
647          if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GoodGammasStepOne->GetEntries())) continue;
648          if(!((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
649             fGoodGammas->Add(PhotonCandidate);
650             hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
651             if(fMCEvent){
652                ProcessTruePhotonCandidates(PhotonCandidate);
653             }
654          }
655          else GoodGammasStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
656       }
657    }
658    if(((AliConversionCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
659       for(Int_t i = 0;i<GoodGammasStepTwo->GetEntries();i++){
660          AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GoodGammasStepTwo->At(i);
661          if(!PhotonCandidate) continue;
662          if(!((AliConversionCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GoodGammasStepTwo,i)) continue;
663          fGoodGammas->Add(PhotonCandidate); // Add gamma to current cut TList
664          hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt()); // Differences to old V0Reader in p_t due to conversion KF->TLorentzVector
665          if(fMCEvent){
666             ProcessTruePhotonCandidates(PhotonCandidate);
667          }
668       }
669    }
670
671    delete GoodGammasStepOne;
672    GoodGammasStepOne = 0x0;
673    delete GoodGammasStepTwo;
674    GoodGammasStepTwo = 0x0;
675
676 }
677 //________________________________________________________________________
678 void AliAnalysisTaskGammaConvV1::ProcessTruePhotonCandidates(AliAODConversionPhoton *TruePhotonCandidate)
679 {
680    // Process True Photons
681    AliStack *MCStack = fMCEvent->Stack();
682    TParticle *posDaughter = TruePhotonCandidate->GetPositiveMCDaughter(MCStack);
683    TParticle *negDaughter = TruePhotonCandidate->GetNegativeMCDaughter(MCStack);
684
685    if(posDaughter == NULL || negDaughter == NULL) return; // One particle does not exist
686    if(posDaughter->GetMother(0) != negDaughter->GetMother(0)){  // Not Same Mother == Combinatorial Bck
687       if(TMath::Abs(posDaughter->GetPdgCode())==11 && TMath::Abs(negDaughter->GetPdgCode())==11)
688          hESDTrueTwoElecCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //Electron Combinatorial
689       else if(TMath::Abs(posDaughter->GetPdgCode())==211 && TMath::Abs(negDaughter->GetPdgCode())==211)
690          hESDTrueTwoPionCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
691       else if( (TMath::Abs(posDaughter->GetPdgCode())==11 && TMath::Abs(negDaughter->GetPdgCode())==211) ||
692                (TMath::Abs(posDaughter->GetPdgCode())==211 && TMath::Abs(negDaughter->GetPdgCode())==11) )
693          hESDTrueElecPionCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
694       else hESDTrueCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
695       return;
696    }
697    else if(posDaughter->GetMother(0) == -1){
698       if(TMath::Abs(posDaughter->GetPdgCode())==11 && TMath::Abs(negDaughter->GetPdgCode())==11)
699          hESDTrueTwoElecCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //Electron Combinatorial
700       else if(TMath::Abs(posDaughter->GetPdgCode())==211 && TMath::Abs(negDaughter->GetPdgCode())==211)
701          hESDTrueTwoPionCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
702       else if( (TMath::Abs(posDaughter->GetPdgCode())==11 && TMath::Abs(negDaughter->GetPdgCode())==211) ||
703                (TMath::Abs(posDaughter->GetPdgCode())==211 && TMath::Abs(negDaughter->GetPdgCode())==11) )
704          hESDTrueElecPionCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
705       else hESDTrueCombPt[fiCut]->Fill(TruePhotonCandidate->Pt()); //At least on Pion Combinatorial
706       return;
707    }
708    if(TMath::Abs(posDaughter->GetPdgCode())!=11 || TMath::Abs(negDaughter->GetPdgCode())!=11) return; //One Particle is not electron
709    if(posDaughter->GetPdgCode()==negDaughter->GetPdgCode()) return; // Same Charge
710    if(posDaughter->GetUniqueID() != 5 || negDaughter->GetUniqueID() !=5) return;// check if the daughters come from a conversion
711
712    TParticle *Photon = TruePhotonCandidate->GetMCParticle(MCStack);
713    if(Photon->GetPdgCode() != 22) return; // Mother is no Photon
714
715    // True Photon
716    hESDTrueConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
717
718    if(posDaughter->GetMother(0) <= MCStack->GetNprimary()){
719       // Count just primary MC Gammas as true --> For Ratio esdtruegamma / mcconvgamma
720       hESDTruePrimaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
721       hESDTruePrimaryConvGammaR[fiCut]->Fill(TruePhotonCandidate->GetConversionRadius());
722       hESDTruePrimaryConvGammaEta[fiCut]->Fill(TruePhotonCandidate->Eta());
723       hESDTruePrimaryConvGammaESDPtMCPt[fiCut]->Fill(TruePhotonCandidate->Pt(),Photon->Pt());
724    }
725    else{
726       hESDTrueSecondaryConvGammaPt[fiCut]->Fill(TruePhotonCandidate->Pt());
727       if(MCStack->Particle(Photon->GetMother(0))->GetPdgCode() == 310){
728          hESDTrueSecondaryConvGammaFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
729       }
730       if(MCStack->Particle(Photon->GetMother(0))->GetMother(0) > -1 &&
731          MCStack->Particle(MCStack->Particle(Photon->GetMother(0))->GetMother(0))->GetPdgCode() == 310){
732          hESDTrueSecondaryConvGammaFromXFromK0sPt[fiCut]->Fill(TruePhotonCandidate->Pt());
733       }
734    }
735 }
736 //________________________________________________________________________
737 void AliAnalysisTaskGammaConvV1::ProcessMCParticles()
738 {
739
740    // Loop over all primary MC particle
741    for(Int_t i = 0; i < fMCStack->GetNprimary(); i++) {
742       TParticle* particle = (TParticle *)fMCStack->Particle(i);
743       if (!particle) continue;
744
745       if(((AliConversionCuts*)fCutArray->At(fiCut))->GetSignalRejection() != 0){
746          if(!((AliConversionCuts*)fCutArray->At(fiCut))->IsParticleFromBGEvent(i, fMCStack)) continue;
747       }
748
749       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kFALSE)){
750          hMCAllGammaPt[fiCut]->Fill(particle->Pt()); // All MC Gamma
751          if(particle->GetMother(0) >-1){ // Meson Decay Gamma
752             switch(fMCStack->Particle(particle->GetMother(0))->GetPdgCode()){
753             case 111: // Pi0
754                hMCDecayGammaPi0Pt[fiCut]->Fill(particle->Pt());
755                break;
756             case 113: // Rho0
757                hMCDecayGammaRhoPt[fiCut]->Fill(particle->Pt());
758                break;
759             case 221: // Eta
760                hMCDecayGammaEtaPt[fiCut]->Fill(particle->Pt());
761                break;
762             case 223: // Omega
763                hMCDecayGammaOmegaPt[fiCut]->Fill(particle->Pt());
764                break;
765             case 331: // Eta'
766                hMCDecayGammaEtapPt[fiCut]->Fill(particle->Pt());
767                break;
768             case 333: // Phi
769                hMCDecayGammaPhiPt[fiCut]->Fill(particle->Pt());
770                break;
771             }
772          }
773       }
774       if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(particle,fMCStack,kTRUE)){
775          hMCConvGammaPt[fiCut]->Fill(particle->Pt());
776          hMCConvGammaR[fiCut]->Fill(((TParticle*)fMCStack->Particle(particle->GetFirstDaughter()))->R());
777          hMCConvGammaEta[fiCut]->Fill(particle->Eta());
778       } // Converted MC Gamma
779       if(fDoMesonAnalysis){
780          if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelectedMC(particle,fMCStack,kFALSE)){
781             if(particle->GetPdgCode() == 111)hMCPi0Pt[fiCut]->Fill(particle->Pt()); // All MC Pi0
782             if(particle->GetPdgCode() == 221)hMCEtaPt[fiCut]->Fill(particle->Pt()); // All MC Eta
783             // Check the acceptance for both gammas
784             if(particle->GetNDaughters() == 2){
785                TParticle* daughter0 = (TParticle*)fMCStack->Particle(particle->GetFirstDaughter());
786                TParticle* daughter1 = (TParticle*)fMCStack->Particle(particle->GetLastDaughter());
787                if(((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter0,fMCStack,kFALSE) &&
788                   ((AliConversionCuts*)fCutArray->At(fiCut))->PhotonIsSelectedMC(daughter1,fMCStack,kFALSE) ){
789                   if(particle->GetPdgCode() == 111)hMCPi0InAccPt[fiCut]->Fill(particle->Pt()); // MC Pi0 with gamma in acc
790                   if(particle->GetPdgCode() == 221)hMCEtaInAccPt[fiCut]->Fill(particle->Pt()); // MC Eta with gamma in acc
791                }
792             }
793          }
794       }
795    }
796 }
797 //________________________________________________________________________
798 void AliAnalysisTaskGammaConvV1::CalculatePi0Candidates(){
799
800    // Conversion Gammas
801    if(fGoodGammas->GetEntries()>1){
802       for(Int_t firstGammaIndex=0;firstGammaIndex<fGoodGammas->GetEntries()-1;firstGammaIndex++){
803          AliAODConversionPhoton *gamma0=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(firstGammaIndex));
804          for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fGoodGammas->GetEntries();secondGammaIndex++){
805             AliAODConversionPhoton *gamma1=dynamic_cast<AliAODConversionPhoton*>(fGoodGammas->At(secondGammaIndex));
806             //Check for same Electron ID
807             if(gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelPositive() ||
808                gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelNegative() ||
809                gamma0->GetTrackLabelNegative() == gamma1->GetTrackLabelPositive() ||
810                gamma0->GetTrackLabelPositive() == gamma1->GetTrackLabelNegative() ) continue;
811
812             AliAODConversionMother *pi0cand = new AliAODConversionMother(gamma0,gamma1);
813             pi0cand->SetLabels(firstGammaIndex,secondGammaIndex);
814
815             if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(pi0cand,kTRUE))){
816                hESDMotherInvMassPt[fiCut]->Fill(pi0cand->M(),pi0cand->Pt());
817                if(pi0cand->GetAlpha()<0.1){
818                   hESDMotherInvMassEalpha[fiCut]->Fill(pi0cand->M(),pi0cand->E());
819                }
820                Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
821                Int_t mbin = 0;
822                if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
823                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
824                } else {
825                   mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
826                }
827                Double_t sparesFill[4] = {pi0cand->M(),pi0cand->Pt(),zbin,mbin};
828                sESDMotherInvMassPtZM[fiCut]->Fill(sparesFill,1);
829                if(fMCEvent){
830                   ProcessTrueMesonCandidates(pi0cand,gamma0,gamma1);
831                }
832             }
833             delete pi0cand;
834             pi0cand=0x0;
835          }
836       }
837    }
838 }
839 //______________________________________________________________________
840 void AliAnalysisTaskGammaConvV1::ProcessTrueMesonCandidates(AliAODConversionMother *Pi0Candidate, AliAODConversionPhoton *TrueGammaCandidate0, AliAODConversionPhoton *TrueGammaCandidate1)
841 {
842    // Process True Mesons
843    AliStack *MCStack = fMCEvent->Stack();
844
845    if(TrueGammaCandidate0->GetV0Index()<fESDEvent->GetNumberOfV0s()){
846       Bool_t isTruePi0 = kFALSE;
847       Bool_t isTrueEta = kFALSE;
848       Int_t gamma0MCLabel = TrueGammaCandidate0->GetMCParticleLabel(MCStack);
849       Int_t gamma0MotherLabel = -1;
850       if(gamma0MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
851          // Daughters Gamma 0
852          TParticle * negativeMC = (TParticle*)TrueGammaCandidate0->GetNegativeMCDaughter(MCStack);
853          TParticle * positiveMC = (TParticle*)TrueGammaCandidate0->GetPositiveMCDaughter(MCStack);
854          TParticle * gammaMC0 = (TParticle*)MCStack->Particle(gamma0MCLabel);
855          if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
856             if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
857                if(gammaMC0->GetPdgCode() == 22){ // ... with Gamma Mother
858                   gamma0MotherLabel=gammaMC0->GetFirstMother();
859                }
860             }
861             if(gammaMC0->GetPdgCode() ==111){ // Conversion but Pi0 Mother
862                gamma0MotherLabel=-111;
863             }
864             if(gammaMC0->GetPdgCode() ==221){ // Conversion but Eta Mother
865                gamma0MotherLabel=-221;
866             }
867          }
868       }
869       if(TrueGammaCandidate1->GetV0Index()<fESDEvent->GetNumberOfV0s()){
870          Int_t gamma1MCLabel = TrueGammaCandidate1->GetMCParticleLabel(MCStack);
871          Int_t gamma1MotherLabel = -1;
872          if(gamma1MCLabel != -1){ // Gamma is Combinatorial; MC Particles don't belong to the same Mother
873             // Daughters Gamma 1
874             TParticle * negativeMC = (TParticle*)TrueGammaCandidate1->GetNegativeMCDaughter(MCStack);
875             TParticle * positiveMC = (TParticle*)TrueGammaCandidate1->GetPositiveMCDaughter(MCStack);
876             TParticle * gammaMC1 = (TParticle*)MCStack->Particle(gamma1MCLabel);
877             if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){  // Electrons ...
878                if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){ // ... From Conversion ...
879                   if(gammaMC1->GetPdgCode() == 22){ // ... with Gamma Mother
880                      gamma1MotherLabel=gammaMC1->GetFirstMother();
881                   }
882                }
883                if(gammaMC1->GetPdgCode() ==111){ // Conversion but Pi0 Mother
884                   gamma1MotherLabel=-111;
885                }
886                if(gammaMC1->GetPdgCode() ==221){ // Conversion but Eta Mother
887                   gamma1MotherLabel=-221;
888                }
889             }
890          }
891          if(gamma0MotherLabel>=0 && gamma0MotherLabel==gamma1MotherLabel){
892             if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 111){
893                isTruePi0=kTRUE;
894             }
895             if(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetPdgCode() == 221){
896                isTrueEta=kTRUE;
897             }
898          }
899          if(isTruePi0 || isTrueEta){ // True Pion or Eta
900             hESDTrueMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
901             if(gamma0MotherLabel > MCStack->GetNprimary()){ // Secondary Meson
902                hESDTrueSecondaryMotherInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
903                if (((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0) >-1){
904                   if(MCStack->Particle(((TParticle*)MCStack->Particle(gamma1MotherLabel))->GetMother(0))->GetPdgCode()==kK0Short){
905                      hESDTrueSecondaryMotherFromK0sInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
906                   }
907                }
908             }
909             if(gamma0MotherLabel <= MCStack->GetNprimary()){ // Only primary pi0 for efficiency calculation
910                hESDTruePrimaryMotherInvMassMCPt[fiCut]->Fill(Pi0Candidate->M(),((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt());
911                if(isTruePi0){ // Only primaries for unfolding
912                   hESDTruePrimaryPi0ESDPtMCPt[fiCut]->Fill(Pi0Candidate->Pt(),((TParticle*)MCStack->Particle(gamma1MotherLabel))->Pt());
913                }
914             }
915          }
916          if(!isTruePi0 && !isTrueEta){ // Background
917             if(gamma0MotherLabel>-1 && gamma1MotherLabel>-1){ // Both Tracks are Photons and have a mother but not Pi0 or Eta
918                hESDTrueBckGGInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
919             } else { // No photon or without mother
920                hESDTrueBckContInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
921             }
922             if((gamma0MotherLabel==-111 || gamma1MotherLabel==-111 || gamma0MotherLabel==-221 || gamma1MotherLabel==-221) ){
923                // Dalitz
924                hESDTrueMotherDalitzInvMassPt[fiCut]->Fill(Pi0Candidate->M(),Pi0Candidate->Pt());
925             }
926          }
927       }
928    }
929 }
930
931 //________________________________________________________________________
932 void AliAnalysisTaskGammaConvV1::CalculateBackground(){
933
934    Int_t zbin= fBGHandler[fiCut]->GetZBinIndex(fESDEvent->GetPrimaryVertex()->GetZ());
935    Int_t mbin = 0;
936
937    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
938       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fNumberOfESDTracks);
939    } else {
940       mbin = fBGHandler[fiCut]->GetMultiplicityBinIndex(fGoodGammas->GetEntries());
941    }
942
943    if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseRotationMethod()){
944
945       for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntries();iCurrent++){
946          AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
947          for(Int_t iCurrent2=iCurrent+1;iCurrent2<fGoodGammas->GetEntries();iCurrent2++){
948             for(Int_t nRandom=0;nRandom<((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NumberOfBGEvents();nRandom++){
949                AliAODConversionPhoton currentEventGoodV02 = *(AliAODConversionPhoton*)(fGoodGammas->At(iCurrent2));
950
951                if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->DoBGProbability()){
952                   AliAODConversionMother *backgroundCandidateProb = new AliAODConversionMother(&currentEventGoodV0,&currentEventGoodV02);
953                   Double_t massBGprob = backgroundCandidateProb->M();
954                   if(massBGprob>0.1 && massBGprob<0.14){
955                      if(fRandom.Rndm()>fBGHandler[fiCut]->GetBGProb(zbin,mbin)){
956                         delete backgroundCandidateProb;
957                         continue;
958                      }
959                   }
960                   delete backgroundCandidateProb;
961                   backgroundCandidateProb = 0x0;
962                }
963
964                RotateParticle(&currentEventGoodV02);
965                AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&currentEventGoodV02);
966                if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
967                   hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
968                   Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
969                   sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
970                }
971                delete backgroundCandidate;
972                backgroundCandidate = 0x0;
973             }
974          }
975       }
976    }else{
977       AliGammaConversionAODBGHandler::GammaConversionVertex *bgEventVertex = NULL;
978
979       if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
980          for(Int_t nEventsInBG=0;nEventsInBG<fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
981             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
982             if(fMoveParticleAccordingToVertex == kTRUE){
983                bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
984             }
985
986             for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntries();iCurrent++){
987                AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
988                for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
989                   AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
990                   if(fMoveParticleAccordingToVertex == kTRUE){
991                      MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
992                   }
993
994                   AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
995                   if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
996                      hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
997                      Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
998                      sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
999                   }
1000                   delete backgroundCandidate;
1001                   backgroundCandidate = 0x0;
1002                }
1003             }
1004          }
1005       }
1006       else{
1007          for(Int_t nEventsInBG=0;nEventsInBG <fBGHandler[fiCut]->GetNBGEvents();nEventsInBG++){
1008             AliGammaConversionAODVector *previousEventV0s = fBGHandler[fiCut]->GetBGGoodV0s(zbin,mbin,nEventsInBG);
1009             if(previousEventV0s){
1010                if(fMoveParticleAccordingToVertex == kTRUE){
1011                   bgEventVertex = fBGHandler[fiCut]->GetBGEventVertex(zbin,mbin,nEventsInBG);
1012                }
1013                for(Int_t iCurrent=0;iCurrent<fGoodGammas->GetEntries();iCurrent++){
1014                   AliAODConversionPhoton currentEventGoodV0 = *(AliAODConversionPhoton*)(fGoodGammas->At(iCurrent));
1015                   for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1016
1017                      AliAODConversionPhoton previousGoodV0 = (AliAODConversionPhoton)(*(previousEventV0s->at(iPrevious)));
1018
1019                      if(fMoveParticleAccordingToVertex == kTRUE){
1020                         MoveParticleAccordingToVertex(&previousGoodV0,bgEventVertex);
1021                      }
1022
1023                      AliAODConversionMother *backgroundCandidate = new AliAODConversionMother(&currentEventGoodV0,&previousGoodV0);
1024
1025                      if((((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->MesonIsSelected(backgroundCandidate,kFALSE))){
1026                         hESDMotherBackInvMassPt[fiCut]->Fill(backgroundCandidate->M(),backgroundCandidate->Pt());
1027                         Double_t sparesFill[4] = {backgroundCandidate->M(),backgroundCandidate->Pt(),zbin,mbin};
1028                         sESDMotherBackInvMassPtZM[fiCut]->Fill(sparesFill,1);
1029                      }
1030                      delete backgroundCandidate;
1031                      backgroundCandidate = 0x0;
1032                   }
1033                }
1034             }
1035          }
1036       }
1037    }
1038 }
1039
1040 //________________________________________________________________________
1041 void AliAnalysisTaskGammaConvV1::RotateParticle(AliAODConversionPhoton *gamma){
1042    Int_t fNDegreesPMBackground= ((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->NDegreesRotation();
1043    Double_t nRadiansPM = fNDegreesPMBackground*TMath::Pi()/180;
1044    Double_t rotationValue = fRandom.Rndm()*2*nRadiansPM + TMath::Pi()-nRadiansPM;
1045    gamma->RotateZ(rotationValue);
1046 }
1047 //________________________________________________________________________
1048 void AliAnalysisTaskGammaConvV1::MoveParticleAccordingToVertex(AliAODConversionPhoton* particle,const AliGammaConversionAODBGHandler::GammaConversionVertex *vertex){
1049    //see header file for documentation
1050
1051    Double_t dx = vertex->fX - fESDEvent->GetPrimaryVertex()->GetX();
1052    Double_t dy = vertex->fY - fESDEvent->GetPrimaryVertex()->GetY();
1053    Double_t dz = vertex->fZ - fESDEvent->GetPrimaryVertex()->GetZ();
1054
1055    Double_t movedPlace[3] = {particle->GetConversionX() - dx,particle->GetConversionY() - dy,particle->GetConversionZ() - dz};
1056    particle->SetConversionPoint(movedPlace);
1057 }
1058 //________________________________________________________________________
1059 void AliAnalysisTaskGammaConvV1::UpdateEventByEventData(){
1060    //see header file for documentation
1061    if(fGoodGammas->GetEntries() >0 ){
1062       if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
1063          fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fNumberOfESDTracks);
1064       }
1065       else{ // means we use #V0s for multiplicity
1066          fBGHandler[fiCut]->AddEvent(fGoodGammas,fESDEvent->GetPrimaryVertex()->GetX(),fESDEvent->GetPrimaryVertex()->GetY(),fESDEvent->GetPrimaryVertex()->GetZ(),fGoodGammas->GetEntries());
1067       }
1068    }
1069 }
1070
1071 //________________________________________________________________________
1072 void AliAnalysisTaskGammaConvV1::CountESDTracks(){
1073
1074    AliESDtrackCuts *EsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
1075    // Using standard function for setting Cuts
1076    Bool_t selectPrimaries=kTRUE;
1077    EsdTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(selectPrimaries);
1078    EsdTrackCuts->SetMaxDCAToVertexZ(2);
1079    EsdTrackCuts->SetEtaRange(-0.8, 0.8);
1080    EsdTrackCuts->SetPtRange(0.15);
1081
1082    fNumberOfESDTracks = 0;
1083    for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1084       AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1085       if(!curTrack) continue;
1086       if(EsdTrackCuts->AcceptTrack(curTrack) ) fNumberOfESDTracks++;
1087    }
1088    delete EsdTrackCuts;
1089    EsdTrackCuts=0x0;
1090
1091    return;
1092 }
1093 //________________________________________________________________________
1094 void AliAnalysisTaskGammaConvV1::Terminate(const Option_t *)
1095 {
1096    fOutputContainer->Add(((AliConversionCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
1097
1098    for(Int_t iCut = 0; iCut<fnCuts;iCut++){
1099       if(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms()){
1100          fCutFolder[iCut]->Add(((AliConversionCuts*)fCutArray->At(iCut))->GetCutHistograms());
1101       }
1102       if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
1103          fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
1104       }
1105    }
1106    fOutputContainer->Print();
1107 }