Cleaning up files, added documentation. Added ntuple for v0s.
[u/mrichter/AliRoot.git] / PWG4 / GammaConv / 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.1                                                            *\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   fWriteNtuple(kFALSE),\r
69   fGammaNtuple(NULL),\r
70   fNeutralMesonNtuple(NULL),\r
71   fTotalNumberOfAddedNtupleEntries(0)\r
72 {\r
73   // Default constructor\r
74   // Common I/O in slot 0\r
75   DefineInput (0, TChain::Class());\r
76   DefineOutput(0, TTree::Class());\r
77 \r
78   // Your private output\r
79   DefineOutput(1, TList::Class());\r
80 }\r
81 \r
82 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
83   AliAnalysisTaskSE(name),\r
84   fV0Reader(NULL),\r
85   fStack(NULL),\r
86   fOutputContainer(0x0),\r
87   fHistograms(NULL),\r
88   fDoMCTruth(kFALSE),\r
89   fMCAllGammas(),\r
90   fMCPi0s(),\r
91   fMCEtas(),\r
92   fMCGammaChic(),\r
93   fKFReconstructedGammas(),\r
94   fElectronMass(-1),\r
95   fGammaMass(-1),\r
96   fPi0Mass(-1),\r
97   fEtaMass(-1),\r
98   fGammaWidth(-1),\r
99   fPi0Width(-1),\r
100   fEtaWidth(-1),\r
101   fCalculateBackground(kFALSE),\r
102   fWriteNtuple(kFALSE),\r
103   fGammaNtuple(NULL),\r
104   fNeutralMesonNtuple(NULL),\r
105   fTotalNumberOfAddedNtupleEntries(0)\r
106 {\r
107   // Common I/O in slot 0\r
108   DefineInput (0, TChain::Class());\r
109   DefineOutput(0, TTree::Class());\r
110   \r
111   // Your private output\r
112   DefineOutput(1, TList::Class());\r
113 }\r
114 \r
115 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
116 {\r
117   // Remove all pointers\r
118  \r
119   if(fOutputContainer){\r
120     fOutputContainer->Clear() ; \r
121     delete fOutputContainer ;\r
122   }\r
123   if(fHistograms){\r
124     delete fHistograms;\r
125   }\r
126   if(fV0Reader){\r
127     delete fV0Reader;\r
128   }\r
129 }\r
130 \r
131 \r
132 void AliAnalysisTaskGammaConversion::Init()\r
133 {\r
134   // Initialization\r
135   AliLog::SetGlobalLogLevel(AliLog::kError);\r
136 }\r
137 \r
138 \r
139 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
140 {\r
141   // Execute analysis for current event\r
142   \r
143   ConnectInputData("");\r
144   \r
145   //clear vectors\r
146   fMCAllGammas.clear();\r
147   fMCPi0s.clear();\r
148   fMCEtas.clear();\r
149   fMCGammaChic.clear();\r
150 \r
151   fKFReconstructedGammas.clear();\r
152 \r
153   //Clear the data in the v0Reader\r
154   fV0Reader->UpdateEventByEventData();\r
155 \r
156   // Process the MC information\r
157   if(fDoMCTruth){\r
158     ProcessMCData();\r
159   }\r
160 \r
161   // Process the v0 information\r
162   ProcessV0s();\r
163 \r
164   //calculate background if flag is set\r
165   if(fCalculateBackground){\r
166     CalculateBackground();\r
167   }\r
168 \r
169   // Process reconstructed gammas\r
170   ProcessGammasForNeutralMesonAnalysis();\r
171 \r
172   PostData(1, fOutputContainer);\r
173   \r
174 }\r
175 \r
176 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
177   // see header file for documentation\r
178 \r
179   if(fV0Reader == NULL){\r
180     // Write warning here cuts and so on are default if this ever happens\r
181   }\r
182   fV0Reader->Initialize();\r
183 }\r
184 \r
185 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
186   // see header file for documentation\r
187   \r
188   fStack = fV0Reader->GetMCStack();\r
189 \r
190   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
191     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
192 \r
193     if (!particle) {\r
194       //print warning here\r
195       continue;\r
196     }\r
197     \r
198     if(particle->Pt()<fV0Reader->GetPtCut()){\r
199       continue;\r
200     }\r
201     \r
202     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut()){\r
203       continue;\r
204     }\r
205 \r
206     if(particle->R()>fV0Reader->GetMaxRCut()){ // cuts on distance from collision point\r
207       continue;\r
208     }\r
209 \r
210     Double_t tmpPhi=particle->Phi();\r
211     if(particle->Phi()> TMath::Pi()){\r
212       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
213     }\r
214 \r
215     \r
216     //process the gammas\r
217     if (particle->GetPdgCode()== 22){\r
218       fMCAllGammas.push_back(particle);\r
219       if(particle->GetMother(0)>-1){ //Means we have a mother\r
220         if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 ){//Checks for a non gamma mother.\r
221           fHistograms->FillHistogram("MC_Gamma_Energy", particle->Energy());\r
222           fHistograms->FillHistogram("MC_Gamma_Pt", particle->Pt());\r
223           fHistograms->FillHistogram("MC_Gamma_Eta", particle->Eta());\r
224           \r
225           fHistograms->FillHistogram("MC_Gamma_Phi", tmpPhi);\r
226 \r
227           //adding the conversion points from all gammas with e+e- daughters\r
228           if(particle->GetNDaughters() >= 2){\r
229             TParticle* daughter0 = NULL;\r
230             TParticle* daughter1 = NULL;\r
231             \r
232             for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
233               TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
234               if(tmpDaughter->GetUniqueID() == 5){\r
235                 if(tmpDaughter->GetPdgCode() == 11){\r
236                   daughter0 = tmpDaughter;\r
237                 }\r
238                 else if(tmpDaughter->GetPdgCode() == -11){\r
239                   daughter1 = tmpDaughter;\r
240                 }\r
241               }\r
242             }\r
243 \r
244             if(daughter0 == NULL || daughter1 == NULL){ // means we do not have two daughters from pair production\r
245               continue;\r
246             }\r
247 \r
248             if(daughter0->R()>fV0Reader->GetMaxRCut() || daughter1->R()>fV0Reader->GetMaxRCut()){\r
249               continue;\r
250             }\r
251             \r
252             if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||\r
253                (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
254 \r
255               // begin Mapping \r
256               Int_t rBin    = fHistograms->GetRBin(daughter0->R());\r
257               Int_t phiBin  = fHistograms->GetPhiBin(daughter0->Phi());\r
258               \r
259               TString nameMCMappingPhiR="";\r
260               nameMCMappingPhiR.Form("MC_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
261               fHistograms->FillHistogram(nameMCMappingPhiR, daughter0->Vz(), particle->Eta());\r
262               \r
263               TString nameMCMappingPhi="";\r
264               nameMCMappingPhi.Form("MC_EP_Mapping-Phi%02d",phiBin);\r
265               fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
266               \r
267               TString nameMCMappingR="";\r
268               nameMCMappingR.Form("MC_EP_Mapping-R%02d",rBin);\r
269               fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
270               \r
271               TString nameMCMappingPhiInR="";\r
272               nameMCMappingPhiInR.Form("MC_EP_Mapping_Phi_R-%02d",rBin);\r
273               fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
274               //end mapping\r
275 \r
276               fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
277               fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
278               fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
279               fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
280             }// end if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode()) == 11 ||....... approx 20 lines above\r
281           }// end if(particle->GetNDaughters() >= 2){\r
282         } // end if( fStack->Particle(particle->GetMother(0))->GetPdgCode() != 22 )\r
283         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
284             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
285             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
286             ){ \r
287           fMCGammaChic.push_back(particle);\r
288         }\r
289       }// end if(particle->GetMother(0)>-1)\r
290       else{//means we have a primary particle\r
291         fHistograms->FillHistogram("MC_DirectGamma_Energy",particle->Energy());\r
292         fHistograms->FillHistogram("MC_DirectGamma_Pt", particle->Pt());\r
293         fHistograms->FillHistogram("MC_DirectGamma_Eta", particle->Eta());\r
294         fHistograms->FillHistogram("MCDirectGammaPhi", tmpPhi);\r
295 \r
296         //adding the conversion points from all gammas with e+e- daughters\r
297         if(particle->GetNDaughters() == 2){\r
298           TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
299           TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
300           if((daughter0->GetPdgCode() == -11 && daughter1->GetPdgCode() == 11) ||\r
301              (daughter0->GetPdgCode() == 11 && daughter1->GetPdgCode() == -11)){\r
302             \r
303             fHistograms->FillHistogram("MC_EP_R",daughter0->R());\r
304             fHistograms->FillHistogram("MC_EP_ZR",daughter0->Vz(),daughter0->R());\r
305             fHistograms->FillHistogram("MC_EP_XY",daughter0->Vx(),daughter0->Vy());\r
306             fHistograms->FillHistogram("MC_EP_OpeningAngle",GetMCOpeningAngle(daughter0, daughter1));\r
307 \r
308           }\r
309         }\r
310       }// end else\r
311     }// end if (particle->GetPdgCode()== 22){\r
312     else if (TMath::Abs(particle->GetPdgCode())== 11){ // Means we have an electron or a positron\r
313       if(particle->GetMother(0)>-1){ // means we have a mother\r
314         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==22 ){ // Means we have a gamma mother\r
315           if(particle->GetPdgCode() == 11){//electron \r
316             fHistograms->FillHistogram("MC_E_Energy", particle->Energy());\r
317             fHistograms->FillHistogram("MC_E_Pt", particle->Pt());\r
318             fHistograms->FillHistogram("MC_E_Eta", particle->Eta());\r
319             fHistograms->FillHistogram("MC_E_Phi", tmpPhi);\r
320           }\r
321           if(particle->GetPdgCode() == -11){//positron \r
322             fHistograms->FillHistogram("MC_P_Energy", particle->Energy());\r
323             fHistograms->FillHistogram("MC_P_Pt", particle->Pt());\r
324             fHistograms->FillHistogram("MC_P_Eta", particle->Eta());\r
325             fHistograms->FillHistogram("MC_P_Phi", tmpPhi);\r
326           }\r
327         }\r
328       }\r
329     } // end else if (TMath::Abs(particle->GetPdgCode())== 11)\r
330     else if(particle->GetNDaughters() == 2){\r
331 \r
332       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
333       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
334       if(daughter0->GetPdgCode() == 22 && daughter1->GetPdgCode() == 22){//check for gamma gamma daughters\r
335         \r
336         if(particle->GetPdgCode()==111){//Pi0\r
337           \r
338           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
339 \r
340           if( iTracks >= fStack->GetNprimary()){\r
341             \r
342             fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
343 \r
344             fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
345             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
346             fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
347             fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
348             fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
349             fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
350             fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
351           }\r
352           else{\r
353             fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());\r
354 \r
355             fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
356             fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
357             fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
358             fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
359             fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
360             fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
361             fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
362           }\r
363         }\r
364         else if(particle->GetPdgCode()==221){//Eta\r
365           fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
366 \r
367           fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
368           fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
369           fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
370           fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
371           fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
372           fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
373           fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
374         }\r
375         \r
376         //the match data should be filled no matter which mother the gamma-gamma comes from\r
377         fHistograms->FillHistogram("MC_Match_Gamma_R", particle->R());\r
378         fHistograms->FillHistogram("MC_Match_Gamma_ZR", particle->Vz(),particle->R());\r
379         fHistograms->FillHistogram("MC_Match_Gamma_XY", particle->Vx(),particle->Vy());\r
380         fHistograms->FillHistogram("MC_Match_Gamma_Mass", particle->GetCalcMass());\r
381         fHistograms->FillHistogram("MC_Match_Gamma_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
382         fHistograms->FillHistogram("MC_Match_Gamma_Energy", particle->Energy());\r
383         fHistograms->FillHistogram("MC_Match_Gamma_Pt", particle->Pt());\r
384         fHistograms->FillHistogram("MC_Match_Gamma_Eta", particle->Eta());\r
385         fHistograms->FillHistogram("MC_Match_Gamma_Phi",tmpPhi);\r
386       }\r
387     }// end else if(particle->GetNDaughters() == 2)\r
388   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
389 } // end ProcessMCData\r
390 \r
391 void AliAnalysisTaskGammaConversion::FillNtuple(){\r
392 \r
393   if(fGammaNtuple == NULL){\r
394     return;\r
395   }\r
396   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
397   for(Int_t i=0;i<numberOfV0s;i++){\r
398     Float_t values[27] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};\r
399     AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
400     Double_t negPID=0;\r
401     Double_t posPID=0;\r
402     fV0Reader->GetPIDProbability(negPID,posPID);\r
403     values[0]=cV0->GetOnFlyStatus();\r
404     values[1]=fV0Reader->CheckForPrimaryVertex();\r
405     values[2]=negPID;\r
406     values[3]=posPID;\r
407     values[4]=fV0Reader->GetX();\r
408     values[5]=fV0Reader->GetY();\r
409     values[6]=fV0Reader->GetZ();\r
410     values[7]=fV0Reader->GetXYRadius();\r
411     values[8]=fV0Reader->GetMotherCandidateNDF();\r
412     values[9]=fV0Reader->GetMotherCandidateChi2();\r
413     values[10]=fV0Reader->GetMotherCandidateEnergy();\r
414     values[11]=fV0Reader->GetMotherCandidateEta();\r
415     values[12]=fV0Reader->GetMotherCandidatePt();\r
416     values[13]=fV0Reader->GetMotherCandidateMass();\r
417     values[14]=fV0Reader->GetMotherCandidateWidth();\r
418     //    values[15]=fV0Reader->GetMotherMCParticle()->Pt();   MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array\r
419     values[16]=fV0Reader->GetOpeningAngle();\r
420     values[17]=fV0Reader->GetNegativeTrackEnergy();\r
421     values[18]=fV0Reader->GetNegativeTrackPt();\r
422     values[19]=fV0Reader->GetNegativeTrackEta();\r
423     values[20]=fV0Reader->GetNegativeTrackPhi();\r
424     values[21]=fV0Reader->GetPositiveTrackEnergy();\r
425     values[22]=fV0Reader->GetPositiveTrackPt();\r
426     values[23]=fV0Reader->GetPositiveTrackEta();\r
427     values[24]=fV0Reader->GetPositiveTrackPhi();\r
428     values[25]=fV0Reader->HasSameMCMother();\r
429     if(values[25] != 0){\r
430       values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
431       values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
432     }\r
433     fTotalNumberOfAddedNtupleEntries++;\r
434     fGammaNtuple->Fill(values);\r
435   }\r
436   fV0Reader->ResetV0IndexNumber();\r
437   \r
438 }\r
439 \r
440 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
441   // see header file for documentation\r
442 \r
443   if(fWriteNtuple == kTRUE){\r
444     FillNtuple();\r
445   }\r
446 \r
447   Int_t nSurvivingV0s=0;\r
448   while(fV0Reader->NextV0()){\r
449     nSurvivingV0s++;\r
450     //-------------------------- filling v0 information -------------------------------------\r
451     fHistograms->FillHistogram("ESD_EP_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
452     fHistograms->FillHistogram("ESD_EP_R", fV0Reader->GetXYRadius());\r
453     fHistograms->FillHistogram("ESD_EP_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
454     fHistograms->FillHistogram("ESD_EP_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
455     \r
456     \r
457     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
458     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
459     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
460     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
461     \r
462     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
463     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
464     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
465     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
466     \r
467     fHistograms->FillHistogram("ESD_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
468     fHistograms->FillHistogram("ESD_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
469     fHistograms->FillHistogram("ESD_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
470     fHistograms->FillHistogram("ESD_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
471 \r
472 \r
473     // begin mapping\r
474     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
475     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
476     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
477     \r
478     TString nameESDMappingPhiR="";\r
479     nameESDMappingPhiR.Form("ESD_EP_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
480     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
481 \r
482     TString nameESDMappingPhi="";\r
483     nameESDMappingPhi.Form("ESD_EP_Mapping-Phi%02d",phiBin);\r
484     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
485 \r
486     TString nameESDMappingR="";\r
487     nameESDMappingR.Form("ESD_EP_Mapping-R%02d",rBin);\r
488     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
489 \r
490     TString nameESDMappingPhiInR="";\r
491     nameESDMappingPhiInR.Form("ESD_EP_Mapping_Phi_R-%02d",rBin);\r
492     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
493     // end mapping\r
494     \r
495     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
496 \r
497     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
498     if(fDoMCTruth){\r
499       if(fV0Reader->HasSameMCMother() == kFALSE){\r
500         continue;\r
501       }\r
502  \r
503       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
504       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
505 \r
506       if(negativeMC->GetPdgCode()!=11 || positiveMC->GetPdgCode()!=-11){\r
507         continue;\r
508       }\r
509       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
510         fHistograms->FillHistogram("ESD_Match_Gamma_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
511         fHistograms->FillHistogram("ESD_Match_Gamma_OpeningAngle", fV0Reader->GetOpeningAngle());\r
512         fHistograms->FillHistogram("ESD_Match_Gamma_Pt", fV0Reader->GetMotherCandidatePt());\r
513         fHistograms->FillHistogram("ESD_Match_Gamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
514         fHistograms->FillHistogram("ESD_Match_Gamma_Eta", fV0Reader->GetMotherCandidateEta());\r
515 \r
516         fHistograms->FillHistogram("ESD_Match_Gamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
517         fHistograms->FillHistogram("ESD_Match_Gamma_Mass", fV0Reader->GetMotherCandidateMass());\r
518         fHistograms->FillHistogram("ESD_Match_Gamma_Width", fV0Reader->GetMotherCandidateWidth());\r
519         fHistograms->FillHistogram("ESD_Match_Gamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
520         fHistograms->FillHistogram("ESD_Match_Gamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
521         fHistograms->FillHistogram("ESD_Match_Gamma_R", fV0Reader->GetXYRadius());\r
522         fHistograms->FillHistogram("ESD_Match_Gamma_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
523 \r
524         //resolution\r
525         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
526         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
527         Double_t resdPt = 0;\r
528         if(mcpt != 0){\r
529           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
530         }\r
531 \r
532         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
533         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
534         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
535         \r
536         Double_t resdZ = 0;\r
537         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
538           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
539         }\r
540         \r
541         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
542         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
543         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
544         \r
545         Double_t resdR = 0;\r
546         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
547           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
548         }\r
549         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
550         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
551         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
552         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
553       }\r
554     }\r
555   }\r
556   fHistograms->FillHistogram("NumberOfSurvivingV0s", nSurvivingV0s);\r
557   fHistograms->FillHistogram("NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
558 }\r
559 \r
560 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
561   // see header file for documentation\r
562 \r
563   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
564     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
565       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
566       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
567 \r
568       \r
569       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
570 \r
571       Double_t massTwoGammaCandidate =0.;\r
572       Double_t widthTwoGammaCandidate = 0.;\r
573       Double_t chi2TwoGammaCandidate =10000.;   \r
574       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
575       if(twoGammaCandidate->GetNDF()>0){\r
576         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
577         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){\r
578 \r
579           TVector3 vectorTwoGammaCandidate(twoGammaCandidate->Px(),twoGammaCandidate->Py(),twoGammaCandidate->Pz());\r
580 \r
581           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);\r
582 \r
583           //Calculating by hand the radius\r
584           Double_t tmpX= twoGammaCandidate->GetX();\r
585           Double_t tmpY= twoGammaCandidate->GetY();\r
586           \r
587           Double_t radiusTwoGammaCandidate = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
588 \r
589           fHistograms->FillHistogram("ESD_TwoGammaCombination_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
590           fHistograms->FillHistogram("ESD_TwoGammaCombination_Energy", twoGammaCandidate->GetE());\r
591           fHistograms->FillHistogram("ESD_TwoGammaCombination_Pt", sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
592           fHistograms->FillHistogram("ESD_TwoGammaCombination_Eta", vectorTwoGammaCandidate.Eta());\r
593           fHistograms->FillHistogram("ESD_TwoGammaCombination_Phi", vectorTwoGammaCandidate.Phi());\r
594           fHistograms->FillHistogram("ESD_TwoGammaCombination_Mass", massTwoGammaCandidate);\r
595           fHistograms->FillHistogram("ESD_TwoGammaCombination_R", radiusTwoGammaCandidate);\r
596           fHistograms->FillHistogram("ESD_TwoGammaCombination_ZR", tmpY, radiusTwoGammaCandidate);\r
597           fHistograms->FillHistogram("ESD_TwoGammaCombination_XY", tmpX, tmpY);\r
598           fHistograms->FillHistogram("InvMass_vs_Pt_Spectra",massTwoGammaCandidate ,sqrt(twoGammaCandidate->GetPx()*twoGammaCandidate->GetPx()+twoGammaCandidate->GetPy()*twoGammaCandidate->GetPy()));\r
599         }\r
600       }\r
601       delete twoGammaCandidate;\r
602     }\r
603   }\r
604 }\r
605 \r
606 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
607   // see header file for documentation\r
608 \r
609   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
610   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
611   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
612     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
613     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
614       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
615 \r
616       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
617 \r
618       Double_t massBG =0.;\r
619       Double_t widthBG = 0.;\r
620       Double_t chi2BG =10000.;  \r
621       backgroundCandidate->GetMass(massBG,widthBG);\r
622       if(backgroundCandidate->GetNDF()>0){\r
623         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
624         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
625 \r
626           TVector3 vectorBGCandidate(backgroundCandidate->Px(),backgroundCandidate->Py(),backgroundCandidate->Pz());\r
627 \r
628           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
629 \r
630           //Calculating by hand the radius (find a better way)\r
631           Double_t tmpX= backgroundCandidate->GetX();\r
632           Double_t tmpY= backgroundCandidate->GetY();\r
633           \r
634           Double_t radiusBG = TMath::Sqrt(tmpX*tmpX + tmpY*tmpY);\r
635 \r
636           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
637           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
638           fHistograms->FillHistogram("ESD_Background_Pt", sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
639           fHistograms->FillHistogram("ESD_Background_Eta", vectorBGCandidate.Eta());\r
640           fHistograms->FillHistogram("ESD_Background_Phi", vectorBGCandidate.Phi());\r
641           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
642           fHistograms->FillHistogram("ESD_Background_R", radiusBG);\r
643           fHistograms->FillHistogram("ESD_Background_ZR", tmpY, radiusBG);\r
644           fHistograms->FillHistogram("ESD_Background_XY", tmpX, tmpY);\r
645           fHistograms->FillHistogram("Background_InvMass_vs_Pt_Spectra",massBG,sqrt(backgroundCandidate->GetPx()*backgroundCandidate->GetPx()+backgroundCandidate->GetPy()*backgroundCandidate->GetPy()));\r
646         }\r
647       }\r
648       delete backgroundCandidate;   \r
649     }\r
650   }\r
651 }\r
652 \r
653 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
654 {\r
655   // Terminate analysis\r
656   //\r
657   AliDebug(1,"Do nothing in Terminate");\r
658 }\r
659 \r
660 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
661 {\r
662   // Create the output container\r
663   if(fOutputContainer != NULL){\r
664     delete fOutputContainer;\r
665     fOutputContainer = NULL;\r
666   }\r
667   if(fOutputContainer == NULL){\r
668     fOutputContainer = new TList();\r
669   }\r
670   \r
671   //Adding the histograms to the output container\r
672   fHistograms->GetOutputContainer(fOutputContainer);\r
673 \r
674   \r
675   if(fWriteNtuple){\r
676     if(fGammaNtuple == NULL){\r
677       fGammaNtuple = new TNtuple("V0ntuple","V0ntuple","OnTheFly:HasVertex:NegPIDProb:PosPIDProb:X:Y:Z:R:MotherCandidateNDF:MotherCandidateChi2:MotherCandidateEnergy:MotherCandidateEta:MotherCandidatePt:MotherCandidateMass:MotherCandidateWidth:MCMotherCandidatePT:EPOpeningAngle:ElectronEnergy:ElectronPt:ElectronEta:ElectronPhi:PositronEnergy:PositronPt:PositronEta:PositronPhi:HasSameMCMother:MotherMCParticlePIDCode",50000);\r
678     }\r
679     if(fNeutralMesonNtuple == NULL){\r
680       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
681     }\r
682     TList * ntupleTList = new TList();\r
683     ntupleTList->SetName("Ntuple");\r
684     ntupleTList->Add((TNtuple*)fGammaNtuple);\r
685     fOutputContainer->Add(ntupleTList);\r
686   }\r
687 \r
688   fOutputContainer->SetName(GetName());\r
689 }\r
690 \r
691 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
692   //helper function\r
693   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
694   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
695   return v3D0.Angle(v3D1);\r
696 }\r