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