bugs 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   AliLog::SetGlobalLogLevel(AliLog::kError);\r
128 }\r
129 \r
130 \r
131 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
132 {\r
133   // Execute analysis for current event\r
134   \r
135   ConnectInputData("");\r
136   \r
137   //clear vectors\r
138   fMCAllGammas.clear();\r
139   fMCPi0s.clear();\r
140   fMCEtas.clear();\r
141   fMCGammaChic.clear();\r
142 \r
143   fKFReconstructedGammas.clear();\r
144 \r
145   //Clear the data in the v0Reader\r
146   fV0Reader->UpdateEventByEventData();\r
147 \r
148   // Process the MC information\r
149   if(fDoMCTruth){\r
150     ProcessMCData();\r
151   }\r
152 \r
153   // Process the v0 information\r
154   ProcessV0s();\r
155 \r
156   //calculate background if flag is set\r
157   if(fCalculateBackground){\r
158     CalculateBackground();\r
159   }\r
160 \r
161   // Process reconstructed gammas\r
162   ProcessGammasForNeutralMesonAnalysis();\r
163 \r
164   PostData(1, fOutputContainer);\r
165   \r
166 }\r
167 \r
168 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
169   // see header file for documentation\r
170 \r
171   if(fV0Reader == NULL){\r
172     // Write warning here cuts and so on are default if this ever happens\r
173   }\r
174   fV0Reader->Initialize();\r
175 }\r
176 \r
177 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
178   // see header file for documentation\r
179   \r
180   fStack = fV0Reader->GetMCStack();\r
181 \r
182   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
183     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
184 \r
185     if (!particle) {\r
186       //print warning here\r
187       continue;\r
188     }\r
189     \r
190     if(particle->Pt()<fV0Reader->GetPtCut()){\r
191       continue;\r
192     }\r
193     \r
194     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){\r
195       continue;\r
196     }\r
197 \r
198     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point\r
199       continue;\r
200     }\r
201 \r
202     Double_t tmpPhi=particle->Phi();\r
203     if(particle->Phi()> TMath::Pi()){\r
204       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
205     }\r
206 \r
207     \r
208     //process the gammas\r
209     if (particle->GetPdgCode()== 22){\r
210       fMCAllGammas.push_back(particle);\r
211       if(particle->GetMother(0)>-1){ //Means we have a mother\r
212         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.\r
213           fHistograms->FillHistogram("MC_Gamma_Energy", particle->Energy());\r
214           fHistograms->FillHistogram("MC_Gamma_Pt", particle->Pt());\r
215           fHistograms->FillHistogram("MC_Gamma_Eta", particle->Eta());\r
216           \r
217           fHistograms->FillHistogram("MC_Gamma_Phi", tmpPhi);\r
218 \r
219           //adding the conversion points from all gammas with e+e- daughters\r
220           if(particle->GetNDaughters() >= 2){\r
221             TParticle* daughter0 = NULL;\r
222             TParticle* daughter1 = NULL;\r
223             \r
224             for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
225               TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
226               if(tmpDaughter->GetUniqueID() == 5){\r
227                 if(tmpDaughter->GetPdgCode() == 11){\r
228                   daughter0 = tmpDaughter;\r
229                 }\r
230                 else if(tmpDaughter->GetPdgCode() == -11){\r
231                   daughter1 = tmpDaughter;\r
232                 }\r
233               }\r
234             }\r
235 \r
236             if(daughter0 == NULL || daughter1 == NULL){ // means we do not have two daughters from pair production\r
237               continue;\r
238             }\r
239 \r
240             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){\r
241               continue;\r
242             }\r
243             \r
244             if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||\r
245                (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
246 \r
247               // begin Mapping \r
248               Int_t rBin    = fHistograms->GetRBin(daughter0->R());\r
249               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());\r
250               \r
251               TString nameMCMappingPhiR="";\r
252               nameMCMappingPhiR.Form("MC_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
253               fHistograms->FillHistogram(nameMCMappingPhiR, daughter0->Vz(), particle->Eta());\r
254               \r
255               TString nameMCMappingPhi="";\r
256               nameMCMappingPhi.Form("MC_EP_Mapping-Phi%02d",phiBin);\r
257               fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
258               \r
259               TString nameMCMappingR="";\r
260               nameMCMappingR.Form("MC_EP_Mapping-R%02d",rBin);\r
261               fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
262               \r
263               TString nameMCMappingPhiInR="";\r
264               nameMCMappingPhiInR.Form("MC_EP_Mapping_Phi_R-%02d",rBin);\r
265               fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
266               //end mapping\r
267 \r
268               fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
269               fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
270               fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
271               fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
272             }// end if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||....... approx 20 lines above\r
273           }// end if(particle->GetNDaughters() >= 2){\r
274         } // end if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 )\r
275         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
276             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
277             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
278             ){ \r
279           fMCGammaChic.push_back(particle);\r
280         }\r
281       }// end if(particle->GetMother(0)>-1)\r
282       else{//means we have a primary particle\r
283         fHistograms->FillHistogram("MC_DirectGamma_Energy",particle->Energy());\r
284         fHistograms->FillHistogram("MC_DirectGamma_Pt", particle->Pt());\r
285         fHistograms->FillHistogram("MC_DirectGamma_Eta", particle->Eta());\r
286         fHistograms->FillHistogram("MCDirectGammaPhi", tmpPhi);\r
287 \r
288         //adding the conversion points from all gammas with e+e- daughters\r
289         if(particle->GetNDaughters() == 2){\r
290           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
291           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
292           if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11) ||\r
293              (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
294             \r
295             fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
296             fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
297             fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
298             fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
299 \r
300           }\r
301         }\r
302       }// end else\r
303     }// end if (particle->GetPdgCode()== 22){\r
304     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron\r
305       if(particle->GetMother(0)>-1){ // means we have a mother\r
306         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother\r
307           if(particle->GetPdgCode() == 11){//electron \r
308             fHistograms->FillHistogram("MC_E_Energy", particle->Energy());\r
309             fHistograms->FillHistogram("MC_E_Pt", particle->Pt());\r
310             fHistograms->FillHistogram("MC_E_Eta", particle->Eta());\r
311             fHistograms->FillHistogram("MC_E_Phi", tmpPhi);\r
312           }\r
313           if(particle->GetPdgCode() == -11){//positron \r
314             fHistograms->FillHistogram("MC_P_Energy", particle->Energy());\r
315             fHistograms->FillHistogram("MC_P_Pt", particle->Pt());\r
316             fHistograms->FillHistogram("MC_P_Eta", particle->Eta());\r
317             fHistograms->FillHistogram("MC_P_Phi", tmpPhi);\r
318           }\r
319         }\r
320       }\r
321     } // end else if (TMath::Abs(particle->GetPdgCode())== 11)\r
322     else if(particle->GetNDaughters() == 2){\r
323 \r
324       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
325       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
326       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters\r
327         \r
328         if(particle->GetPdgCode()==111){//Pi0\r
329           \r
330           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
331 \r
332           if( iTracks >= fStack->GetNprimary()){\r
333             \r
334             fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
335 \r
336             fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
337             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
338             fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
339             fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
340             fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
341             fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
342             fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
343           }\r
344           else{\r
345             fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());\r
346 \r
347             fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
348             fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
349             fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
350             fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
351             fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
352             fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
353             fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
354           }\r
355         }\r
356         else if(particle->GetPdgCode()==221){//Eta\r
357           fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
358 \r
359           fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
360           fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
361           fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
362           fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
363           fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
364           fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
365           fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
366         }\r
367         \r
368         //the match data should be filled no matter which mother the gamma-gamma comes from\r
369         fHistograms->FillHistogram("MC_Match_Gamma_R", particle->R());\r
370         fHistograms->FillHistogram("MC_Match_Gamma_ZR", particle->Vz(),particle->R());\r
371         fHistograms->FillHistogram("MC_Match_Gamma_XY", particle->Vx(),particle->Vy());\r
372         fHistograms->FillHistogram("MC_Match_Gamma_Mass", particle->GetCalcMass());\r
373         fHistograms->FillHistogram("MC_Match_Gamma_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
374         fHistograms->FillHistogram("MC_Match_Gamma_Energy", particle->Energy());\r
375         fHistograms->FillHistogram("MC_Match_Gamma_Pt", particle->Pt());\r
376         fHistograms->FillHistogram("MC_Match_Gamma_Eta", particle->Eta());\r
377         fHistograms->FillHistogram("MC_Match_Gamma_Phi",tmpPhi);\r
378       }\r
379     }// end else if(particle->GetNDaughters() == 2)\r
380   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
381 } // end ProcessMCData\r
382 \r
383 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
384   // see header file for documentation\r
385 \r
386   Int_t nSurvivingV0s=0;\r
387   while(fV0Reader->NextV0()){\r
388     nSurvivingV0s++;\r
389     //-------------------------- filling v0 information -------------------------------------\r
390     fHistograms->FillHistogram("ESD_EP_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
391     fHistograms->FillHistogram("ESD_EP_R", fV0Reader->GetXYRadius());\r
392     fHistograms->FillHistogram("ESD_EP_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
393     fHistograms->FillHistogram("ESD_EP_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
394     \r
395     \r
396     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
397     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
398     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
399     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
400     \r
401     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
402     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
403     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
404     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
405     \r
406     fHistograms->FillHistogram("ESD_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
407     fHistograms->FillHistogram("ESD_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
408     fHistograms->FillHistogram("ESD_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
409     fHistograms->FillHistogram("ESD_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
410 \r
411 \r
412     // begin mapping\r
413     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
414     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
415     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
416     \r
417     TString nameESDMappingPhiR="";\r
418     nameESDMappingPhiR.Form("ESD_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
419     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
420 \r
421     TString nameESDMappingPhi="";\r
422     nameESDMappingPhi.Form("ESD_EP_Mapping-Phi%02d",phiBin);\r
423     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
424 \r
425     TString nameESDMappingR="";\r
426     nameESDMappingR.Form("ESD_EP_Mapping-R%02d",rBin);\r
427     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
428 \r
429     TString nameESDMappingPhiInR="";\r
430     nameESDMappingPhiInR.Form("ESD_EP_Mapping_Phi_R-%02d",rBin);\r
431     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
432     // end mapping\r
433     \r
434     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
435 \r
436     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
437     if(fDoMCTruth){\r
438       if(fV0Reader->HasSameMCMother() == kFALSE){\r
439         continue;\r
440       }\r
441  \r
442       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
443       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
444 \r
445       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){\r
446         continue;\r
447       }\r
448       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
449         fHistograms->FillHistogram("ESD_Match_Gamma_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
450         fHistograms->FillHistogram("ESD_Match_Gamma_OpeningAngle", fV0Reader->GetOpeningAngle());\r
451         fHistograms->FillHistogram("ESD_Match_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
452         fHistograms->FillHistogram("ESD_Match_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
453         fHistograms->FillHistogram("ESD_Match_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
454 \r
455         fHistograms->FillHistogram("ESD_Match_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
456         fHistograms->FillHistogram("ESD_Match_Gamma_Mass", fV0Reader->GetMotherCandidateMass());\r
457         fHistograms->FillHistogram("ESD_Match_Gamma_Width", fV0Reader->GetMotherCandidateWidth());\r
458         fHistograms->FillHistogram("ESD_Match_Gamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
459         fHistograms->FillHistogram("ESD_Match_Gamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
460         fHistograms->FillHistogram("ESD_Match_Gamma_R", fV0Reader->GetXYRadius());\r
461         fHistograms->FillHistogram("ESD_Match_Gamma_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
462 \r
463         //resolution\r
464         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
465         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
466         Double_t resdPt = 0;\r
467         if(mcpt != 0){\r
468           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
469         }\r
470 \r
471         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
472         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
473         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
474         \r
475         Double_t resdZ = 0;\r
476         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
477           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
478         }\r
479         \r
480         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
481         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
482         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
483         \r
484         Double_t resdR = 0;\r
485         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
486           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
487         }\r
488         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
489         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
490         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
491         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
492       }\r
493     }\r
494   }\r
495   fHistograms->FillHistogram("NumberOfSurvivingV0s", nSurvivingV0s);\r
496   fHistograms->FillHistogram("NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
497 }\r
498 \r
499 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
500   // see header file for documentation\r
501 \r
502   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
503     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
504       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
505       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
506 \r
507       \r
508       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
509 \r
510       Double_t massTwoGammaCandidate =0.;\r
511       Double_t widthTwoGammaCandidate = 0.;\r
512       Double_t chi2TwoGammaCandidate =10000.;   \r
513       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
514       if(twoGammaCandidate->GetNDF()>0){\r
515         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
516         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){\r
517 \r
518           TVector3 vectorTwoGammaCandidate(twoGammaCandidate->Px(),twoGammaCandidate->Py(),twoGammaCandidate->Pz());\r
519 \r
520           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);\r
521 \r
522           //Calculating by hand the radius\r
523           Double_t tmpX= twoGammaCandidate->GetX();\r
524           Double_t tmpY= twoGammaCandidate->GetY();\r
525           \r
526           Double_t radiusTwoGammaCandidate = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
527 \r
528           fHistograms->FillHistogram("ESD_TwoGammaCombination_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
529           fHistograms->FillHistogram("ESD_TwoGammaCombination_Energy", twoGammaCandidate->GetE());\r
530           fHistograms->FillHistogram("ESD_TwoGammaCombination_Pt", sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
531           fHistograms->FillHistogram("ESD_TwoGammaCombination_Eta", vectorTwoGammaCandidate.Eta());\r
532           fHistograms->FillHistogram("ESD_TwoGammaCombination_Phi", vectorTwoGammaCandidate.Phi());\r
533           fHistograms->FillHistogram("ESD_TwoGammaCombination_Mass", massTwoGammaCandidate);\r
534           fHistograms->FillHistogram("ESD_TwoGammaCombination_R", radiusTwoGammaCandidate);\r
535           fHistograms->FillHistogram("ESD_TwoGammaCombination_ZR", tmpY, radiusTwoGammaCandidate);\r
536           fHistograms->FillHistogram("ESD_TwoGammaCombination_XY", tmpX, tmpY);\r
537           fHistograms->FillHistogram("InvMass_vs_Pt_Spectra",massTwoGammaCandidate ,sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
538         }\r
539       }\r
540       delete twoGammaCandidate;\r
541     }\r
542   }\r
543 }\r
544 \r
545 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
546   // see header file for documentation\r
547 \r
548   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
549   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
550   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
551     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
552     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
553       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
554 \r
555       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
556 \r
557       Double_t massBG =0.;\r
558       Double_t widthBG = 0.;\r
559       Double_t chi2BG =10000.;  \r
560       backgroundCandidate->GetMass(massBG,widthBG);\r
561       if(backgroundCandidate->GetNDF()>0){\r
562         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
563         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
564 \r
565           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());\r
566 \r
567           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
568 \r
569           //Calculating by hand the radius (find a better way)\r
570           Double_t tmpX= backgroundCandidate->GetX();\r
571           Double_t tmpY= backgroundCandidate->GetY();\r
572           \r
573           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
574 \r
575           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
576           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
577           fHistograms->FillHistogram("ESD_Background_Pt", sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
578           fHistograms->FillHistogram("ESD_Background_Eta", vectorBGCandidate.Eta());\r
579           fHistograms->FillHistogram("ESD_Background_Phi", vectorBGCandidate.Phi());\r
580           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
581           fHistograms->FillHistogram("ESD_Background_R", radiusBG);\r
582           fHistograms->FillHistogram("ESD_Background_ZR", tmpY, radiusBG);\r
583           fHistograms->FillHistogram("ESD_Background_XY", tmpX, tmpY);\r
584           fHistograms->FillHistogram("Background_InvMass_vs_Pt_Spectra",massBG,sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
585         }\r
586       }\r
587       delete backgroundCandidate;   \r
588     }\r
589   }\r
590 }\r
591 \r
592 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
593 {\r
594   // Terminate analysis\r
595   //\r
596   AliDebug(1,"Do nothing in Terminate");\r
597 }\r
598 \r
599 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
600 {\r
601   // Create the output container\r
602   if(fOutputContainer != NULL){\r
603     delete fOutputContainer;\r
604     fOutputContainer = NULL;\r
605   }\r
606   if(fOutputContainer == NULL){\r
607     fOutputContainer = new TList();\r
608   }\r
609   fHistograms->GetOutputContainer(fOutputContainer);\r
610   fOutputContainer->SetName(GetName());\r
611 }\r
612 \r
613 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
614   //helper function\r
615   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
616   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
617   return v3D0.Angle(v3D1);\r
618 }\r