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