covertity fix 24879
[u/mrichter/AliRoot.git] / PWGGA / GammaConv / AliAnalysisTaskGammaConvFlow.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: Andrea Dubla, Redmer Alexander Bertens, Friederike Bock            *
5  *                                                                        *
6  * Permission to use, copy, modify and distribute this software and its   *
7  * documentation strictly for non-commercial purposes is hereby granted   *
8  * without fee, provided that the above copyright notice appears in all   *
9  * copies and that both the copyright notice and this permission notice   *
10  * appear in the supporting documentation. The authors make no claims     *
11  * about the suitability of this software for any purpose. It is              *
12  * provided "as is" without express or implied warranty.                      *
13  **************************************************************************/
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 "TBranch.h"
24 #include "TFile.h"
25 #include "TH1F.h"
26 #include "TH2F.h"
27 #include "TH3F.h"
28 #include "THnSparse.h"
29 #include "TCanvas.h"
30 #include "TNtuple.h"
31 #include "AliAnalysisTask.h"
32 #include "AliAnalysisManager.h"
33 #include "AliESDEvent.h"
34 #include "AliESDInputHandler.h"
35 #include "AliMCEventHandler.h"
36 #include "AliMCEvent.h"
37 #include "AliMCParticle.h"
38 #include "AliCentrality.h"
39 #include "AliESDVZERO.h"
40 #include "AliESDpid.h"
41 #include "AliAnalysisTaskGammaConvFlow.h"
42 #include "AliVParticle.h"
43 #include "AliESDtrack.h"
44 #include "AliESDtrackCuts.h"
45 #include "AliKFVertex.h"
46 #include "AliV0ReaderV1.h"
47 #include "AliGenCocktailEventHeader.h"
48 #include "AliConversionAODBGHandlerRP.h"
49 #include "AliAODMCParticle.h"
50 #include "AliAODMCHeader.h"
51 #include "AliEventplane.h"
52
53 #include "AliFlowCandidateTrack.h"
54 #include "AliFlowTrackCuts.h"
55 #include "AliFlowEventSimple.h"
56 #include "AliFlowCommonConstants.h"
57 #include "AliFlowEvent.h"
58 #include "AliFlowTrack.h"
59
60 class AliFlowTrackCuts;
61
62 using namespace std;
63
64 ClassImp(AliAnalysisTaskGammaConvFlow)
65
66 //________________________________________________________________________
67 AliAnalysisTaskGammaConvFlow::AliAnalysisTaskGammaConvFlow(): AliAnalysisTaskSE(),
68 fV0Reader(NULL),
69 fBGHandler(NULL),
70 fBGHandlerRP(NULL),
71 fInputEvent(NULL),
72
73 fCutFolder(NULL),
74 fESDList(NULL),
75 fBackList(NULL),
76 fMotherList(NULL),
77 fPhotonDCAList(NULL),
78 fMesonDCAList(NULL),
79 //fTrueList(NULL),
80 //fMCList(NULL),
81 fHeaderNameList(NULL),
82 fOutputContainer(0),
83 fReaderGammas(NULL),
84 fGammaCandidates(NULL),
85 fEventCutArray(NULL),
86 fEventCuts(NULL),
87 fCutArray(NULL),
88 fConversionCuts(NULL),
89 fMesonCutArray(NULL),
90 fMesonCuts(NULL),
91 hESDConvGammaPt(NULL),
92 hESDConvGammaR(NULL),
93 hESDConvGammaEta(NULL),
94 //tESDConvGammaPtDcazCat(NULL),
95 fPtGamma(0),
96 fDCAzPhoton(0),
97 fRConvPhoton(0),
98 fEtaPhoton(0),
99 iCatPhoton(0),
100 iPhotonMCInfo(0),
101 hESDMotherInvMassPt(NULL),
102 //sESDMotherInvMassPtZM(NULL),
103 hESDMotherBackInvMassPt(NULL),
104 //sESDMotherBackInvMassPtZM(NULL),
105 hESDMotherInvMassEalpha(NULL),
106 hESDMotherPi0PtY(NULL),
107 hESDMotherEtaPtY(NULL),
108 hESDMotherPi0PtAlpha(NULL),
109 hESDMotherEtaPtAlpha(NULL),
110 hESDMotherPi0PtOpenAngle(NULL),
111 hESDMotherEtaPtOpenAngle(NULL),
112
113 hNEvents(NULL),
114 hNGoodESDTracks(NULL),
115 hNGammaCandidates(NULL),
116 hNGoodESDTracksVsNGammaCanditates(NULL),
117 hNV0Tracks(NULL),
118 hEtaShift(NULL),
119 tESDMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
120 fInvMass(0),
121 fPt(0),
122 fDCAzGammaMin(0),
123 fDCAzGammaMax(0),
124 iFlag(0),
125 iMesonMCInfo(0),
126 fEventPlaneAngle(-100),
127 fRandom(0),
128 fnGammaCandidates(0),
129 fUnsmearedPx(NULL),
130 fUnsmearedPy(NULL),
131 fUnsmearedPz(NULL),
132 fUnsmearedE(NULL),
133 fMCStackPos(NULL),
134 fMCStackNeg(NULL),
135 fESDArrayPos(NULL),
136 fESDArrayNeg(NULL),
137 fnCuts(0),
138 fiCut(0),
139 fMoveParticleAccordingToVertex(kTRUE),
140 fIsHeavyIon(0),
141 fDoMesonAnalysis(kTRUE),
142 fDoMesonQA(0),
143 fDoPhotonQA(0),
144 fIsFromMBHeader(kTRUE),
145 fhistoEPVZ(NULL),
146 fDebug(0),
147 fCutsRP(0),
148 fNullCuts(0), 
149 fFlowEvent(0)
150 {
151         // DefineOutput(1, TList::Class());
152         // DefineOutput(2, AliFlowEventSimple::Class());
153 }
154
155 //________________________________________________________________________
156 AliAnalysisTaskGammaConvFlow::AliAnalysisTaskGammaConvFlow(const char *name):
157 AliAnalysisTaskSE(name),
158 fV0Reader(NULL),
159 fBGHandler(NULL),
160 fBGHandlerRP(NULL),
161 fInputEvent(NULL),
162 fCutFolder(NULL),
163 fESDList(NULL),
164 fBackList(NULL),
165 fMotherList(NULL),
166 fPhotonDCAList(NULL),
167 fMesonDCAList(NULL),
168 //fTrueList(NULL),
169 //fMCList(NULL),
170 fHeaderNameList(NULL),
171 fOutputContainer(0),
172 fReaderGammas(NULL),
173 fGammaCandidates(NULL),
174 fEventCutArray(NULL),
175 fEventCuts(NULL),
176 fCutArray(NULL),
177 fConversionCuts(NULL),
178 fMesonCutArray(NULL),
179 fMesonCuts(NULL),
180 hESDConvGammaPt(NULL),
181 hESDConvGammaR(NULL),
182 hESDConvGammaEta(NULL),
183 //tESDConvGammaPtDcazCat(NULL),
184 fPtGamma(0),
185 fDCAzPhoton(0),
186 fRConvPhoton(0),
187 fEtaPhoton(0),
188 iCatPhoton(0),
189 iPhotonMCInfo(0),
190 hESDMotherInvMassPt(NULL),
191 //sESDMotherInvMassPtZM(NULL),
192 hESDMotherBackInvMassPt(NULL),
193 //sESDMotherBackInvMassPtZM(NULL),
194 hESDMotherInvMassEalpha(NULL),
195 hESDMotherPi0PtY(NULL),
196 hESDMotherEtaPtY(NULL),
197 hESDMotherPi0PtAlpha(NULL),
198 hESDMotherEtaPtAlpha(NULL),
199 hESDMotherPi0PtOpenAngle(NULL),
200 hESDMotherEtaPtOpenAngle(NULL),
201 hNEvents(NULL),
202 hNGoodESDTracks(NULL),
203 hNGammaCandidates(NULL),
204 hNGoodESDTracksVsNGammaCanditates(NULL),
205 hNV0Tracks(NULL),
206 hEtaShift(NULL),
207 tESDMesonsInvMassPtDcazMinDcazMaxFlag(NULL),
208 fInvMass(0),
209 fPt(0),
210 fDCAzGammaMin(0),
211 fDCAzGammaMax(0),
212 iFlag(0),
213 iMesonMCInfo(0),
214 fEventPlaneAngle(-100),
215 fRandom(0),
216 fnGammaCandidates(0),
217 fUnsmearedPx(NULL),
218 fUnsmearedPy(NULL),
219 fUnsmearedPz(NULL),
220 fUnsmearedE(NULL),
221 fMCStackPos(NULL),
222 fMCStackNeg(NULL),
223 fESDArrayPos(NULL),
224 fESDArrayNeg(NULL),
225 fnCuts(0),
226 fiCut(0),
227 fMoveParticleAccordingToVertex(kTRUE),
228 fIsHeavyIon(0),
229 fDoMesonAnalysis(kTRUE),
230 fDoMesonQA(0),
231 fDoPhotonQA(0),
232 fIsFromMBHeader(kTRUE),
233 fhistoEPVZ(NULL),
234 fDebug(0),
235 fCutsRP(0), 
236 fNullCuts(0), 
237 fFlowEvent(0)
238
239 {
240     // Define output slots here
241     DefineOutput(1, TList::Class());
242     DefineOutput(2, AliFlowEventSimple::Class());
243
244 }
245 //___________________________________________________________
246 AliAnalysisTaskGammaConvFlow::~AliAnalysisTaskGammaConvFlow()
247 {
248         if(fGammaCandidates){
249                 delete fGammaCandidates;
250                 fGammaCandidates = 0x0;
251         }
252         if(fBGHandler){
253                 delete[] fBGHandler;
254                 fBGHandler = 0x0;
255         }
256         if(fBGHandlerRP){
257                 delete[] fBGHandlerRP;
258                 fBGHandlerRP = 0x0;
259         }
260         
261         if (fFlowEvent) delete fFlowEvent;
262
263 }
264
265 //________________________________________________________________________
266 void AliAnalysisTaskGammaConvFlow::UserCreateOutputObjects(){
267         
268         // Create histograms
269         if(fOutputContainer != NULL){
270                 delete fOutputContainer;
271                 fOutputContainer = NULL;
272         }
273         if(fOutputContainer == NULL){
274                 fOutputContainer = new TList();
275                 fOutputContainer->SetOwner(kTRUE);
276         }
277         
278         //========================= again flow setting==========================
279         //Create histograms
280         //----------hfe initialising begin---------
281         fNullCuts = new AliFlowTrackCuts("null_cuts");
282         
283         AliFlowCommonConstants* cc = AliFlowCommonConstants::GetMaster();
284         cc->SetNbinsMult(10000);
285         cc->SetMultMin(0);
286         cc->SetMultMax(10000);
287         
288         cc->SetNbinsPt(100);
289         cc->SetPtMin(0);
290         cc->SetPtMax(20);
291         
292         cc->SetNbinsPhi(180);
293         cc->SetPhiMin(0.0);
294         cc->SetPhiMax(TMath::TwoPi());
295         
296         cc->SetNbinsEta(30);
297         cc->SetEtaMin(-8.0);
298         cc->SetEtaMax(+8.0);
299         
300         cc->SetNbinsQ(500);
301         cc->SetQMin(0.0);
302         cc->SetQMax(3.0);
303
304         // Array of current cut's gammas
305         fGammaCandidates = new TList();
306         
307         fCutFolder = new TList*[fnCuts];
308         fESDList = new TList*[fnCuts];
309         fBackList = new TList*[fnCuts];
310         fMotherList = new TList*[fnCuts];
311         hNEvents = new TH1I*[fnCuts];
312         hNGoodESDTracks = new TH1I*[fnCuts];
313         hNGammaCandidates = new TH1I*[fnCuts];
314         hNGoodESDTracksVsNGammaCanditates = new TH2F*[fnCuts];
315         hNV0Tracks = new TH1I*[fnCuts];
316         hEtaShift = new TProfile*[fnCuts];
317         hESDConvGammaPt = new TH1F*[fnCuts];
318         
319         if (fDoPhotonQA == 2){
320                 fPhotonDCAList = new TList*[fnCuts];
321 //      tESDConvGammaPtDcazCat = new TTree*[fnCuts];
322         }
323         if (fDoPhotonQA > 0){
324                 hESDConvGammaR = new TH1F*[fnCuts];
325                 hESDConvGammaEta = new TH1F*[fnCuts];
326         }
327         
328         if(fDoMesonAnalysis){
329                 hESDMotherInvMassPt = new TH2F*[fnCuts];
330                 hESDMotherBackInvMassPt = new TH2F*[fnCuts];
331                 hESDMotherInvMassEalpha = new TH2F*[fnCuts];
332                 if (fDoMesonQA == 2){
333                         fMesonDCAList = new TList*[fnCuts];
334                         tESDMesonsInvMassPtDcazMinDcazMaxFlag = new TTree*[fnCuts];
335                 }
336                 if (fDoMesonQA > 0){
337                         hESDMotherPi0PtY =  new TH2F*[fnCuts];
338                         hESDMotherEtaPtY =  new TH2F*[fnCuts];
339                         hESDMotherPi0PtAlpha =  new TH2F*[fnCuts];
340                         hESDMotherEtaPtAlpha =  new TH2F*[fnCuts];
341                         hESDMotherPi0PtOpenAngle =  new TH2F*[fnCuts];
342                         hESDMotherEtaPtOpenAngle =  new TH2F*[fnCuts];
343                 }
344         }
345         
346         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
347                 
348                 TString cutstringEvent  = ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber();
349                 TString cutstringPhoton = ((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutNumber();
350                 TString cutstringMeson = "NoMesonCut";
351                 if(fDoMesonAnalysis)cutstringMeson = ((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutNumber();
352                 
353                 fCutFolder[iCut] = new TList();
354                 fCutFolder[iCut]->SetName(Form("Cut Number %s_%s_%s",cutstringEvent.Data() ,cutstringPhoton.Data(),cutstringMeson.Data()));
355                 fCutFolder[iCut]->SetOwner(kTRUE);
356                 fOutputContainer->Add(fCutFolder[iCut]);
357                 fESDList[iCut] = new TList();
358                 fESDList[iCut]->SetName(Form("%s_%s_%s ESD histograms",cutstringEvent.Data() ,cutstringPhoton.Data(),cutstringMeson.Data()));
359                 fESDList[iCut]->SetOwner(kTRUE);
360                 fCutFolder[iCut]->Add(fESDList[iCut]);
361                 
362                 hNEvents[iCut] = new TH1I("NEvents","NEvents",10,-0.5,9.5);
363                 hNEvents[iCut]->GetXaxis()->SetBinLabel(1,"Accepted");
364                 hNEvents[iCut]->GetXaxis()->SetBinLabel(2,"Centrality");
365                 hNEvents[iCut]->GetXaxis()->SetBinLabel(3,"Missing MC");
366                 if (((AliConvEventCuts*)fEventCutArray->At(iCut))->IsSpecialTrigger() > 1 ){
367                         TString TriggerNames = "Not Trigger: ";
368                         TriggerNames = TriggerNames+ ( (AliConvEventCuts*)fEventCutArray->At(iCut))->GetSpecialTriggerName();
369                         hNEvents[iCut]->GetXaxis()->SetBinLabel(4,TriggerNames.Data());
370                 } else {
371                         hNEvents[iCut]->GetXaxis()->SetBinLabel(4,"Trigger");
372                 }
373                 hNEvents[iCut]->GetXaxis()->SetBinLabel(5,"Vertex Z");
374                 hNEvents[iCut]->GetXaxis()->SetBinLabel(6,"Cont. Vertex");
375                 hNEvents[iCut]->GetXaxis()->SetBinLabel(7,"Pile-Up");
376                 hNEvents[iCut]->GetXaxis()->SetBinLabel(8,"no SDD");
377                 hNEvents[iCut]->GetXaxis()->SetBinLabel(9,"no V0AND");
378                 hNEvents[iCut]->GetXaxis()->SetBinLabel(10,"EMCAL problem");
379                 fESDList[iCut]->Add(hNEvents[iCut]);
380                 
381                 if(fIsHeavyIon == 1) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",4000,0,4000);
382                 else if(fIsHeavyIon == 2) hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",400,0,400);
383                 else hNGoodESDTracks[iCut] = new TH1I("GoodESDTracks","GoodESDTracks",200,0,200);
384                 fESDList[iCut]->Add(hNGoodESDTracks[iCut]);
385                 if(fIsHeavyIon == 1) hNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",100,0,100);
386                 else if(fIsHeavyIon == 2) hNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
387                 else hNGammaCandidates[iCut] = new TH1I("GammaCandidates","GammaCandidates",50,0,50);
388                 fESDList[iCut]->Add(hNGammaCandidates[iCut]);
389                 if(fIsHeavyIon == 1) hNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",4000,0,4000,100,0,100);
390                 else if(fIsHeavyIon == 2) hNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",400,0,400,50,0,50);
391                 else hNGoodESDTracksVsNGammaCanditates[iCut] = new TH2F("GoodESDTracksVsGammaCandidates","GoodESDTracksVsGammaCandidates",200,0,200,50,0,50);
392                 fESDList[iCut]->Add(hNGoodESDTracksVsNGammaCanditates[iCut]);
393                 
394                 
395                 if(fIsHeavyIon == 1) hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",30000,0,30000);
396                 else if(fIsHeavyIon == 2) hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",2500,0,2500);
397                 else hNV0Tracks[iCut] = new TH1I("V0 Multiplicity","V0 Multiplicity",1500,0,1500);
398                 fESDList[iCut]->Add(hNV0Tracks[iCut]);
399                 hEtaShift[iCut] = new TProfile("Eta Shift","Eta Shift",1, -0.5,0.5);
400                 fESDList[iCut]->Add(hEtaShift[iCut]);
401                 hESDConvGammaPt[iCut] = new TH1F("ESD_ConvGamma_Pt","ESD_ConvGamma_Pt",250,0,25);
402                 fESDList[iCut]->Add(hESDConvGammaPt[iCut]);
403                 
404                 if (fDoPhotonQA == 2){
405                         fPhotonDCAList[iCut] = new TList();
406                         fPhotonDCAList[iCut]->SetName(Form("%s_%s_%s Photon DCA tree",cutstringEvent.Data(),cutstringPhoton.Data(),cutstringMeson.Data()));
407                         fPhotonDCAList[iCut]->SetOwner(kTRUE);
408                         fCutFolder[iCut]->Add(fPhotonDCAList[iCut]);
409                         
410 //            tESDConvGammaPtDcazCat[iCut] = new TTree("ESD_ConvGamma_Pt_Dcaz_R_Eta","ESD_ConvGamma_Pt_Dcaz_R_Eta_Cat");
411 //            tESDConvGammaPtDcazCat[iCut]->Branch("Pt",&fPtGamma,"fPtGamma/F");
412 //            tESDConvGammaPtDcazCat[iCut]->Branch("DcaZPhoton",&fDCAzPhoton,"fDCAzPhoton/F");
413                         //          tESDConvGammaPtDcazCat[iCut]->Branch("R",&fRConvPhoton,"fRConvPhoton/F");
414                         //          tESDConvGammaPtDcazCat[iCut]->Branch("Eta",&fEtaPhoton,"fEtaPhoton/F");
415                         
416                 //    tESDConvGammaPtDcazCat[iCut]->Branch("cat",&iCatPhoton,"iCatPhoton/b");
417
418         //        fPhotonDCAList[iCut]->Add(tESDConvGammaPtDcazCat[iCut]);
419                 }
420                 
421                 if (fDoPhotonQA > 0){
422                         hESDConvGammaR[iCut] = new TH1F("ESD_ConvGamma_R","ESD_ConvGamma_R",800,0,200);
423                         fESDList[iCut]->Add(hESDConvGammaR[iCut]);
424                         hESDConvGammaEta[iCut] = new TH1F("ESD_ConvGamma_Eta","ESD_ConvGamma_Eta",2000,-2,2);
425                         fESDList[iCut]->Add(hESDConvGammaEta[iCut]);
426                 }
427                 
428                 if(fDoMesonAnalysis){
429                         hESDMotherInvMassPt[iCut] = new TH2F("ESD_Mother_InvMass_Pt","ESD_Mother_InvMass_Pt",800,0,0.8,250,0,25);
430                         fESDList[iCut]->Add(hESDMotherInvMassPt[iCut]);
431                         hESDMotherBackInvMassPt[iCut] = new TH2F("ESD_Background_InvMass_Pt","ESD_Background_InvMass_Pt",800,0,0.8,250,0,25);
432                         fESDList[iCut]->Add(hESDMotherBackInvMassPt[iCut]);
433                         hESDMotherInvMassEalpha[iCut] = new TH2F("ESD_Mother_InvMass_vs_E_alpha","ESD_Mother_InvMass_vs_E_alpha",800,0,0.8,250,0,25);
434                         fESDList[iCut]->Add(hESDMotherInvMassEalpha[iCut]);
435                         if (fDoMesonQA == 2){
436                                 fMesonDCAList[iCut] = new TList();
437                                 fMesonDCAList[iCut]->SetName(Form("%s_%s_%s Meson DCA tree",cutstringEvent.Data() ,cutstringPhoton.Data(),cutstringMeson.Data()));
438                                 fMesonDCAList[iCut]->SetOwner(kTRUE);
439                                 fCutFolder[iCut]->Add(fMesonDCAList[iCut]);
440                                 
441                                 tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut] = new TTree("ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag","ESD_Mesons_InvMass_Pt_DcazMin_DcazMax_Flag");
442                                 tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("InvMass",&fInvMass,"fInvMass/F");
443                                 tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("Pt",&fPt,"fPt/F");
444                                 tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMin",&fDCAzGammaMin,"fDCAzGammaMin/F");
445                                 tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]->Branch("DcaZMax",&fDCAzGammaMax,"fDCAzGammaMax/F");
446                                 
447                                 fMesonDCAList[iCut]->Add(tESDMesonsInvMassPtDcazMinDcazMaxFlag[iCut]);
448                                 
449                         }
450                         if (fDoMesonQA > 0 ){
451                                 hESDMotherPi0PtY[iCut] = new TH2F("ESD_MotherPi0_Pt_Y","ESD_MotherPi0_Pt_Y",150,0.03,15.,150,-1.5,1.5);
452                                 SetLogBinningXTH2(hESDMotherPi0PtY[iCut]);
453                                 fESDList[iCut]->Add(hESDMotherPi0PtY[iCut]);
454                                 hESDMotherEtaPtY[iCut] = new TH2F("ESD_MotherEta_Pt_Y","ESD_MotherEta_Pt_Y",150,0.03,15.,150,-1.5,1.5);
455                                 SetLogBinningXTH2(hESDMotherEtaPtY[iCut]);
456                                 fESDList[iCut]->Add(hESDMotherEtaPtY[iCut]);
457                                 hESDMotherPi0PtAlpha[iCut] = new TH2F("ESD_MotherPi0_Pt_Alpha","ESD_MotherPi0_Pt_Alpha",150,0.03,15.,100,0,1);
458                                 SetLogBinningXTH2(hESDMotherPi0PtAlpha[iCut]);
459                                 fESDList[iCut]->Add(hESDMotherPi0PtAlpha[iCut]);
460                                 hESDMotherEtaPtAlpha[iCut] = new TH2F("ESD_MotherEta_Pt_Alpha","ESD_MotherEta_Pt_Alpha",150,0.03,15.,100,0,1);
461                                 SetLogBinningXTH2(hESDMotherEtaPtAlpha[iCut]);
462                                 fESDList[iCut]->Add(hESDMotherEtaPtAlpha[iCut]);
463                                 hESDMotherPi0PtOpenAngle[iCut] = new TH2F("ESD_MotherPi0_Pt_OpenAngle","ESD_MotherPi0_Pt_OpenAngle",150,0.03,15.,100,0,TMath::Pi());
464                                 SetLogBinningXTH2(hESDMotherPi0PtOpenAngle[iCut]);
465                                 fESDList[iCut]->Add(hESDMotherPi0PtOpenAngle[iCut]);
466                                 hESDMotherEtaPtOpenAngle[iCut] = new TH2F("ESD_MotherEta_Pt_OpenAngle","ESD_MotherEta_Pt_OpenAngle",150,0.03,15.,100,0,TMath::Pi());
467                                 SetLogBinningXTH2(hESDMotherEtaPtOpenAngle[iCut]);
468                                 fESDList[iCut]->Add(hESDMotherEtaPtOpenAngle[iCut]);
469                         }
470                         
471                         
472                 }
473                 
474                 
475         }
476 //     if(fDoMesonAnalysis){
477 //         InitBack(); // Init Background Handler
478 //     }
479         
480
481         
482         fV0Reader=(AliV0ReaderV1*)AliAnalysisManager::GetAnalysisManager()->GetTask("V0ReaderV1");
483         if(!fV0Reader){printf("Error: No V0 Reader");return;} // GetV0Reader
484         
485         if(fV0Reader)
486                 if((AliConvEventCuts*)fV0Reader->GetEventCuts())
487                         if(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms())
488                                 fOutputContainer->Add(((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetCutHistograms());
489         
490         if(fV0Reader)
491                 if((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())
492                         if(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms())
493                                 fOutputContainer->Add(((AliConversionPhotonCuts*)fV0Reader->GetConversionCuts())->GetCutHistograms());
494         
495         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
496                 if(!((AliConvEventCuts*)fEventCutArray->At(iCut))) continue;
497                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms()){
498                         fCutFolder[iCut]->Add(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutHistograms());
499                 }
500                 if(!((AliConversionPhotonCuts*)fCutArray->At(iCut))) continue;
501                 if(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms()){
502                         fCutFolder[iCut]->Add(((AliConversionPhotonCuts*)fCutArray->At(iCut))->GetCutHistograms());
503                 }
504                 if(fDoMesonAnalysis){
505                         if(!((AliConversionMesonCuts*)fMesonCutArray->At(iCut))) continue;
506                         if(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms()){
507                                 fCutFolder[iCut]->Add(((AliConversionMesonCuts*)fMesonCutArray->At(iCut))->GetCutHistograms());
508                         }
509                 }
510         }
511         
512         fhistoEPVZ = new TH1D("EPVZ", "EPVZ", 60, -TMath::Pi()/2, TMath::Pi()/2);
513     fOutputContainer->Add(fhistoEPVZ);
514
515
516         
517         PostData(1, fOutputContainer);
518
519         fFlowEvent = new AliFlowEvent(10000);
520         PostData(2, fFlowEvent);
521         
522 }
523
524 //_____________________________________________________________________________
525 Bool_t AliAnalysisTaskGammaConvFlow::Notify()
526 {
527         for(Int_t iCut = 0; iCut<fnCuts;iCut++){
528                 if(!((AliConvEventCuts*)fEventCutArray->At(iCut))->GetDoEtaShift()){
529                         hEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
530                         continue; // No Eta Shift requested, continue
531                 }
532                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift() == 0.0){ // Eta Shift requested but not set, get shift automatically
533                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCorrectEtaShiftFromPeriod(fV0Reader->GetPeriodName());
534                         hEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
535                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
536                         continue;
537                 }
538                 else{
539                         printf(" Gamma Conversion Task %s :: Eta Shift Manually Set to %f \n\n",
540                                 (((AliConvEventCuts*)fEventCutArray->At(iCut))->GetCutNumber()).Data(),((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift());
541                         hEtaShift[iCut]->Fill(0.,(((AliConvEventCuts*)fEventCutArray->At(iCut))->GetEtaShift()));
542                         ((AliConvEventCuts*)fEventCutArray->At(iCut))->DoEtaShift(kFALSE); // Eta Shift Set, make sure that it is called only once
543                 }
544         }
545         return kTRUE;
546 }
547 //_____________________________________________________________________________
548 void AliAnalysisTaskGammaConvFlow::UserExec(Option_t *)
549 {
550         //
551         // Called for each event
552         //
553         Int_t eventQuality = ((AliConvEventCuts*)fV0Reader->GetEventCuts())->GetEventQuality();
554         if(eventQuality == 2 || eventQuality == 3){// Event Not Accepted due to MC event missing or wrong trigger for V0ReaderV1
555                 for(Int_t iCut = 0; iCut<fnCuts; iCut++){
556                         hNEvents[iCut]->Fill(eventQuality);
557                 }
558                 return;
559         }
560         
561         fInputEvent = InputEvent();
562
563         fReaderGammas = fV0Reader->GetReconstructedGammas(); // Gammas from default Cut
564         
565         // ------------------- BeginEvent ----------------------------
566         
567         AliEventplane *EventPlane = fInputEvent->GetEventplane();
568         if(fIsHeavyIon ==1)fEventPlaneAngle = EventPlane->GetEventplane("V0",fInputEvent,2);
569         else fEventPlaneAngle=0.0;
570         
571         SetNullCuts(fInputEvent);
572         PrepareFlowEvent(fInputEvent->GetNumberOfTracks(),fFlowEvent);    //Calculate event plane Qvector and EP resolution for inclusive
573
574         
575
576         
577 //      for(Int_t iCut = 0; iCut<fnCuts; iCut++){
578                 Int_t iCut = 0;
579                 fiCut = iCut;
580                 Int_t eventNotAccepted =
581                 ((AliConvEventCuts*)fEventCutArray->At(iCut))
582                 ->IsEventAcceptedByCut(fV0Reader->GetEventCuts(),fInputEvent,fMCEvent,fIsHeavyIon,kFALSE);
583                 if(eventNotAccepted){
584                         // cout << "event rejected due to wrong trigger: " <<eventNotAccepted << endl;
585                         hNEvents[iCut]->Fill(eventNotAccepted); // Check Centrality, PileUp, SDD and V0AND --> Not Accepted => eventQuality = 1
586                         return;
587                 }
588                 
589                 if(eventQuality != 0){// Event Not Accepted
590                         // cout << "event rejected due to: " <<eventQuality << endl;
591                         hNEvents[iCut]->Fill(eventQuality);
592                         return;
593                 }
594                                 
595                 hNEvents[iCut]->Fill(eventQuality); // Should be 0 here
596                 hNGoodESDTracks[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks());
597                 if(((AliConvEventCuts*)fEventCutArray->At(iCut))->IsHeavyIon() == 2)    hNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A());
598                 else hNV0Tracks[iCut]->Fill(fInputEvent->GetVZEROData()->GetMTotV0A()+fInputEvent->GetVZEROData()->GetMTotV0C());
599                 
600                 ProcessPhotonCandidates(); // Process this cuts gammas
601                 ProcessPhotonCandidatesforV2(); // Process this cuts gammas and do v2
602
603                 hNGammaCandidates[iCut]->Fill(fGammaCandidates->GetEntries());
604                 hNGoodESDTracksVsNGammaCanditates[iCut]->Fill(fV0Reader->GetNumberOfPrimaryTracks(),fGammaCandidates->GetEntries());
605
606                 fGammaCandidates->Clear(); // delete this cuts good gammas
607 //      }
608         
609         fhistoEPVZ->Fill(fEventPlaneAngle);
610
611         PostData(1, fOutputContainer);
612         PostData(2, fFlowEvent);
613
614 }
615 //________________________________________________________________________
616 void AliAnalysisTaskGammaConvFlow::ProcessPhotonCandidates()
617 {
618     Int_t nV0 = 0;
619     TList *GammaCandidatesStepOne = new TList();
620     TList *GammaCandidatesStepTwo = new TList();
621     // Loop over Photon Candidates allocated by ReaderV1
622     for(Int_t i = 0; i < fReaderGammas->GetEntriesFast(); i++){
623         AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) fReaderGammas->At(i);
624         if(!PhotonCandidate) continue;
625         fIsFromMBHeader = kTRUE;
626
627         
628         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->PhotonIsSelected(PhotonCandidate,fInputEvent)) continue;
629         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->InPlaneOutOfPlaneCut(PhotonCandidate->GetPhotonPhi(),fEventPlaneAngle)) continue;
630         if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
631            !((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
632             fGammaCandidates->Add(PhotonCandidate); // if no second loop is required add to events good gammas
633             
634             if(fIsFromMBHeader){
635                 hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
636                 if (fDoPhotonQA > 0){
637                     hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
638                     hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
639                 }
640             }
641
642             
643             if (fIsFromMBHeader && fDoPhotonQA == 2){
644                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
645                     fPtGamma = PhotonCandidate->Pt();
646                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
647                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
648                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
649                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
650                    // tESDConvGammaPtDcazCat[fiCut]->Fill();
651                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
652                     fPtGamma = PhotonCandidate->Pt();
653                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
654                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
655                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
656                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
657                   //  tESDConvGammaPtDcazCat[fiCut]->Fill();
658                 }
659             }
660         } else if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){ // if Shared Electron cut is enabled, Fill array, add to step one
661             ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->FillElectonLabelArray(PhotonCandidate,nV0);
662             nV0++;
663             GammaCandidatesStepOne->Add(PhotonCandidate);
664         } else if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut() &&
665                   ((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // shared electron is disabled, step one not needed -> step two
666             GammaCandidatesStepTwo->Add(PhotonCandidate);
667         }
668     }
669     
670     
671     
672     if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseElecSharingCut()){
673         for(Int_t i = 0;i<GammaCandidatesStepOne->GetEntries();i++){
674             AliAODConversionPhoton *PhotonCandidate= (AliAODConversionPhoton*) GammaCandidatesStepOne->At(i);
675             if(!PhotonCandidate) continue;
676             fIsFromMBHeader = kTRUE;
677
678             if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectSharedElectronV0s(PhotonCandidate,i,GammaCandidatesStepOne->GetEntries())) continue;
679             if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){ // To Colse v0s cut diabled, step two not needed
680                 fGammaCandidates->Add(PhotonCandidate);
681                 if(fIsFromMBHeader){
682                     hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
683                     if (fDoPhotonQA > 0){
684                         hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
685                         hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
686                     }
687                 }
688             }
689
690                 
691                 GammaCandidatesStepTwo->Add(PhotonCandidate); // Close v0s cut enabled -> add to list two
692             
693             if (fIsFromMBHeader && fDoPhotonQA == 2){
694                 if (fIsHeavyIon ==1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
695                     fPtGamma = PhotonCandidate->Pt();
696                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
697                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
698                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
699                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
700                 //    tESDConvGammaPtDcazCat[fiCut]->Fill();
701                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
702                     fPtGamma = PhotonCandidate->Pt();
703                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
704                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
705                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
706                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
707                 //    tESDConvGammaPtDcazCat[fiCut]->Fill();
708                 }
709             }
710         }
711     }
712     if(((AliConversionPhotonCuts*)fCutArray->At(fiCut))->UseToCloseV0sCut()){
713         for(Int_t i = 0;i<GammaCandidatesStepTwo->GetEntries();i++){
714             AliAODConversionPhoton* PhotonCandidate = (AliAODConversionPhoton*) GammaCandidatesStepTwo->At(i);
715             if(!PhotonCandidate) continue;
716             fIsFromMBHeader = kTRUE;
717
718             if(!((AliConversionPhotonCuts*)fCutArray->At(fiCut))->RejectToCloseV0s(PhotonCandidate,GammaCandidatesStepTwo,i)) continue;
719             fGammaCandidates->Add(PhotonCandidate); // Add gamma to current cut TList
720             if(fIsFromMBHeader){
721                 hESDConvGammaPt[fiCut]->Fill(PhotonCandidate->Pt());
722                 if (fDoPhotonQA > 0){
723                     hESDConvGammaR[fiCut]->Fill(PhotonCandidate->GetConversionRadius());
724                     hESDConvGammaEta[fiCut]->Fill(PhotonCandidate->Eta());
725                 }
726             }
727
728             if (fIsFromMBHeader){
729                 if (fIsHeavyIon == 1 && PhotonCandidate->Pt() > 0.399 && PhotonCandidate->Pt() < 12.){
730                     fPtGamma = PhotonCandidate->Pt();
731                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
732                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
733                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
734                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
735                //     tESDConvGammaPtDcazCat[fiCut]->Fill();
736                 } else if ( PhotonCandidate->Pt() > 0.299 && PhotonCandidate->Pt() < 16.){
737                     fPtGamma = PhotonCandidate->Pt();
738                     fDCAzPhoton = PhotonCandidate->GetDCAzToPrimVtx();
739                     fRConvPhoton = PhotonCandidate->GetConversionRadius();
740                     fEtaPhoton = PhotonCandidate->GetPhotonEta();
741                     iCatPhoton = PhotonCandidate->GetPhotonQuality();
742                   //  tESDConvGammaPtDcazCat[fiCut]->Fill();
743                 }
744             }
745         }
746     }
747     
748     delete GammaCandidatesStepOne;
749     GammaCandidatesStepOne = 0x0;
750     delete GammaCandidatesStepTwo;
751     GammaCandidatesStepTwo = 0x0;
752     
753 }
754 //________________________________________________________________________
755 void AliAnalysisTaskGammaConvFlow::UpdateEventByEventData(){
756     //see header file for documentation
757     if(fGammaCandidates->GetEntries() >0 ){
758         if(((AliConversionMesonCuts*)fMesonCutArray->At(fiCut))->UseTrackMultiplicity()){
759             fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fV0Reader->GetNumberOfPrimaryTracks(),fEventPlaneAngle);
760         }
761         else{ // means we use #V0s for multiplicity
762             fBGHandler[fiCut]->AddEvent(fGammaCandidates,fInputEvent->GetPrimaryVertex()->GetX(),fInputEvent->GetPrimaryVertex()->GetY(),fInputEvent->GetPrimaryVertex()->GetZ(),fGammaCandidates->GetEntries(),fEventPlaneAngle);
763         }
764     }
765 }
766 //________________________________________________________________________
767 void AliAnalysisTaskGammaConvFlow::SetLogBinningXTH2(TH2* histoRebin){
768     TAxis *axisafter = histoRebin->GetXaxis(); 
769     Int_t bins = axisafter->GetNbins();
770     Double_t from = axisafter->GetXmin();
771     Double_t to = axisafter->GetXmax();
772     Double_t *newbins = new Double_t[bins+1];
773     newbins[0] = from;
774     Double_t factor = TMath::Power(to/from, 1./bins);
775     for(Int_t i=1; i<=bins; ++i) newbins[i] = factor * newbins[i-1];
776     axisafter->Set(bins, newbins);
777     delete [] newbins;
778 }
779
780 //________________________________________________________________________
781 void AliAnalysisTaskGammaConvFlow::Terminate(const Option_t *)
782 {
783     
784     //fOutputContainer->Print(); // Will crash on GRID
785 }
786
787 //________________________________________________________________________
788 void AliAnalysisTaskGammaConvFlow::ProcessPhotonCandidatesforV2()
789 {
790
791         // Loop over Photon Candidates allocated by ReaderV1    
792         for(Int_t i = 0; i < fGammaCandidates->GetEntries(); i++){
793                 
794                 AliAODConversionPhoton *gammaForv2=dynamic_cast<AliAODConversionPhoton*>(fGammaCandidates->At(i));
795                 if (gammaForv2 == NULL) return;
796                 AliFlowTrack *sTrack = new AliFlowTrack();
797                 sTrack->SetForRPSelection(kFALSE);
798                 sTrack->SetForPOISelection(kTRUE);
799                 sTrack->SetMass(263732);
800                 sTrack->SetPt(gammaForv2->Pt());
801                 sTrack->SetPhi(gammaForv2->GetPhotonPhi());
802                 sTrack->SetEta(gammaForv2->GetPhotonEta());
803
804         /*     for(int iRPs=0; iRPs!=fFlowEvent->NumberOfTracks(); ++iRPs)
805                 {
806                         //   cout << " no of rps " << iRPs << endl;
807                         AliFlowTrack *iRP = dynamic_cast<AliFlowTrack*>(fFlowEvent->GetTrack( iRPs ));
808                         if (!iRP) continue;
809                         if (!iRP->InRPSelection()) continue;
810                         if( sTrack->GetID() == iRP->GetID())
811                         {
812                                 if(fDebug) printf(" was in RP set");
813                                 //  cout << sTrack->GetID() <<"   ==  " << iRP->GetID() << " was in RP set====REMOVED" <<endl;
814                                 iRP->SetForRPSelection(kFALSE);
815                                 // fFlowEvent->SetNumberOfRPs(fFlowEvent->GetNumberOfRPs() - 1);
816                         }
817                 } //end of for loop on RPs*/
818                 fFlowEvent->InsertTrack(((AliFlowTrack*) sTrack));
819                 fFlowEvent->SetNumberOfPOIs(fFlowEvent->GetNumberOfPOIs()+1);
820         }
821 }
822
823 //_____________________________________________________________________________
824 template <typename T> void AliAnalysisTaskGammaConvFlow::SetNullCuts(T* event)
825 {
826         // Set null cuts
827         if (fDebug) cout << " fCutsRP " << fCutsRP << endl;
828         fCutsRP->SetEvent(event, MCEvent());
829         fNullCuts->SetParamType(AliFlowTrackCuts::kGlobal);
830         fNullCuts->SetPtRange(+1, -1); // select nothing QUICK
831         fNullCuts->SetEtaRange(+1, -1); // select nothing VZERO
832         fNullCuts->SetEvent(event, MCEvent());
833 }
834
835 //_____________________________________________________________________________
836 void AliAnalysisTaskGammaConvFlow::PrepareFlowEvent(Int_t iMulti, AliFlowEvent *FlowEv) const
837 {
838    //Prepare flow events
839     FlowEv->ClearFast();
840     FlowEv->Fill(fCutsRP, fNullCuts);
841     FlowEv->SetReferenceMultiplicity(iMulti);
842     FlowEv->DefineDeadZone(0, 0, 0, 0);
843     //  FlowEv->TagSubeventsInEta(-0.7, 0, 0, 0.7);
844 }
845
846
847
848
849
850