Changed to new naming scheme for histograms, added functionality to the analysis...
[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   fIsTrueReconstructedGammas(),\r
61   electronv1(),\r
62   electronv2(),\r
63   fElectronMass(-1),\r
64   fGammaMass(-1),\r
65   fPi0Mass(-1),\r
66   fEtaMass(-1),\r
67   fGammaWidth(-1),\r
68   fPi0Width(-1),\r
69   fEtaWidth(-1),\r
70   fMinOpeningAngleGhostCut(0.),\r
71   fCalculateBackground(kFALSE),\r
72   fWriteNtuple(kFALSE),\r
73   fGammaNtuple(NULL),\r
74   fNeutralMesonNtuple(NULL),\r
75   fTotalNumberOfAddedNtupleEntries(0)\r
76 {\r
77   // Default constructor\r
78   // Common I/O in slot 0\r
79   DefineInput (0, TChain::Class());\r
80   DefineOutput(0, TTree::Class());\r
81         \r
82   // Your private output\r
83   DefineOutput(1, TList::Class());\r
84 }\r
85 \r
86 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
87   AliAnalysisTaskSE(name),\r
88   fV0Reader(NULL),\r
89   fStack(NULL),\r
90   fOutputContainer(0x0),\r
91   fHistograms(NULL),\r
92   fDoMCTruth(kFALSE),\r
93   fMCAllGammas(),\r
94   fMCPi0s(),\r
95   fMCEtas(),\r
96   fMCGammaChic(),\r
97   fKFReconstructedGammas(),\r
98   fIsTrueReconstructedGammas(),\r
99   electronv1(),\r
100   electronv2(),\r
101   fElectronMass(-1),\r
102   fGammaMass(-1),\r
103   fPi0Mass(-1),\r
104   fEtaMass(-1),\r
105   fGammaWidth(-1),\r
106   fPi0Width(-1),\r
107   fEtaWidth(-1),\r
108   fMinOpeningAngleGhostCut(0.),\r
109   fCalculateBackground(kFALSE),\r
110   fWriteNtuple(kFALSE),\r
111   fGammaNtuple(NULL),\r
112   fNeutralMesonNtuple(NULL),\r
113   fTotalNumberOfAddedNtupleEntries(0)\r
114 {\r
115   // Common I/O in slot 0\r
116   DefineInput (0, TChain::Class());\r
117   DefineOutput(0, TTree::Class());\r
118         \r
119   // Your private output\r
120   DefineOutput(1, TList::Class());\r
121 }\r
122 \r
123 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
124 {\r
125   // Remove all pointers\r
126         \r
127   if(fOutputContainer){\r
128     fOutputContainer->Clear() ; \r
129     delete fOutputContainer ;\r
130   }\r
131   if(fHistograms){\r
132     delete fHistograms;\r
133   }\r
134   if(fV0Reader){\r
135     delete fV0Reader;\r
136   }\r
137 }\r
138 \r
139 \r
140 void AliAnalysisTaskGammaConversion::Init()\r
141 {\r
142   // Initialization\r
143   AliLog::SetGlobalLogLevel(AliLog::kError);\r
144 }\r
145 \r
146 \r
147 void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
148 {\r
149   // Execute analysis for current event\r
150         \r
151   ConnectInputData("");\r
152         \r
153   //clear vectors\r
154   fMCAllGammas.clear();\r
155   fMCPi0s.clear();\r
156   fMCEtas.clear();\r
157   fMCGammaChic.clear();\r
158         \r
159   fKFReconstructedGammas.clear();\r
160   fIsTrueReconstructedGammas.clear();\r
161   electronv1.clear();\r
162   electronv2.clear();\r
163         \r
164   //Clear the data in the v0Reader\r
165   fV0Reader->UpdateEventByEventData();\r
166 \r
167   \r
168   // Process the MC information\r
169   if(fDoMCTruth){\r
170     ProcessMCData();\r
171   }\r
172   \r
173   //Process the v0 information with no cuts\r
174   ProcessV0sNoCut();\r
175   \r
176   // Process the v0 information\r
177   ProcessV0s();\r
178   \r
179   //calculate background if flag is set\r
180   if(fCalculateBackground){\r
181     CalculateBackground();\r
182   }\r
183   \r
184   // Process reconstructed gammas\r
185   ProcessGammasForNeutralMesonAnalysis();\r
186   \r
187   PostData(1, fOutputContainer);\r
188         \r
189 }\r
190 \r
191 void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
192   // see header file for documentation\r
193         \r
194   if(fV0Reader == NULL){\r
195     // Write warning here cuts and so on are default if this ever happens\r
196   }\r
197   fV0Reader->Initialize();\r
198 }\r
199 \r
200 \r
201 \r
202 void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
203   // see header file for documentation\r
204         \r
205   fStack = fV0Reader->GetMCStack();\r
206 \r
207   if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
208     return; // aborts if the primary vertex does not have contributors.\r
209   }\r
210 \r
211   for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
212     TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
213                 \r
214     if (!particle) {\r
215       //print warning here\r
216       continue;\r
217     }\r
218 \r
219     if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() )    continue;\r
220                                         \r
221     if(particle->R()>fV0Reader->GetMaxRCut())   continue; // cuts on distance from collision point\r
222                 \r
223     Double_t tmpPhi=particle->Phi();\r
224                 \r
225     if(particle->Phi()> TMath::Pi()){\r
226       tmpPhi = particle->Phi()-(2*TMath::Pi());\r
227     }\r
228                 \r
229     Double_t rapidity;\r
230     if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
231       rapidity=0;\r
232     }\r
233     else{\r
234       rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
235     }   \r
236                 \r
237     //process the gammas\r
238     if (particle->GetPdgCode() == 22){\r
239                         \r
240       if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
241         continue; // no photon as mothers!\r
242       }\r
243 \r
244       if(particle->GetMother(0) >= fStack->GetNprimary()){\r
245         continue; // the gamma has a mother, and it is not a primary particle\r
246       }\r
247 \r
248       fMCAllGammas.push_back(particle);\r
249                         \r
250       fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
251       fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
252       fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
253       fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
254       fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
255                         \r
256                         \r
257       if(particle->GetMother(0) < 0){   // direct gamma\r
258         fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
259         fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
260         fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
261         fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
262         fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);                                \r
263       }\r
264                         \r
265                         \r
266       // looking for conversion (electron + positron from pairbuilding (= 5) )\r
267       TParticle* ePos = NULL;\r
268       TParticle* eNeg = NULL;\r
269                         \r
270       if(particle->GetNDaughters() >= 2){\r
271         for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
272           TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
273           if(tmpDaughter->GetUniqueID() == 5){\r
274             if(tmpDaughter->GetPdgCode() == 11){\r
275               eNeg = tmpDaughter;\r
276             }\r
277             else if(tmpDaughter->GetPdgCode() == -11){\r
278               ePos = tmpDaughter;\r
279             }\r
280           }\r
281         }\r
282       }\r
283                         \r
284                         \r
285       if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
286         continue;\r
287       }\r
288                         \r
289                         \r
290       Double_t ePosPhi = ePos->Phi();\r
291       if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
292                         \r
293       Double_t eNegPhi = eNeg->Phi();\r
294       if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
295                         \r
296                         \r
297       if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
298         continue; // no reconstruction below the Pt cut\r
299       }\r
300                                         \r
301       if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
302         continue;\r
303       } \r
304                                 \r
305       if(ePos->R()>fV0Reader->GetMaxRCut()){\r
306         continue; // cuts on distance from collision point\r
307       }\r
308       \r
309       \r
310       if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  > ePos->R()){\r
311         continue;               // line cut to exclude regions where we do not reconstruct\r
312       }         \r
313                 \r
314       fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
315       fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
316       fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
317       fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
318       fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
319       fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
320                         \r
321       fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
322       fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
323       fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
324       fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
325                         \r
326       fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
327       fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
328       fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
329       fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
330                         \r
331                         \r
332                         \r
333       //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
334                         \r
335       // begin Mapping \r
336       Int_t rBin    = fHistograms->GetRBin(ePos->R());\r
337       Int_t phiBin  = fHistograms->GetPhiBin(particle->Phi());\r
338                         \r
339       TString nameMCMappingPhiR="";\r
340       nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
341       fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
342                         \r
343       TString nameMCMappingPhi="";\r
344       nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
345       fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
346                         \r
347       TString nameMCMappingR="";\r
348       nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
349       fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
350                         \r
351       TString nameMCMappingPhiInR="";\r
352       nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
353       fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
354       //end mapping\r
355                         \r
356       fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
357       fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
358       fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
359       fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
360                         \r
361       //cout << "mapping is done" << endl;\r
362                         \r
363                         \r
364       if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
365         fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
366         fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
367         fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
368         fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
369         fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
370                                 \r
371       } // end direct gamma\r
372       else{   // mother exits \r
373         if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
374             fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
375             fStack->Particle(particle->GetMother(0))->GetPdgCode()==445  //chic2\r
376             ){ \r
377           fMCGammaChic.push_back(particle);\r
378         }\r
379       }  // end if mother exits\r
380     } // end if particle is a photon\r
381                 \r
382     if(particle->GetNDaughters() == 2){\r
383                         \r
384       TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
385       TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
386                         \r
387       if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
388                         \r
389                         \r
390                         \r
391       // check for conversions now -> have to pass eta and line cut!\r
392       Bool_t daughter0Electron = kFALSE;\r
393       Bool_t daughter0Positron = kFALSE;\r
394       Bool_t daughter1Electron = kFALSE;\r
395       Bool_t daughter1Positron = kFALSE;\r
396                         \r
397                         \r
398                         \r
399       if(daughter0->GetNDaughters() >= 2){\r
400         for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
401           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
402           if(tmpDaughter->GetUniqueID() == 5){\r
403             if(tmpDaughter->GetPdgCode() == 11){\r
404               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
405                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
406                   daughter0Electron = kTRUE;\r
407                 }\r
408                                                                 \r
409               }\r
410             }\r
411             else if(tmpDaughter->GetPdgCode() == -11){\r
412               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
413                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
414                   daughter0Positron = kTRUE;\r
415                 }\r
416                                                                 \r
417               }\r
418                                                         \r
419             }\r
420           }\r
421         }\r
422       }\r
423                         \r
424                         \r
425         \r
426       if(daughter1->GetNDaughters() >= 2){\r
427         for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
428           TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
429           if(tmpDaughter->GetUniqueID() == 5){\r
430             if(tmpDaughter->GetPdgCode() == 11){\r
431               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
432                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
433                   daughter1Electron = kTRUE;\r
434                 }\r
435                                                                 \r
436               }\r
437             }\r
438             else if(tmpDaughter->GetPdgCode() == -11){\r
439               if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
440                 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue()  < tmpDaughter->R() ){\r
441                   daughter1Positron = kTRUE;\r
442                 }\r
443                                                                 \r
444               }\r
445                                                         \r
446             }\r
447           }\r
448         }\r
449       }\r
450                         \r
451                                                                                                 \r
452                         \r
453                         \r
454       if(particle->GetPdgCode()==111){     //Pi0\r
455         if( iTracks >= fStack->GetNprimary()){\r
456           fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
457           fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
458           fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
459           fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
460           fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
461           fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
462           fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
463           fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
464           fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
465                                         \r
466           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
467             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
468             fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
469             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
470               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
471               fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
472             }\r
473           }\r
474         }\r
475         else{\r
476           fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());    \r
477           fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
478           fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
479           fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
480           fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
481           fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
482           fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
483           fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
484           fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
485                                         \r
486           if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
487             fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
488             fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
489             if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
490               fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
491               fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
492             }\r
493           }\r
494         }\r
495       }\r
496                         \r
497       if(particle->GetPdgCode()==221){   //Eta\r
498         fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
499         fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
500         fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
501         fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
502         fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
503         fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
504         fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
505         fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
506         fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
507                                 \r
508         if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
509           fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
510           fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
511           if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
512             fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
513             fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
514           }\r
515                                         \r
516         }\r
517                                 \r
518       }\r
519                         \r
520       // all motherparticles with 2 gammas as daughters\r
521       fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
522       fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
523       fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
524       fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
525       fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
526       fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
527       fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
528       fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
529       fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
530       fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
531       fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());                 \r
532       if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
533         fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
534         fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
535         fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());                      \r
536         if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
537           fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
538           fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
539           fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());                  \r
540 \r
541         }\r
542                                 \r
543                                 \r
544       }\r
545                         \r
546       //cout << "mother histos are filled" << endl;\r
547                         \r
548     } // end if(particle->GetNDaughters() == 2)\r
549                 \r
550   }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
551         \r
552   //cout << "right before the end of processMCdata" << endl;\r
553         \r
554 } // end ProcessMCData\r
555 \r
556 \r
557 \r
558 void AliAnalysisTaskGammaConversion::FillNtuple(){\r
559         \r
560   if(fGammaNtuple == NULL){\r
561     return;\r
562   }\r
563   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
564   for(Int_t i=0;i<numberOfV0s;i++){\r
565     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
566     AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
567     Double_t negPID=0;\r
568     Double_t posPID=0;\r
569     fV0Reader->GetPIDProbability(negPID,posPID);\r
570     values[0]=cV0->GetOnFlyStatus();\r
571     values[1]=fV0Reader->CheckForPrimaryVertex();\r
572     values[2]=negPID;\r
573     values[3]=posPID;\r
574     values[4]=fV0Reader->GetX();\r
575     values[5]=fV0Reader->GetY();\r
576     values[6]=fV0Reader->GetZ();\r
577     values[7]=fV0Reader->GetXYRadius();\r
578     values[8]=fV0Reader->GetMotherCandidateNDF();\r
579     values[9]=fV0Reader->GetMotherCandidateChi2();\r
580     values[10]=fV0Reader->GetMotherCandidateEnergy();\r
581     values[11]=fV0Reader->GetMotherCandidateEta();\r
582     values[12]=fV0Reader->GetMotherCandidatePt();\r
583     values[13]=fV0Reader->GetMotherCandidateMass();\r
584     values[14]=fV0Reader->GetMotherCandidateWidth();\r
585     //    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
586     values[16]=fV0Reader->GetOpeningAngle();\r
587     values[17]=fV0Reader->GetNegativeTrackEnergy();\r
588     values[18]=fV0Reader->GetNegativeTrackPt();\r
589     values[19]=fV0Reader->GetNegativeTrackEta();\r
590     values[20]=fV0Reader->GetNegativeTrackPhi();\r
591     values[21]=fV0Reader->GetPositiveTrackEnergy();\r
592     values[22]=fV0Reader->GetPositiveTrackPt();\r
593     values[23]=fV0Reader->GetPositiveTrackEta();\r
594     values[24]=fV0Reader->GetPositiveTrackPhi();\r
595     values[25]=fV0Reader->HasSameMCMother();\r
596     if(values[25] != 0){\r
597       values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
598       values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
599     }\r
600     fTotalNumberOfAddedNtupleEntries++;\r
601     fGammaNtuple->Fill(values);\r
602   }\r
603   fV0Reader->ResetV0IndexNumber();\r
604         \r
605 }\r
606 \r
607 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
608 \r
609   Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
610   for(Int_t i=0;i<numberOfV0s;i++){\r
611     /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
612 \r
613     if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
614       return;\r
615     }\r
616     \r
617     if(fDoMCTruth){\r
618       \r
619       if(fV0Reader->HasSameMCMother() == kFALSE){\r
620         continue;\r
621       }\r
622                 \r
623       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
624       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
625 \r
626       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
627         continue;\r
628       }\r
629       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
630         continue;\r
631       }\r
632         \r
633       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
634       \r
635         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
636         fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
637         fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                               \r
638         fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
639         fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
640         fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
641         fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
642         fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
643         fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
644         fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
645         \r
646         fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
647         fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
648         fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
649         fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
650         \r
651         /*\r
652           ESD_NoCutConvGamma_Pt\r
653           ESD_NoCutConvGamma_Energy\r
654           ESD_NoCutConvGamma_Eta\r
655           ESD_NoCutConvGamma_Phi\r
656           ESD_NoCutConvGamma_Mass\r
657           ESD_NoCutConvGamma_Width\r
658           ESD_NoCutConvGamma_Chi2\r
659           ESD_NoCutConvGamma_NDF\r
660           ESD_NoCutConvGamma_PtvsEta\r
661           ESD_NoCutConversion_XY\r
662           ESD_NoCutConversion_R\r
663           ESD_NoCutConversion_ZR\r
664           ESD_NoCutConversion_OpeningAngle\r
665         */\r
666       }\r
667     }\r
668   }\r
669   fV0Reader->ResetV0IndexNumber();\r
670 }\r
671 \r
672 void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
673   // see header file for documentation\r
674         \r
675   if(fWriteNtuple == kTRUE){\r
676     FillNtuple();\r
677   }\r
678         \r
679   Int_t nSurvivingV0s=0;\r
680   while(fV0Reader->NextV0()){\r
681     nSurvivingV0s++;\r
682                 \r
683                 \r
684     //-------------------------- filling v0 information -------------------------------------\r
685     fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
686     fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
687     fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
688     fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());    \r
689                 \r
690     fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
691     fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
692     fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
693     fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
694                 \r
695     fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
696     fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
697     fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
698     fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
699                 \r
700     fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
701     fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
702     fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
703     fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
704     fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
705     fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
706     fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
707     fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
708     fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
709     fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
710                 \r
711                 \r
712     // begin mapping\r
713     Int_t rBin    = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
714     Int_t phiBin  = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
715     Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
716                 \r
717     TString nameESDMappingPhiR="";\r
718     nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
719     fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
720                 \r
721     TString nameESDMappingPhi="";\r
722     nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
723     fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
724                 \r
725     TString nameESDMappingR="";\r
726     nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
727     fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);  \r
728                 \r
729     TString nameESDMappingPhiInR="";\r
730     nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
731     fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
732     // end mapping\r
733                 \r
734     fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
735     electronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
736     electronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
737 \r
738                 \r
739     //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
740     if(fDoMCTruth){\r
741                         \r
742       if(fV0Reader->HasSameMCMother() == kFALSE){\r
743         fIsTrueReconstructedGammas.push_back(kFALSE);\r
744         continue;\r
745       }\r
746       \r
747                 \r
748       TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
749       TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
750 \r
751       if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
752         fIsTrueReconstructedGammas.push_back(kFALSE);\r
753         continue;\r
754       }\r
755       if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
756         fIsTrueReconstructedGammas.push_back(kFALSE);\r
757         continue;\r
758       }\r
759         \r
760       if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
761         fIsTrueReconstructedGammas.push_back(kTRUE);\r
762                                 \r
763         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
764         fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
765         fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());                                \r
766         fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
767         fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
768         fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
769         fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
770         fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
771         fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
772         fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
773         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
774         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
775         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
776         fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
777         \r
778         fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
779         fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
780         fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
781         fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
782 \r
783         \r
784         /*\r
785           fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
786           fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
787           fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
788           fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
789         */\r
790 \r
791 \r
792 \r
793         //resolution\r
794         Double_t mcpt   = fV0Reader->GetMotherMCParticle()->Pt();\r
795         Double_t esdpt  = fV0Reader->GetMotherCandidatePt();\r
796         Double_t resdPt = 0;\r
797         if(mcpt != 0){\r
798           resdPt = ((esdpt - mcpt)/mcpt)*100;\r
799         }\r
800                                 \r
801         fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
802         fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
803         fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
804                                 \r
805         Double_t resdZ = 0;\r
806         if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
807           resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
808         }\r
809                                 \r
810         fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
811         fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
812         fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
813                                 \r
814         Double_t resdR = 0;\r
815         if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
816           resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
817         }\r
818 \r
819         fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
820         fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
821         fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
822         fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
823       }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
824       else{\r
825         fIsTrueReconstructedGammas.push_back(kFALSE);\r
826       }\r
827     }//if(fDoMCTruth)\r
828   }//while(fV0Reader->NextV0)\r
829   fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
830   fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
831         \r
832   //cout << "nearly at the end of doMCTruth" << endl;\r
833         \r
834 }\r
835 \r
836 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
837   // see header file for documentation\r
838         \r
839   for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
840     for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
841                         \r
842       AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
843       AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
844       \r
845       if(electronv1[firstGammaIndex]==electronv1[secondGammaIndex] || electronv1[firstGammaIndex]==electronv2[secondGammaIndex]){\r
846         continue;\r
847       }\r
848       if(electronv2[firstGammaIndex]==electronv1[secondGammaIndex] || electronv2[firstGammaIndex]==electronv2[secondGammaIndex]){\r
849         continue;\r
850       }\r
851 \r
852       /*\r
853         if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
854         continue;\r
855         }\r
856       */\r
857 \r
858       AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
859                         \r
860       Double_t massTwoGammaCandidate = 0.;\r
861       Double_t widthTwoGammaCandidate = 0.;\r
862       Double_t chi2TwoGammaCandidate =10000.;   \r
863       twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
864       if(twoGammaCandidate->GetNDF()>0){\r
865         chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
866                                 \r
867         if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){                                      \r
868                                         \r
869           TVector3 MomentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
870           TVector3 SpaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
871                                         \r
872           Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);                                 \r
873           Double_t rapidity;\r
874           if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
875             rapidity=0;\r
876           }\r
877           else{\r
878             rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
879           }\r
880                                         \r
881           if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue;   // minimum opening angle to avoid using ghosttracks\r
882 \r
883           fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
884           fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
885           fHistograms->FillHistogram("ESD_Mother_Pt", MomentumVectorTwoGammaCandidate.Pt());\r
886           fHistograms->FillHistogram("ESD_Mother_Eta", MomentumVectorTwoGammaCandidate.Eta());\r
887           fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);                                     \r
888           fHistograms->FillHistogram("ESD_Mother_Phi", SpaceVectorTwoGammaCandidate.Phi());\r
889           fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
890           fHistograms->FillHistogram("ESD_Mother_R", SpaceVectorTwoGammaCandidate.Pt());    // Pt in Space == R!!!\r
891           fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), SpaceVectorTwoGammaCandidate.Pt());\r
892           fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
893           fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,MomentumVectorTwoGammaCandidate.Pt());\r
894         }\r
895       }\r
896       delete twoGammaCandidate;\r
897                         \r
898       //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
899                         \r
900                         \r
901     }\r
902   }\r
903 }\r
904 \r
905 void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
906   // see header file for documentation\r
907         \r
908   vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
909   vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
910   for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
911     AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
912     for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
913       AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
914 \r
915       AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
916                         \r
917       Double_t massBG =0.;\r
918       Double_t widthBG = 0.;\r
919       Double_t chi2BG =10000.;  \r
920       backgroundCandidate->GetMass(massBG,widthBG);\r
921       if(backgroundCandidate->GetNDF()>0){\r
922         chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
923         if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
924                                         \r
925           TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
926           TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
927                                         \r
928           Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
929 \r
930           Double_t rapidity;\r
931           if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
932           else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
933 \r
934                                         \r
935                                         \r
936                                         \r
937           if(openingAngleBG < fMinOpeningAngleGhostCut ) continue;   // minimum opening angle to avoid using ghosttracks\r
938                         \r
939                                         \r
940           fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
941           fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
942           fHistograms->FillHistogram("ESD_Background_Pt",  MomentumVectorbackgroundCandidate.Pt());\r
943           fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
944           fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
945           fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
946           fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
947           fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt());  // Pt in Space == R!!!!\r
948           fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
949           fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
950           fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
951         }\r
952       }\r
953       delete backgroundCandidate;   \r
954       //cout << "nearly at the end of background" << endl;\r
955                         \r
956     }\r
957   }\r
958 }\r
959 \r
960 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
961 {\r
962   // Terminate analysis\r
963   //\r
964   AliDebug(1,"Do nothing in Terminate");\r
965 }\r
966 \r
967 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
968 {\r
969   // Create the output container\r
970   if(fOutputContainer != NULL){\r
971     delete fOutputContainer;\r
972     fOutputContainer = NULL;\r
973   }\r
974   if(fOutputContainer == NULL){\r
975     fOutputContainer = new TList();\r
976   }\r
977         \r
978   //Adding the histograms to the output container\r
979   fHistograms->GetOutputContainer(fOutputContainer);\r
980         \r
981         \r
982   if(fWriteNtuple){\r
983     if(fGammaNtuple == NULL){\r
984       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
985     }\r
986     if(fNeutralMesonNtuple == NULL){\r
987       fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
988     }\r
989     TList * ntupleTList = new TList();\r
990     ntupleTList->SetName("Ntuple");\r
991     ntupleTList->Add((TNtuple*)fGammaNtuple);\r
992     fOutputContainer->Add(ntupleTList);\r
993   }\r
994         \r
995   fOutputContainer->SetName(GetName());\r
996 }\r
997 \r
998 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* daughter0, TParticle* daughter1) const{\r
999   //helper function\r
1000   TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
1001   TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
1002   return v3D0.Angle(v3D1);\r
1003 }\r