coding violations corrected
[u/mrichter/AliRoot.git] / PWG4 / PartCorr / AliAnalysisTaskGammaConversion.cxx
1 /**************************************************************************\r
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
3  *                                                                        *\r
4  * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt                        *\r
5  * Version 1.0                                                            *\r
6  *                                                                        *\r
7  * Permission to use, copy, modify and distribute this software and its   *\r
8  * documentation strictly for non-commercial purposes is hereby granted   *\r
9  * without fee, provided that the above copyright notice appears in all   *\r
10  * copies and that both the copyright notice and this permission notice   *\r
11  * appear in the supporting documentation. The authors make no claims     *\r
12  * about the suitability of this software for any purpose. It is          *\r
13  * provided "as is" without express or implied warranty.                  *\r
14  **************************************************************************/\r
15 \r
16 ////////////////////////////////////////////////\r
17 //--------------------------------------------- \r
18 // Class used to do analysis on conversion pairs\r
19 //---------------------------------------------\r
20 ////////////////////////////////////////////////\r
21 \r
22 // root\r
23 #include <TChain.h>\r
24 \r
25 // analysis\r
26 #include "AliAnalysisTaskGammaConversion.h"\r
27 #include "AliStack.h"\r
28 #include "AliLog.h"\r
29 #include <vector>\r
30 \r
31 class AliKFVertex;\r
32 class AliAODHandler;\r
33 class AliAODEvent;\r
34 class ALiESDEvent;\r
35 class AliMCEvent;\r
36 class AliMCEventHandler;\r
37 class AliESDInputHandler;\r
38 class AliAnalysisManager;\r
39 class Riostream;\r
40 class TFile;\r
41 class TInterpreter;\r
42 class TSystem;\r
43 class TROOT;\r
44 \r
45 ClassImp(AliAnalysisTaskGammaConversion)\r
46 \r
47 \r
48 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():\r
49   AliAnalysisTaskSE(),\r
50   fV0Reader(NULL),\r
51   fStack(NULL),\r
52   fOutputContainer(NULL),\r
53   fHistograms(NULL),\r
54   fDoMCTruth(kFALSE),\r
55   fMCAllGammas(),\r
56   fMCPi0s(),\r
57   fMCEtas(),\r
58   fMCGammaChic(),\r
59   fKFReconstructedGammas(),\r
60   fElectronMass(-1),\r
61   fGammaMass(-1),\r
62   fPi0Mass(-1),\r
63   fEtaMass(-1),\r
64   fGammaWidth(-1),\r
65   fPi0Width(-1),\r
66   fEtaWidth(-1),\r
67   fCalculateBackground(kFALSE)\r
68 {\r
69   // Default constructor\r
70   // Common I/O in slot 0\r
71   DefineInput (0, TChain::Class());\r
72   DefineOutput(0, TTree::Class());\r
73 \r
74   // Your private output\r
75   DefineOutput(1, TList::Class());\r
76 }\r
77 \r
78 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
79   AliAnalysisTaskSE(name),\r
80   fV0Reader(NULL),\r
81   fStack(NULL),\r
82   fOutputContainer(0x0),\r
83   fHistograms(NULL),\r
84   fDoMCTruth(kFALSE),\r
85   fMCAllGammas(),\r
86   fMCPi0s(),\r
87   fMCEtas(),\r
88   fMCGammaChic(),\r
89   fKFReconstructedGammas(),\r
90   fElectronMass(-1),\r
91   fGammaMass(-1),\r
92   fPi0Mass(-1),\r
93   fEtaMass(-1),\r
94   fGammaWidth(-1),\r
95   fPi0Width(-1),\r
96   fEtaWidth(-1),\r
97   fCalculateBackground(kFALSE)\r
98 {\r
99   // Common I/O in slot 0\r
100   DefineInput (0, TChain::Class());\r
101   DefineOutput(0, TTree::Class());\r
102   \r
103   // Your private output\r
104   DefineOutput(1, TList::Class());\r
105 }\r
106 \r
107 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
108 {\r
109   // Remove all pointers\r
110  \r
111   if(fOutputContainer){\r
112     fOutputContainer->Clear() ; \r
113     delete fOutputContainer ;\r
114   }\r
115   if(fHistograms){\r
116     delete fHistograms;\r
117   }\r
118   if(fV0Reader){\r
119     delete fV0Reader;\r
120   }\r
121 }\r
122 \r
123 \r
124 void AliAnalysisTaskGammaConversion::Init()\r
125 {\r
126   // Initialization\r
127 }\r
128 \r
129 \r
130 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
131 {\r
132   // Execute analysis for current event\r
133   \r
134   ConnectInputData("");\r
135   \r
136   //clear vectors\r
137   fMCAllGammas.clear();\r
138   fMCPi0s.clear();\r
139   fMCEtas.clear();\r
140   fMCGammaChic.clear();\r
141 \r
142   fKFReconstructedGammas.clear();\r
143 \r
144   //Clear the data in the v0Reader\r
145   fV0Reader->UpdateEventByEventData();\r
146 \r
147   // Process the MC information\r
148   if(fDoMCTruth){\r
149     ProcessMCData();\r
150   }\r
151 \r
152   // Process the v0 information\r
153   ProcessV0s();\r
154 \r
155   //calculate background if flag is set\r
156   if(fCalculateBackground){\r
157     CalculateBackground();\r
158   }\r
159 \r
160   // Process reconstructed gammas\r
161   ProcessGammasForNeutralMesonAnalysis();\r
162 \r
163   PostData(1, fOutputContainer);\r
164   \r
165 }\r
166 \r
167 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
168   // see header file for documentation\r
169 \r
170   if(fV0Reader == NULL){\r
171     // Write warning here cuts and so on are default if this ever happens\r
172   }\r
173   fV0Reader->Initialize();\r
174 }\r
175 \r
176 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
177   // see header file for documentation\r
178   \r
179   fStack = fV0Reader->GetMCStack();\r
180 \r
181   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
182     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
183 \r
184     if (!particle) {\r
185       //print warning here\r
186       continue;\r
187     }\r
188     \r
189     if(particle->Pt()<fV0Reader->GetPtCut()){\r
190       continue;\r
191     }\r
192     \r
193     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){\r
194       continue;\r
195     }\r
196 \r
197     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point\r
198       continue;\r
199     }\r
200 \r
201     Double_t tmpPhi=particle->Phi();\r
202     if(particle->Phi()> TMath::Pi()){\r
203       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
204     }\r
205 \r
206     \r
207     //process the gammas\r
208     if (particle->GetPdgCode()== 22){\r
209       fMCAllGammas.push_back(particle);\r
210       if(particle->GetMother(0)>-1){ //Means we have a mother\r
211         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.\r
212           fHistograms->FillHistogram("MC_Gamma_Energy", particle->Energy());\r
213           fHistograms->FillHistogram("MC_Gamma_Pt", particle->Pt());\r
214           fHistograms->FillHistogram("MC_Gamma_Eta", particle->Eta());\r
215           \r
216           fHistograms->FillHistogram("MC_Gamma_Phi", tmpPhi);\r
217 \r
218           //adding the conversion points from all gammas with e+e- daughters\r
219           if(particle->GetNDaughters() == 2){\r
220             TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
221             TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
222             \r
223             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){\r
224               continue;\r
225             }\r
226             \r
227             if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||\r
228                (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
229 \r
230               // begin Mapping \r
231               Int_t rBin    = fHistograms->GetRBin(daughter0->R());\r
232               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());\r
233               \r
234               TString nameMCMappingPhiR="";\r
235               nameMCMappingPhiR.Form("MC_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
236               fHistograms->FillHistogram(nameMCMappingPhiR, daughter0->Vz(), particle->Eta());\r
237               \r
238               TString nameMCMappingPhi="";\r
239               nameMCMappingPhi.Form("MC_EP_Mapping-Phi%02d",phiBin);\r
240               fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
241               \r
242               TString nameMCMappingR="";\r
243               nameMCMappingR.Form("MC_EP_Mapping-R%02d",rBin);\r
244               fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
245               \r
246               TString nameMCMappingPhiInR="";\r
247               nameMCMappingPhiInR.Form("MC_EP_Mapping_Phi_vs_R_R-%02d",rBin);\r
248               fHistograms->FillHistogram(nameMCMappingPhiInR, daughter0->R(), tmpPhi);\r
249               //end mapping\r
250 \r
251               fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
252               fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
253               fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
254               fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
255             }\r
256           }\r
257         }\r
258         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
259             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
260             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
261         ){ \r
262           fMCGammaChic.push_back(particle);\r
263          }\r
264       }\r
265       else{//means we have a primary particle\r
266         fHistograms->FillHistogram("MC_DirectGamma_Energy",particle->Energy());\r
267         fHistograms->FillHistogram("MC_DirectGamma_Pt", particle->Pt());\r
268         fHistograms->FillHistogram("MC_DirectGamma_Eta", particle->Eta());\r
269         fHistograms->FillHistogram("MCDirectGammaPhi", tmpPhi);\r
270 \r
271         //adding the conversion points from all gammas with e+e- daughters\r
272         if(particle->GetNDaughters() == 2){\r
273           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
274           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
275           if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11) ||\r
276              (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
277             \r
278             fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
279             fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
280             fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
281             fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
282 \r
283           }\r
284         }\r
285 \r
286       }\r
287     }\r
288     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron\r
289       if(particle->GetMother(0)>-1){ // means we have a mother\r
290         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother\r
291           if(particle->GetPdgCode() == 11){//electron \r
292             fHistograms->FillHistogram("MC_E_Energy", particle->Energy());\r
293             fHistograms->FillHistogram("MC_E_Pt", particle->Pt());\r
294             fHistograms->FillHistogram("MC_E_Eta", particle->Eta());\r
295             fHistograms->FillHistogram("MC_E_Phi", tmpPhi);\r
296           }\r
297           if(particle->GetPdgCode() == -11){//positron \r
298             fHistograms->FillHistogram("MC_P_Energy", particle->Energy());\r
299             fHistograms->FillHistogram("MC_P_Pt", particle->Pt());\r
300             fHistograms->FillHistogram("MC_P_Eta", particle->Eta());\r
301             fHistograms->FillHistogram("MC_P_Phi", tmpPhi);\r
302           }\r
303         }\r
304       }\r
305     }\r
306     else if(particle->GetNDaughters() == 2){\r
307 \r
308       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
309       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
310       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters\r
311         \r
312         if(particle->GetPdgCode()==111){//Pi0\r
313           \r
314           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
315 \r
316           if( iTracks >= fStack->GetNprimary()){\r
317             \r
318             fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
319 \r
320             fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
321             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
322             fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
323             fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
324             fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
325             fHistograms->FillHistogram("MC_Pi0_Secondaries_OpeningAngle_Gamma", GetMCOpeningAngle(daughter0,daughter1));\r
326             fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
327           }\r
328           else{\r
329             fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());\r
330 \r
331             fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
332             fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
333             fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
334             fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
335             fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
336             fHistograms->FillHistogram("MC_Pi0_OpeningAngle_Gamma", GetMCOpeningAngle(daughter0,daughter1));\r
337             fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
338           }\r
339         }\r
340         else if(particle->GetPdgCode()==221){//Eta\r
341           fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
342 \r
343           fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
344           fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
345           fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
346           fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
347           fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
348           fHistograms->FillHistogram("MC_Eta_OpeningAngle_Gamma", GetMCOpeningAngle(daughter0,daughter1));\r
349           fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
350         }\r
351         \r
352         //the match data should be filled no matter which mother the gamma-gamma comes from\r
353         fHistograms->FillHistogram("MC_Match_Gamma_R", particle->R());\r
354         fHistograms->FillHistogram("MC_Match_Gamma_ZR", particle->Vz(),particle->R());\r
355         fHistograms->FillHistogram("MC_Match_Gamma_XY", particle->Vx(),particle->Vy());\r
356         fHistograms->FillHistogram("MC_Match_Gamma_Mass", particle->GetCalcMass());\r
357         fHistograms->FillHistogram("MC_Match_Gamma_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
358         fHistograms->FillHistogram("MC_Match_Gamma_Energy", particle->Energy());\r
359         fHistograms->FillHistogram("MC_Match_Gamma_Pt", particle->Pt());\r
360         fHistograms->FillHistogram("MC_Match_Gamma_Eta", particle->Eta());\r
361         fHistograms->FillHistogram("MC_Match_Gamma_Phi",tmpPhi);\r
362       }\r
363     }\r
364   }\r
365 }\r
366 \r
367 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
368   // see header file for documentation\r
369 \r
370   Int_t nSurvivingV0s=0;\r
371   while(fV0Reader->NextV0()){\r
372     nSurvivingV0s++;\r
373     //-------------------------- filling v0 information -------------------------------------\r
374     fHistograms->FillHistogram("ESD_EP_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
375     fHistograms->FillHistogram("ESD_EP_R", fV0Reader->GetXYRadius());\r
376     fHistograms->FillHistogram("ESD_EP_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
377     fHistograms->FillHistogram("ESD_EP_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
378     \r
379     \r
380     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
381     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
382     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
383     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
384     \r
385     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
386     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
387     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
388     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
389     \r
390     fHistograms->FillHistogram("ESD_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
391     fHistograms->FillHistogram("ESD_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
392     fHistograms->FillHistogram("ESD_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
393     fHistograms->FillHistogram("ESD_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
394 \r
395 \r
396     // begin mapping\r
397     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
398     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
399     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
400     \r
401     TString nameESDMappingPhiR="";\r
402     nameESDMappingPhiR.Form("ESD_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
403     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
404 \r
405     TString nameESDMappingPhi="";\r
406     nameESDMappingPhi.Form("ESD_EP_Mapping-Phi%02d",phiBin);\r
407     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
408 \r
409     TString nameESDMappingR="";\r
410     nameESDMappingR.Form("ESD_EP_Mapping-R%02d",rBin);\r
411     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
412 \r
413     TString nameESDMappingPhiInR="";\r
414     nameESDMappingPhiInR.Form("ESD_EP_Mapping_Phi_vs_R_R-%02d",rBin);\r
415     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetXYRadius(), fV0Reader->GetMotherCandidatePhi());\r
416     // end mapping\r
417     \r
418     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
419 \r
420     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
421     if(fDoMCTruth){\r
422       if(fV0Reader->HasSameMCMother() == kFALSE){\r
423         continue;\r
424       }\r
425 \r
426       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
427       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
428       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){\r
429         continue;\r
430       }\r
431 \r
432       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
433         fHistograms->FillHistogram("ESD_Match_Gamma_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
434         fHistograms->FillHistogram("ESD_Match_Gamma_OpeningAngle", fV0Reader->GetOpeningAngle());\r
435         fHistograms->FillHistogram("ESD_Match_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
436         fHistograms->FillHistogram("ESD_Match_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
437         fHistograms->FillHistogram("ESD_Match_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
438 \r
439         fHistograms->FillHistogram("ESD_Match_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
440         fHistograms->FillHistogram("ESD_Match_Gamma_Mass", fV0Reader->GetMotherCandidateMass());\r
441         fHistograms->FillHistogram("ESD_Match_Gamma_Width", fV0Reader->GetMotherCandidateWidth());\r
442         fHistograms->FillHistogram("ESD_Match_Gamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
443         fHistograms->FillHistogram("ESD_Match_Gamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
444         fHistograms->FillHistogram("ESD_Match_Gamma_R", fV0Reader->GetXYRadius());\r
445         fHistograms->FillHistogram("ESD_Match_Gamma_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
446 \r
447         //resolution\r
448         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
449         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
450         Double_t resdPt = 0;\r
451         if(mcpt != 0){\r
452           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
453         }\r
454 \r
455         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
456         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
457         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
458         \r
459         Double_t resdZ = 0;\r
460         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
461           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
462         }\r
463         \r
464         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
465         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
466         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
467         \r
468         Double_t resdR = 0;\r
469         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
470           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
471         }\r
472         fHistograms->FillHistogram("Resolutiond_R", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
473         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
474         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
475         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
476       }\r
477     }\r
478   }\r
479   fHistograms->FillHistogram("NumberOfSurvivingV0s", nSurvivingV0s);\r
480   fHistograms->FillHistogram("NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
481 }\r
482 \r
483 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
484   // see header file for documentation\r
485 \r
486   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
487     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
488       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
489       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
490 \r
491       \r
492       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
493 \r
494       Double_t massTwoGammaCandidate =0.;\r
495       Double_t widthTwoGammaCandidate = 0.;\r
496       Double_t chi2TwoGammaCandidate =10000.;   \r
497       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
498       if(twoGammaCandidate->GetNDF()>0){\r
499         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
500         //      if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fChi2Cut){//TODO  find this out\r
501         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<10000){//TODO  find this out see line above\r
502 \r
503           TVector3 vectorTwoGammaCandidate(twoGammaCandidate->Px(),twoGammaCandidate->Py(),twoGammaCandidate->Pz());\r
504 \r
505           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);\r
506 \r
507           //Calculating by hand the radius\r
508           Double_t tmpX= twoGammaCandidate->GetX();\r
509           Double_t tmpY= twoGammaCandidate->GetY();\r
510           \r
511           Double_t radiusTwoGammaCandidate = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
512 \r
513           fHistograms->FillHistogram("ESD_TwoGammaCombination_OpeningAngleGamma", openingAngleTwoGammaCandidate);\r
514           fHistograms->FillHistogram("ESD_TwoGammaCombination_Energy", twoGammaCandidate->GetE());\r
515           fHistograms->FillHistogram("ESD_TwoGammaCombination_Pt", sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
516           fHistograms->FillHistogram("ESD_TwoGammaCombination_Eta", vectorTwoGammaCandidate.Eta());\r
517           fHistograms->FillHistogram("ESD_TwoGammaCombination_Phi", vectorTwoGammaCandidate.Phi());\r
518           fHistograms->FillHistogram("ESD_TwoGammaCombination_Mass", massTwoGammaCandidate);\r
519           fHistograms->FillHistogram("ESD_TwoGammaCombination_R", radiusTwoGammaCandidate);\r
520           fHistograms->FillHistogram("ESD_TwoGammaCombination_ZR", tmpY, radiusTwoGammaCandidate);\r
521           fHistograms->FillHistogram("ESD_TwoGammaCombination_XY", tmpX, tmpY);   \r
522           fHistograms->FillHistogram("InvMass_vs_Pt__Spectra",massTwoGammaCandidate ,sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
523         }\r
524       }\r
525       delete twoGammaCandidate;\r
526     }\r
527   }\r
528 \r
529 }\r
530 \r
531 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
532   // see header file for documentation\r
533 \r
534   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
535   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
536   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
537     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
538     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
539       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
540 \r
541       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
542 \r
543       Double_t massBG =0.;\r
544       Double_t widthBG = 0.;\r
545       Double_t chi2BG =10000.;  \r
546       backgroundCandidate->GetMass(massBG,widthBG);\r
547       if(backgroundCandidate->GetNDF()>0){\r
548         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
549         //      if(chi2Pi0>0 && chi2Pi0<fChi2Cut){//TODO  find this out\r
550         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2Cut()){//TODO  find this out see line above\r
551 \r
552           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());\r
553 \r
554           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
555 \r
556           //Calculating by hand the radius (find a better way)\r
557           Double_t tmpX= backgroundCandidate->GetX();\r
558           Double_t tmpY= backgroundCandidate->GetY();\r
559           \r
560           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
561 \r
562           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
563           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
564           fHistograms->FillHistogram("ESD_Background_Pt", sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
565           fHistograms->FillHistogram("ESD_Background_Eta", vectorBGCandidate.Eta());\r
566           fHistograms->FillHistogram("ESD_Background_Phi", vectorBGCandidate.Phi());\r
567           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
568           fHistograms->FillHistogram("ESD_Background_R", radiusBG);\r
569           fHistograms->FillHistogram("ESD_Background_ZR", tmpY, radiusBG);\r
570           fHistograms->FillHistogram("ESD_Background_XY", tmpX, tmpY);\r
571         }\r
572       }\r
573       delete backgroundCandidate;   \r
574     }\r
575   }\r
576 }\r
577 \r
578 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
579 {\r
580   // Terminate analysis\r
581   //\r
582   AliDebug(1,"Do nothing in Terminate");\r
583 }\r
584 \r
585 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
586 {\r
587   // Create the output container\r
588   if(fOutputContainer != NULL){\r
589     delete fOutputContainer;\r
590     fOutputContainer = NULL;\r
591   }\r
592   if(fOutputContainer == NULL){\r
593     fOutputContainer = new TList();\r
594   }\r
595   fHistograms->GetOutputContainer(fOutputContainer);\r
596   fOutputContainer->SetName(GetName()) ;  \r
597 }\r
598 \r
599 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
600   //helper function\r
601   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
602   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
603   return v3D0.Angle(v3D1);\r
604 }\r
605 \r