git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
1. Adding sorting of the Points in point array
[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
15e046fa 29#include "TNtuple.h"\r
3c538586 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
32a6d407 49AliAnalysisTaskSE(),\r
3c538586 50 fV0Reader(NULL),\r
51 fStack(NULL),\r
00a6a31a 52 fESDEvent(NULL), \r
3c538586 53 fOutputContainer(NULL),\r
54 fHistograms(NULL),\r
55 fDoMCTruth(kFALSE),\r
56 fMCAllGammas(),\r
57 fMCPi0s(),\r
58 fMCEtas(),\r
59 fMCGammaChic(),\r
60 fKFReconstructedGammas(),\r
32a6d407 61 fIsTrueReconstructedGammas(),\r
15e046fa 62 fElectronv1(),\r
63 fElectronv2(),\r
00a6a31a 64 fCurrentEventPosElectron(),\r
65 fPreviousEventPosElectron(),\r
66 fCurrentEventNegElectron(),\r
67 fPreviousEventNegElectron(),\r
68 fKFReconstructedGammasCut(), \r
69 fPreviousEventTLVNegElectron(),\r
70 fPreviousEventTLVPosElectron(), \r
3c538586 71 fElectronMass(-1),\r
72 fGammaMass(-1),\r
73 fPi0Mass(-1),\r
74 fEtaMass(-1),\r
75 fGammaWidth(-1),\r
76 fPi0Width(-1),\r
77 fEtaWidth(-1),\r
32a6d407 78 fMinOpeningAngleGhostCut(0.),\r
a19c3402 79 fCalculateBackground(kFALSE),\r
80 fWriteNtuple(kFALSE),\r
81 fGammaNtuple(NULL),\r
82 fNeutralMesonNtuple(NULL),\r
83 fTotalNumberOfAddedNtupleEntries(0)\r
3c538586 84{\r
85 // Default constructor\r
86 // Common I/O in slot 0\r
87 DefineInput (0, TChain::Class());\r
88 DefineOutput(0, TTree::Class());\r
32a6d407 89 \r
3c538586 90 // Your private output\r
91 DefineOutput(1, TList::Class());\r
92}\r
93\r
94AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
95 AliAnalysisTaskSE(name),\r
96 fV0Reader(NULL),\r
97 fStack(NULL),\r
00a6a31a 98 fESDEvent(NULL), \r
3c538586 99 fOutputContainer(0x0),\r
100 fHistograms(NULL),\r
101 fDoMCTruth(kFALSE),\r
102 fMCAllGammas(),\r
103 fMCPi0s(),\r
104 fMCEtas(),\r
105 fMCGammaChic(),\r
106 fKFReconstructedGammas(),\r
32a6d407 107 fIsTrueReconstructedGammas(),\r
15e046fa 108 fElectronv1(),\r
109 fElectronv2(),\r
00a6a31a 110 fCurrentEventPosElectron(),\r
111 fPreviousEventPosElectron(),\r
112 fCurrentEventNegElectron(),\r
113 fPreviousEventNegElectron(),\r
114 fKFReconstructedGammasCut(), \r
115 fPreviousEventTLVNegElectron(),\r
116 fPreviousEventTLVPosElectron(),\r
3c538586 117 fElectronMass(-1),\r
118 fGammaMass(-1),\r
119 fPi0Mass(-1),\r
120 fEtaMass(-1),\r
121 fGammaWidth(-1),\r
122 fPi0Width(-1),\r
123 fEtaWidth(-1),\r
32a6d407 124 fMinOpeningAngleGhostCut(0.),\r
a19c3402 125 fCalculateBackground(kFALSE),\r
126 fWriteNtuple(kFALSE),\r
127 fGammaNtuple(NULL),\r
128 fNeutralMesonNtuple(NULL),\r
129 fTotalNumberOfAddedNtupleEntries(0)\r
3c538586 130{\r
131 // Common I/O in slot 0\r
132 DefineInput (0, TChain::Class());\r
133 DefineOutput(0, TTree::Class());\r
32a6d407 134 \r
3c538586 135 // Your private output\r
136 DefineOutput(1, TList::Class());\r
137}\r
138\r
139AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
140{\r
141 // Remove all pointers\r
32a6d407 142 \r
3c538586 143 if(fOutputContainer){\r
144 fOutputContainer->Clear() ; \r
145 delete fOutputContainer ;\r
146 }\r
147 if(fHistograms){\r
148 delete fHistograms;\r
149 }\r
150 if(fV0Reader){\r
151 delete fV0Reader;\r
152 }\r
153}\r
154\r
155\r
156void AliAnalysisTaskGammaConversion::Init()\r
157{\r
158 // Initialization\r
159 AliLog::SetGlobalLogLevel(AliLog::kError);\r
160}\r
161\r
162\r
163void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
164{\r
165 // Execute analysis for current event\r
32a6d407 166 \r
3c538586 167 ConnectInputData("");\r
32a6d407 168 \r
3c538586 169 //clear vectors\r
170 fMCAllGammas.clear();\r
171 fMCPi0s.clear();\r
172 fMCEtas.clear();\r
173 fMCGammaChic.clear();\r
32a6d407 174 \r
3c538586 175 fKFReconstructedGammas.clear();\r
32a6d407 176 fIsTrueReconstructedGammas.clear();\r
15e046fa 177 fElectronv1.clear();\r
178 fElectronv2.clear();\r
00a6a31a 179 fCurrentEventPosElectron.clear();\r
180 fCurrentEventNegElectron.clear(); \r
181 fKFReconstructedGammasCut.clear(); \r
32a6d407 182 \r
3c538586 183 //Clear the data in the v0Reader\r
184 fV0Reader->UpdateEventByEventData();\r
185\r
32a6d407 186 \r
3c538586 187 // Process the MC information\r
188 if(fDoMCTruth){\r
189 ProcessMCData();\r
190 }\r
32a6d407 191 \r
192 //Process the v0 information with no cuts\r
193 ProcessV0sNoCut();\r
194 \r
3c538586 195 // Process the v0 information\r
196 ProcessV0s();\r
32a6d407 197 \r
3c538586 198 //calculate background if flag is set\r
199 if(fCalculateBackground){\r
200 CalculateBackground();\r
201 }\r
32a6d407 202 \r
3c538586 203 // Process reconstructed gammas\r
204 ProcessGammasForNeutralMesonAnalysis();\r
f78c0011 205\r
206 CheckV0Efficiency();\r
3c538586 207 \r
00a6a31a 208 //Process reconstructed gammas electrons for Chi_c Analysis\r
209 ProcessGammaElectronsForChicAnalysis();\r
210\r
32a6d407 211 PostData(1, fOutputContainer);\r
212 \r
3c538586 213}\r
214\r
215void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
216 // see header file for documentation\r
32a6d407 217 \r
3c538586 218 if(fV0Reader == NULL){\r
219 // Write warning here cuts and so on are default if this ever happens\r
220 }\r
221 fV0Reader->Initialize();\r
222}\r
223\r
32a6d407 224\r
225\r
3c538586 226void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
227 // see header file for documentation\r
32a6d407 228 \r
3c538586 229 fStack = fV0Reader->GetMCStack();\r
230\r
32a6d407 231 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
232 return; // aborts if the primary vertex does not have contributors.\r
233 }\r
234\r
3c538586 235 for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
236 TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
32a6d407 237 \r
3c538586 238 if (!particle) {\r
239 //print warning here\r
240 continue;\r
241 }\r
3c538586 242\r
00a6a31a 243 ///////////////////////Begin Chic Analysis/////////////////////////////\r
244\r
245\r
246 if(particle->GetPdgCode() == 443){//Is JPsi\r
247\r
248 if(particle->GetNDaughters()==2){\r
249 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&\r
250 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){\r
251 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
252 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
253 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)\r
254 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance\r
255\r
256 if( TMath::Abs(daug0->Eta()) < 0.9){\r
257 if(daug0->GetPdgCode() == -11)\r
258 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
259 else\r
260 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
261\r
262 }\r
263 if(TMath::Abs(daug1->Eta()) < 0.9){\r
264 if(daug1->GetPdgCode() == -11)\r
265 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
266 else\r
267 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
268 }\r
269 }\r
270 }\r
271 }\r
272 // const int CHI_C0 = 10441;\r
273 // const int CHI_C1 = 20443;\r
274 // const int CHI_C2 = 445\r
275 if(particle->GetPdgCode() == 22){//gamma from JPsi\r
276 if(particle->GetMother(0) > -1){\r
277 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||\r
278 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||\r
279 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){\r
280 if(TMath::Abs(particle->Eta()) < 1.2)\r
281 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance\r
282 }\r
283 }\r
284 }\r
285 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){\r
286 if( particle->GetNDaughters() == 2){\r
287 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
288 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
289\r
290 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){\r
291 if( daug0->GetPdgCode() == 443){\r
292 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());\r
293 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());\r
294 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
295 fHistograms->FillTable("Table_Electrons",18);\r
296\r
297 }//if\r
298 else if (daug1->GetPdgCode() == 443){\r
299 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());\r
300 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());\r
301 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
302 fHistograms->FillTable("Table_Electrons",18);\r
303 }//else if\r
304 }//gamma o Jpsi\r
305 }//GetNDaughters\r
306 }\r
307\r
308\r
309 /////////////////////End Chic Analysis////////////////////////////\r
310\r
311\r
32a6d407 312 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;\r
313 \r
314 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point\r
315 \r
3c538586 316 Double_t tmpPhi=particle->Phi();\r
32a6d407 317 \r
3c538586 318 if(particle->Phi()> TMath::Pi()){\r
319 tmpPhi = particle->Phi()-(2*TMath::Pi());\r
320 }\r
32a6d407 321 \r
322 Double_t rapidity;\r
323 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
324 rapidity=0;\r
325 }\r
326 else{\r
327 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
328 } \r
329 \r
3c538586 330 //process the gammas\r
32a6d407 331 if (particle->GetPdgCode() == 22){\r
332 \r
333 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
334 continue; // no photon as mothers!\r
335 }\r
3c538586 336\r
32a6d407 337 if(particle->GetMother(0) >= fStack->GetNprimary()){\r
338 continue; // the gamma has a mother, and it is not a primary particle\r
339 }\r
3c538586 340\r
32a6d407 341 fMCAllGammas.push_back(particle);\r
342 \r
343 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
344 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
345 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
346 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
347 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
348 \r
349 \r
350 if(particle->GetMother(0) < 0){ // direct gamma\r
351 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
352 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
353 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
354 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
355 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity); \r
356 }\r
357 \r
358 \r
359 // looking for conversion (electron + positron from pairbuilding (= 5) )\r
360 TParticle* ePos = NULL;\r
361 TParticle* eNeg = NULL;\r
362 \r
363 if(particle->GetNDaughters() >= 2){\r
364 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
365 TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
366 if(tmpDaughter->GetUniqueID() == 5){\r
367 if(tmpDaughter->GetPdgCode() == 11){\r
368 eNeg = tmpDaughter;\r
369 }\r
370 else if(tmpDaughter->GetPdgCode() == -11){\r
371 ePos = tmpDaughter;\r
3c538586 372 }\r
32a6d407 373 }\r
374 }\r
375 }\r
376 \r
377 \r
378 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
379 continue;\r
380 }\r
381 \r
382 \r
383 Double_t ePosPhi = ePos->Phi();\r
384 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
385 \r
386 Double_t eNegPhi = eNeg->Phi();\r
387 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
388 \r
389 \r
390 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
391 continue; // no reconstruction below the Pt cut\r
392 }\r
393 \r
394 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
395 continue;\r
396 } \r
397 \r
398 if(ePos->R()>fV0Reader->GetMaxRCut()){\r
399 continue; // cuts on distance from collision point\r
400 }\r
401 \r
402 \r
403 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){\r
404 continue; // line cut to exclude regions where we do not reconstruct\r
405 } \r
406 \r
407 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
408 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
409 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
410 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
411 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
412 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
413 \r
414 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
415 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
416 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
417 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
418 \r
419 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
420 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
421 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
422 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
423 \r
424 \r
425 \r
426 //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
427 \r
428 // begin Mapping \r
429 Int_t rBin = fHistograms->GetRBin(ePos->R());\r
430 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());\r
431 \r
432 TString nameMCMappingPhiR="";\r
433 nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
434 fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
435 \r
436 TString nameMCMappingPhi="";\r
437 nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
438 fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
439 \r
440 TString nameMCMappingR="";\r
441 nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
442 fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
443 \r
444 TString nameMCMappingPhiInR="";\r
445 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
446 fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
447 //end mapping\r
448 \r
449 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
450 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
451 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
452 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
453 \r
454 //cout << "mapping is done" << endl;\r
455 \r
456 \r
457 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
458 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
459 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
460 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
461 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
462 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
463 \r
464 } // end direct gamma\r
465 else{ // mother exits \r
3c538586 466 if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
467 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
468 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2\r
469 ){ \r
470 fMCGammaChic.push_back(particle);\r
471 }\r
32a6d407 472 } // end if mother exits\r
473 } // end if particle is a photon\r
474 \r
475 if(particle->GetNDaughters() == 2){\r
476 \r
477 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
478 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
479 \r
480 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
481 \r
482 \r
483 \r
484 // check for conversions now -> have to pass eta and line cut!\r
485 Bool_t daughter0Electron = kFALSE;\r
486 Bool_t daughter0Positron = kFALSE;\r
487 Bool_t daughter1Electron = kFALSE;\r
488 Bool_t daughter1Positron = kFALSE;\r
489 \r
490 \r
491 \r
492 if(daughter0->GetNDaughters() >= 2){\r
493 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
494 TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
495 if(tmpDaughter->GetUniqueID() == 5){\r
496 if(tmpDaughter->GetPdgCode() == 11){\r
497 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
498 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
499 daughter0Electron = kTRUE;\r
500 }\r
501 \r
502 }\r
503 }\r
504 else if(tmpDaughter->GetPdgCode() == -11){\r
505 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
506 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
507 daughter0Positron = kTRUE;\r
508 }\r
509 \r
510 }\r
511 \r
512 }\r
3c538586 513 }\r
514 }\r
32a6d407 515 }\r
516 \r
517 \r
518 \r
519 if(daughter1->GetNDaughters() >= 2){\r
520 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
521 TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
522 if(tmpDaughter->GetUniqueID() == 5){\r
523 if(tmpDaughter->GetPdgCode() == 11){\r
524 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
525 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
526 daughter1Electron = kTRUE;\r
527 }\r
528 \r
529 }\r
530 }\r
531 else if(tmpDaughter->GetPdgCode() == -11){\r
532 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
533 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
534 daughter1Positron = kTRUE;\r
535 }\r
536 \r
537 }\r
538 \r
539 }\r
3c538586 540 }\r
541 }\r
542 }\r
32a6d407 543 \r
544 \r
545 \r
546 \r
547 if(particle->GetPdgCode()==111){ //Pi0\r
548 if( iTracks >= fStack->GetNprimary()){\r
549 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
550 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
3c538586 551 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
32a6d407 552 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
553 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
554 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
555 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
556 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
557 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
558 \r
559 if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
560 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
561 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
562 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
563 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
564 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
565 }\r
3c538586 566 }\r
32a6d407 567 }\r
568 else{\r
569 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta()); \r
570 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
571 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
572 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
573 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
574 fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
575 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
576 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
577 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
578 \r
579 if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
580 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
581 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
582 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
583 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
584 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
585 }\r
3c538586 586 }\r
587 }\r
32a6d407 588 }\r
589 \r
590 if(particle->GetPdgCode()==221){ //Eta\r
591 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
592 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
593 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
594 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
595 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
596 fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
597 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
598 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
599 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
600 \r
601 if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
602 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
603 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
604 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
605 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
606 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
607 }\r
608 \r
3c538586 609 }\r
32a6d407 610 \r
611 }\r
612 \r
613 // all motherparticles with 2 gammas as daughters\r
614 fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
615 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
616 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
617 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
618 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
619 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
620 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
621 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
622 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
623 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
624 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt()); \r
625 if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
626 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
627 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
628 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt()); \r
629 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
630 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
631 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
632 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt()); \r
633\r
634 }\r
635 \r
636 \r
3c538586 637 }\r
32a6d407 638 \r
639 //cout << "mother histos are filled" << endl;\r
640 \r
641 } // end if(particle->GetNDaughters() == 2)\r
642 \r
3c538586 643 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
32a6d407 644 \r
645 //cout << "right before the end of processMCdata" << endl;\r
646 \r
3c538586 647} // end ProcessMCData\r
648\r
a19c3402 649\r
32a6d407 650\r
651void AliAnalysisTaskGammaConversion::FillNtuple(){\r
15e046fa 652 //Fills the ntuple with the different values\r
653\r
a19c3402 654 if(fGammaNtuple == NULL){\r
655 return;\r
656 }\r
657 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
658 for(Int_t i=0;i<numberOfV0s;i++){\r
659 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
660 AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
661 Double_t negPID=0;\r
662 Double_t posPID=0;\r
663 fV0Reader->GetPIDProbability(negPID,posPID);\r
664 values[0]=cV0->GetOnFlyStatus();\r
665 values[1]=fV0Reader->CheckForPrimaryVertex();\r
666 values[2]=negPID;\r
667 values[3]=posPID;\r
668 values[4]=fV0Reader->GetX();\r
669 values[5]=fV0Reader->GetY();\r
670 values[6]=fV0Reader->GetZ();\r
671 values[7]=fV0Reader->GetXYRadius();\r
672 values[8]=fV0Reader->GetMotherCandidateNDF();\r
673 values[9]=fV0Reader->GetMotherCandidateChi2();\r
674 values[10]=fV0Reader->GetMotherCandidateEnergy();\r
675 values[11]=fV0Reader->GetMotherCandidateEta();\r
676 values[12]=fV0Reader->GetMotherCandidatePt();\r
677 values[13]=fV0Reader->GetMotherCandidateMass();\r
678 values[14]=fV0Reader->GetMotherCandidateWidth();\r
679 // 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
680 values[16]=fV0Reader->GetOpeningAngle();\r
681 values[17]=fV0Reader->GetNegativeTrackEnergy();\r
682 values[18]=fV0Reader->GetNegativeTrackPt();\r
683 values[19]=fV0Reader->GetNegativeTrackEta();\r
684 values[20]=fV0Reader->GetNegativeTrackPhi();\r
685 values[21]=fV0Reader->GetPositiveTrackEnergy();\r
686 values[22]=fV0Reader->GetPositiveTrackPt();\r
687 values[23]=fV0Reader->GetPositiveTrackEta();\r
688 values[24]=fV0Reader->GetPositiveTrackPhi();\r
689 values[25]=fV0Reader->HasSameMCMother();\r
690 if(values[25] != 0){\r
691 values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
692 values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
693 }\r
694 fTotalNumberOfAddedNtupleEntries++;\r
695 fGammaNtuple->Fill(values);\r
696 }\r
697 fV0Reader->ResetV0IndexNumber();\r
32a6d407 698 \r
699}\r
700\r
701void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
15e046fa 702 // Process all the V0's without applying any cuts to it\r
32a6d407 703\r
704 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
705 for(Int_t i=0;i<numberOfV0s;i++){\r
706 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
707\r
708 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
709 return;\r
710 }\r
711 \r
712 if(fDoMCTruth){\r
713 \r
714 if(fV0Reader->HasSameMCMother() == kFALSE){\r
715 continue;\r
716 }\r
717 \r
718 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
719 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
720\r
721 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
722 continue;\r
723 }\r
724 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
725 continue;\r
726 }\r
727 \r
728 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
729 \r
730 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
731 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
732 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta()); \r
733 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
734 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
735 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
736 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
737 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
738 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
739 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
740 \r
741 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
742 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
743 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
744 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
745 \r
746 /*\r
747 ESD_NoCutConvGamma_Pt\r
748 ESD_NoCutConvGamma_Energy\r
749 ESD_NoCutConvGamma_Eta\r
750 ESD_NoCutConvGamma_Phi\r
751 ESD_NoCutConvGamma_Mass\r
752 ESD_NoCutConvGamma_Width\r
753 ESD_NoCutConvGamma_Chi2\r
754 ESD_NoCutConvGamma_NDF\r
755 ESD_NoCutConvGamma_PtvsEta\r
756 ESD_NoCutConversion_XY\r
757 ESD_NoCutConversion_R\r
758 ESD_NoCutConversion_ZR\r
759 ESD_NoCutConversion_OpeningAngle\r
760 */\r
761 }\r
762 }\r
763 }\r
764 fV0Reader->ResetV0IndexNumber();\r
a19c3402 765}\r
766\r
3c538586 767void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
768 // see header file for documentation\r
32a6d407 769 \r
a19c3402 770 if(fWriteNtuple == kTRUE){\r
771 FillNtuple();\r
772 }\r
32a6d407 773 \r
3c538586 774 Int_t nSurvivingV0s=0;\r
775 while(fV0Reader->NextV0()){\r
776 nSurvivingV0s++;\r
32a6d407 777 \r
778 \r
3c538586 779 //-------------------------- filling v0 information -------------------------------------\r
32a6d407 780 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
781 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
782 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
783 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle()); \r
784 \r
3c538586 785 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
786 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
787 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
788 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
32a6d407 789 \r
3c538586 790 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
791 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
792 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
793 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
32a6d407 794 \r
795 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
796 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
797 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
798 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
799 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
800 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
801 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
802 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
803 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
804 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
805 \r
806 \r
3c538586 807 // begin mapping\r
808 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
809 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
810 Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
32a6d407 811 \r
3c538586 812 TString nameESDMappingPhiR="";\r
32a6d407 813 nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
3c538586 814 fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
32a6d407 815 \r
3c538586 816 TString nameESDMappingPhi="";\r
32a6d407 817 nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
3c538586 818 fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
32a6d407 819 \r
3c538586 820 TString nameESDMappingR="";\r
32a6d407 821 nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
3c538586 822 fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta); \r
32a6d407 823 \r
3c538586 824 TString nameESDMappingPhiInR="";\r
32a6d407 825 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
3c538586 826 fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
827 // end mapping\r
32a6d407 828 \r
3c538586 829 fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
15e046fa 830 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
831 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
3c538586 832\r
32a6d407 833 \r
3c538586 834 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
835 if(fDoMCTruth){\r
32a6d407 836 \r
3c538586 837 if(fV0Reader->HasSameMCMother() == kFALSE){\r
32a6d407 838 fIsTrueReconstructedGammas.push_back(kFALSE);\r
3c538586 839 continue;\r
840 }\r
32a6d407 841 \r
842 \r
3c538586 843 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
844 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
845\r
32a6d407 846 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
847 fIsTrueReconstructedGammas.push_back(kFALSE);\r
3c538586 848 continue;\r
849 }\r
32a6d407 850 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
851 fIsTrueReconstructedGammas.push_back(kFALSE);\r
852 continue;\r
853 }\r
854 \r
3c538586 855 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
32a6d407 856 fIsTrueReconstructedGammas.push_back(kTRUE);\r
857 \r
858 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
859 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
860 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta()); \r
861 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
862 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
863 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
864 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
865 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
866 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
867 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
868 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
869 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
870 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
871 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
872 \r
873 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
874 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
875 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
876 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
877\r
878 \r
879 /*\r
880 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
881 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
882 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
883 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
884 */\r
885\r
886\r
3c538586 887\r
888 //resolution\r
889 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();\r
890 Double_t esdpt = fV0Reader->GetMotherCandidatePt();\r
891 Double_t resdPt = 0;\r
892 if(mcpt != 0){\r
893 resdPt = ((esdpt - mcpt)/mcpt)*100;\r
894 }\r
32a6d407 895 \r
3c538586 896 fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
897 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
898 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
32a6d407 899 \r
3c538586 900 Double_t resdZ = 0;\r
901 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
902 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
903 }\r
32a6d407 904 \r
3c538586 905 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
906 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
907 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
32a6d407 908 \r
3c538586 909 Double_t resdR = 0;\r
910 if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
911 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
912 }\r
32a6d407 913\r
3c538586 914 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
915 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
916 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
917 fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
32a6d407 918 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
919 else{\r
920 fIsTrueReconstructedGammas.push_back(kFALSE);\r
3c538586 921 }\r
32a6d407 922 }//if(fDoMCTruth)\r
923 }//while(fV0Reader->NextV0)\r
924 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
925 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
926 \r
927 //cout << "nearly at the end of doMCTruth" << endl;\r
928 \r
3c538586 929}\r
930\r
931void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
932 // see header file for documentation\r
32a6d407 933 \r
3c538586 934 for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
935 for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
32a6d407 936 \r
3c538586 937 AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
938 AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
3c538586 939 \r
15e046fa 940 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
32a6d407 941 continue;\r
942 }\r
15e046fa 943 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
32a6d407 944 continue;\r
945 }\r
946\r
947 /*\r
948 if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
949 continue;\r
950 }\r
951 */\r
3c538586 952\r
32a6d407 953 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
954 \r
955 Double_t massTwoGammaCandidate = 0.;\r
3c538586 956 Double_t widthTwoGammaCandidate = 0.;\r
957 Double_t chi2TwoGammaCandidate =10000.; \r
958 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
959 if(twoGammaCandidate->GetNDF()>0){\r
960 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
32a6d407 961 \r
962 if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){ \r
963 \r
15e046fa 964 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
965 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
32a6d407 966 \r
967 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1); \r
968 Double_t rapidity;\r
969 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
970 rapidity=0;\r
971 }\r
972 else{\r
973 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
974 }\r
975 \r
976 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue; // minimum opening angle to avoid using ghosttracks\r
977\r
978 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
979 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
15e046fa 980 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());\r
981 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());\r
32a6d407 982 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity); \r
15e046fa 983 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());\r
32a6d407 984 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
15e046fa 985 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!\r
986 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());\r
32a6d407 987 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
15e046fa 988 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());\r
f78c0011 989 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);\r
3c538586 990 }\r
991 }\r
992 delete twoGammaCandidate;\r
32a6d407 993 \r
994 //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
995 \r
996 \r
3c538586 997 }\r
998 }\r
999}\r
1000\r
1001void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
1002 // see header file for documentation\r
32a6d407 1003 \r
3c538586 1004 vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
1005 vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
1006 for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
1007 AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
1008 for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
1009 AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
1010\r
1011 AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
32a6d407 1012 \r
3c538586 1013 Double_t massBG =0.;\r
1014 Double_t widthBG = 0.;\r
1015 Double_t chi2BG =10000.; \r
1016 backgroundCandidate->GetMass(massBG,widthBG);\r
1017 if(backgroundCandidate->GetNDF()>0){\r
1018 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
1019 if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
32a6d407 1020 \r
1021 TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
1022 TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
1023 \r
3c538586 1024 Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
1025\r
32a6d407 1026 Double_t rapidity;\r
1027 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
1028 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
3c538586 1029\r
32a6d407 1030 \r
1031 \r
1032 \r
1033 if(openingAngleBG < fMinOpeningAngleGhostCut ) continue; // minimum opening angle to avoid using ghosttracks\r
1034 \r
1035 \r
3c538586 1036 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
1037 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
32a6d407 1038 fHistograms->FillHistogram("ESD_Background_Pt", MomentumVectorbackgroundCandidate.Pt());\r
1039 fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
1040 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
1041 fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
3c538586 1042 fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
32a6d407 1043 fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!\r
1044 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
1045 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
1046 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
f78c0011 1047 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);\r
3c538586 1048 }\r
1049 }\r
1050 delete backgroundCandidate; \r
32a6d407 1051 //cout << "nearly at the end of background" << endl;\r
1052 \r
3c538586 1053 }\r
1054 }\r
1055}\r
1056\r
1057void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
1058{\r
1059 // Terminate analysis\r
1060 //\r
1061 AliDebug(1,"Do nothing in Terminate");\r
1062}\r
1063\r
1064void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
1065{\r
1066 // Create the output container\r
1067 if(fOutputContainer != NULL){\r
1068 delete fOutputContainer;\r
1069 fOutputContainer = NULL;\r
1070 }\r
1071 if(fOutputContainer == NULL){\r
1072 fOutputContainer = new TList();\r
1073 }\r
32a6d407 1074 \r
a19c3402 1075 //Adding the histograms to the output container\r
3c538586 1076 fHistograms->GetOutputContainer(fOutputContainer);\r
32a6d407 1077 \r
1078 \r
a19c3402 1079 if(fWriteNtuple){\r
1080 if(fGammaNtuple == NULL){\r
1081 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
1082 }\r
1083 if(fNeutralMesonNtuple == NULL){\r
1084 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
1085 }\r
1086 TList * ntupleTList = new TList();\r
1087 ntupleTList->SetName("Ntuple");\r
1088 ntupleTList->Add((TNtuple*)fGammaNtuple);\r
1089 fOutputContainer->Add(ntupleTList);\r
1090 }\r
32a6d407 1091 \r
3c538586 1092 fOutputContainer->SetName(GetName());\r
1093}\r
1094\r
15e046fa 1095Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{\r
3c538586 1096 //helper function\r
1097 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
1098 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
1099 return v3D0.Angle(v3D1);\r
1100}\r
f78c0011 1101\r
1102void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){\r
613bcce9 1103 // see header file for documentation\r
f78c0011 1104\r
1105 vector<Int_t> indexOfGammaParticle;\r
1106\r
1107 fStack = fV0Reader->GetMCStack();\r
1108\r
1109 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
1110 return; // aborts if the primary vertex does not have contributors.\r
1111 }\r
1112\r
1113 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {\r
1114 TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
1115 if(particle->GetPdgCode()==22){ //Gamma\r
1116 if(particle->GetNDaughters() >= 2){\r
1117 TParticle* electron=NULL;\r
1118 TParticle* positron=NULL; \r
1119 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
1120 TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
1121 if(tmpDaughter->GetUniqueID() == 5){\r
1122 if(tmpDaughter->GetPdgCode() == 11){\r
1123 electron = tmpDaughter;\r
1124 }\r
1125 else if(tmpDaughter->GetPdgCode() == -11){\r
1126 positron = tmpDaughter;\r
1127 }\r
1128 }\r
1129 }\r
1130 if(electron!=NULL && positron!=0){\r
1131 if(electron->R()<160){\r
1132 indexOfGammaParticle.push_back(iTracks);\r
1133 }\r
1134 }\r
1135 }\r
1136 }\r
1137 }\r
1138\r
1139 Int_t nFoundGammas=0;\r
1140 Int_t nNotFoundGammas=0;\r
1141\r
1142 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
1143 for(Int_t i=0;i<numberOfV0s;i++){\r
1144 fV0Reader->GetV0(i);\r
1145 \r
1146 if(fV0Reader->HasSameMCMother() == kFALSE){\r
1147 continue;\r
1148 }\r
1149 \r
1150 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
1151 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
1152 \r
1153 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
1154 continue;\r
1155 }\r
1156 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
1157 continue;\r
1158 }\r
1159 \r
1160 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
1161 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();\r
1162 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){\r
1163 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){\r
1164 nFoundGammas++;\r
1165 }\r
1166 else{\r
1167 nNotFoundGammas++;\r
1168 }\r
1169 }\r
1170 }\r
1171 }\r
1172 // cout<<"Found: "<<nFoundGammas<<" of: "<<indexOfGammaParticle.size()<<endl;\r
1173}\r
00a6a31a 1174\r
1175\r
1176void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){\r
613bcce9 1177 // see header file for documantation\r
00a6a31a 1178\r
1179 fESDEvent = fV0Reader->GetESDEvent();\r
1180\r
1181\r
613bcce9 1182 vector <AliESDtrack*> vESDeNegTemp(0);\r
1183 vector <AliESDtrack*> vESDePosTemp(0);\r
1184 vector <AliESDtrack*> vESDxNegTemp(0);\r
1185 vector <AliESDtrack*> vESDxPosTemp(0);\r
1186 vector <AliESDtrack*> vESDeNegNoJPsi(0);\r
1187 vector <AliESDtrack*> vESDePosNoJPsi(0); \r
00a6a31a 1188\r
1189\r
1190\r
1191 fHistograms->FillTable("Table_Electrons",0);//Count number of Events\r
1192\r
1193 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){\r
1194 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);\r
1195\r
1196 if(!curTrack){\r
1197 //print warning here\r
1198 continue;\r
1199 }\r
1200\r
1201 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;\r
1202 double r[3];curTrack->GetConstrainedXYZ(r);\r
1203\r
1204 TVector3 rXYZ(r);\r
1205\r
1206 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks\r
1207\r
1208 Bool_t flagKink = kTRUE;\r
1209 Bool_t flagTPCrefit = kTRUE;\r
1210 Bool_t flagTRDrefit = kTRUE;\r
1211 Bool_t flagITSrefit = kTRUE;\r
1212 Bool_t flagTRDout = kTRUE;\r
1213 Bool_t flagVertex = kTRUE;\r
1214\r
1215\r
1216 //Cuts ---------------------------------------------------------------\r
1217\r
1218 if(curTrack->GetKinkIndex(0) > 0){\r
1219 fHistograms->FillHistogram("Table_Electrons",5);//Count kink\r
1220 flagKink = kFALSE;\r
1221 }\r
1222\r
1223 ULong_t trkStatus = curTrack->GetStatus();\r
1224\r
1225 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);\r
1226\r
1227 if(!tpcRefit){\r
1228 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit\r
1229 flagTPCrefit = kFALSE;\r
1230 }\r
1231\r
1232 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);\r
1233 if(!itsRefit){\r
1234 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit\r
1235 flagITSrefit = kFALSE;\r
1236 }\r
1237\r
1238 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);\r
1239\r
1240 if(!trdRefit){\r
1241 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit\r
1242 flagTRDrefit = kFALSE;\r
1243 }\r
1244\r
1245 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);\r
1246\r
1247 if(!trdOut) {\r
1248 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout\r
1249 flagTRDout = kFALSE;\r
1250 }\r
1251\r
1252 double nsigmaToVxt = GetSigmaToVertex(curTrack);\r
1253\r
1254 if(nsigmaToVxt > 3){\r
1255 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3\r
1256 flagVertex = kFALSE;\r
1257 }\r
1258\r
1259 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;\r
1260 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts\r
1261\r
1262\r
1263 Stat_t pid, weight;\r
1264 GetPID(curTrack, pid, weight);\r
1265\r
1266 if(pid!=0){\r
1267 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0\r
1268 }\r
1269\r
1270 if(pid == 0){\r
1271 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0\r
1272 }\r
1273\r
1274\r
1275\r
1276\r
1277 Int_t labelMC = TMath::Abs(curTrack->GetLabel());\r
1278 TParticle* curParticle = fStack->Particle(labelMC);\r
1279\r
1280\r
1281\r
1282\r
1283 TLorentzVector curElec;\r
1284 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);\r
1285\r
1286\r
1287\r
1288\r
1289 if(curTrack->GetSign() > 0){\r
1290\r
613bcce9 1291 vESDxPosTemp.push_back(curTrack);\r
00a6a31a 1292\r
1293 if( pid == 0){\r
1294\r
1295 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1296 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());\r
1297 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1298 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1299 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
613bcce9 1300 vESDePosTemp.push_back(curTrack);\r
00a6a31a 1301\r
1302\r
1303\r
1304 }\r
1305\r
1306 }\r
1307 else {\r
613bcce9 1308 vESDxNegTemp.push_back(curTrack);\r
00a6a31a 1309\r
1310 if( pid == 0){\r
1311\r
1312 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
1313 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());\r
1314 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
1315 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
1316 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
613bcce9 1317 vESDeNegTemp.push_back(curTrack);\r
00a6a31a 1318\r
1319\r
1320\r
1321\r
1322 }\r
1323 }\r
1324\r
1325 }\r
1326\r
1327\r
1328 Bool_t ePosJPsi = kFALSE;\r
1329 Bool_t eNegJPsi = kFALSE; \r
1330 Bool_t ePosPi0 = kFALSE;\r
1331 Bool_t eNegPi0 = kFALSE;\r
1332 \r
1333 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;\r
1334 \r
613bcce9 1335 for(UInt_t iNeg=0; iNeg < vESDeNegTemp.size(); iNeg++){\r
1336 if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)\r
1337 if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){\r
1338 Int_t labelMother = fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);\r
00a6a31a 1339 TParticle* partMother = fStack ->Particle(labelMother);\r
1340 if (partMother->GetPdgCode() == 111){\r
1341 ieNegPi0 = iNeg;\r
1342 eNegPi0 = kTRUE;\r
1343 }\r
1344 if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1345 fHistograms->FillTable("Table_Electrons",14);\r
1346 ieNegJPsi = iNeg;\r
1347 eNegJPsi = kTRUE;\r
1348 }\r
1349 else{ \r
613bcce9 1350 vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);\r
00a6a31a 1351 // cout<<"ESD No Positivo JPsi "<<endl;\r
1352 }\r
1353\r
1354 }\r
1355 } \r
1356\r
613bcce9 1357 for(UInt_t iPos=0; iPos < vESDePosTemp.size(); iPos++){\r
1358 if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)\r
1359 if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){\r
1360 Int_t labelMother = fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0);\r
00a6a31a 1361 TParticle* partMother = fStack ->Particle(labelMother);\r
1362 if (partMother->GetPdgCode() == 111){\r
1363 iePosPi0 = iPos;\r
1364 ePosPi0 = kTRUE;\r
1365 }\r
1366 if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
1367 fHistograms->FillTable("Table_Electrons",15);\r
1368 iePosJPsi = iPos;\r
1369 ePosJPsi = kTRUE;\r
1370 }\r
1371 else{\r
613bcce9 1372 vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);\r
00a6a31a 1373 // cout<<"ESD No Negativo JPsi "<<endl;\r
1374 }\r
1375\r
1376 }\r
1377 }\r
1378 \r
1379 if( eNegJPsi && ePosJPsi ){\r
1380 TVector3 tempeNegV,tempePosV;\r
613bcce9 1381 tempeNegV.SetXYZ(vESDeNegTemp[ieNegJPsi]->Px(),vESDeNegTemp[ieNegJPsi]->Py(),vESDeNegTemp[ieNegJPsi]->Pz()); \r
1382 tempePosV.SetXYZ(vESDePosTemp[iePosJPsi]->Px(),vESDePosTemp[iePosJPsi]->Py(),vESDePosTemp[iePosJPsi]->Pz());\r
00a6a31a 1383 fHistograms->FillTable("Table_Electrons",16);\r
1384 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV)); \r
613bcce9 1385 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegJPsi]->GetLabel())),\r
1386 fStack->Particle(TMath::Abs(vESDePosTemp[iePosJPsi]->GetLabel())))); \r
00a6a31a 1387 }\r
1388 \r
1389 if( eNegPi0 && ePosPi0 ){\r
1390 TVector3 tempeNegV,tempePosV;\r
613bcce9 1391 tempeNegV.SetXYZ(vESDeNegTemp[ieNegPi0]->Px(),vESDeNegTemp[ieNegPi0]->Py(),vESDeNegTemp[ieNegPi0]->Pz());\r
1392 tempePosV.SetXYZ(vESDePosTemp[iePosPi0]->Px(),vESDePosTemp[iePosPi0]->Py(),vESDePosTemp[iePosPi0]->Pz());\r
00a6a31a 1393 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));\r
613bcce9 1394 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegPi0]->GetLabel())),\r
1395 fStack->Particle(TMath::Abs(vESDePosTemp[iePosPi0]->GetLabel())))); \r
00a6a31a 1396 }\r
1397 \r
1398\r
613bcce9 1399 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));\r
00a6a31a 1400\r
613bcce9 1401 CleanWithAngleCuts(vESDeNegTemp,vESDePosTemp,fKFReconstructedGammas);\r
00a6a31a 1402 \r
613bcce9 1403 vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);\r
1404 vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);\r
00a6a31a 1405\r
1406\r
613bcce9 1407 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);\r
00a6a31a 1408\r
1409 \r
1410\r
1411\r
1412 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
1413\r
1414\r
613bcce9 1415 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);\r
1416 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));\r
00a6a31a 1417\r
1418 \r
1419\r
1420 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",\r
613bcce9 1421 fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
00a6a31a 1422\r
1423 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",\r
613bcce9 1424 fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
00a6a31a 1425\r
1426 //BackGround\r
1427\r
1428 //Like Sign e+e-\r
613bcce9 1429 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);\r
1430 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);\r
1431 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);\r
1432 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);\r
00a6a31a 1433\r
1434 // Like Sign e+e- no JPsi\r
613bcce9 1435 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));\r
1436 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));\r
00a6a31a 1437\r
1438 //Mixed Event\r
1439\r
1440 if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){\r
1441 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",\r
1442 fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);\r
613bcce9 1443 fPreviousEventTLVNegElectron = vCurrentTLVeNeg;\r
1444 fPreviousEventTLVPosElectron = vCurrentTLVePos;\r
00a6a31a 1445\r
1446 }\r
1447\r
1448 /*\r
1449 //Photons P\r
1450 Double_t vtx[3];\r
1451 vtx[0]=0;vtx[1]=0;vtx[2]=0;\r
1452 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){\r
1453\r
1454 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;\r
1455\r
1456\r
1457\r
1458 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();\r
1459 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;\r
1460\r
1461 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;\r
1462\r
1463 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )\r
1464\r
1465 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());\r
1466\r
1467\r
1468 }\r
1469\r
1470\r
1471 */\r
1472\r
1473\r
1474}\r
1475\r
613bcce9 1476void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){\r
1477 //see header file for documentation\r
1478 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){\r
1479 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
1480 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));\r
00a6a31a 1481 }\r
1482 }\r
1483}\r
1484void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){\r
613bcce9 1485 //see header file for documentation\r
00a6a31a 1486 for( UInt_t n=0; n < eNeg.size(); n++){\r
1487\r
1488 TLorentzVector en = eNeg.at(n);\r
1489 for (UInt_t p=0; p < ePos.size(); p++){\r
1490 TLorentzVector ep = ePos.at(p);\r
1491 TLorentzVector np = ep + en;\r
1492 fHistograms->FillHistogram(histoName.Data(),np.M());\r
1493 }\r
1494 }\r
1495\r
1496}\r
1497\r
1498void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,\r
613bcce9 1499 vector <TLorentzVector> tlVeNeg,vector<TLorentzVector> tlVePos)\r
00a6a31a 1500{\r
613bcce9 1501 //see header file for documentation\r
00a6a31a 1502\r
613bcce9 1503 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++ ){\r
00a6a31a 1504\r
613bcce9 1505 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
00a6a31a 1506\r
613bcce9 1507 TLorentzVector xy = tlVePos[iPos] + tlVeNeg[iNeg];\r
00a6a31a 1508\r
1509 for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){\r
1510\r
613bcce9 1511 AliKFParticle * gammaCandidate = &fKFGammas[iGam];\r
00a6a31a 1512 TLorentzVector g;\r
1513\r
613bcce9 1514 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);\r
00a6a31a 1515 TLorentzVector xyg = xy + g;\r
1516 fHistograms->FillHistogram(histoMass.Data(),xyg.M());\r
1517 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));\r
1518 }\r
1519 }\r
1520 }\r
1521\r
1522}\r
1523void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)\r
1524{\r
613bcce9 1525 // see header file for documentation\r
00a6a31a 1526 for(UInt_t i=0; i < e.size(); i++)\r
1527 {\r
1528 for (UInt_t j=i+1; j < e.size(); j++)\r
1529 {\r
1530 TLorentzVector ee = e[i] + e[j];\r
1531\r
1532 fHistograms->FillHistogram(hBg.Data(),ee.M());\r
1533 }\r
1534 }\r
1535}\r
1536\r
1537\r
613bcce9 1538void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> negativeElectrons,\r
1539 vector <AliESDtrack*> positiveElectrons, vector <AliKFParticle> gammas){\r
1540 // see header file for documentation\r
00a6a31a 1541\r
613bcce9 1542 UInt_t sizeN = negativeElectrons.size();\r
1543 UInt_t sizeP = positiveElectrons.size();\r
1544 UInt_t sizeG = gammas.size();\r
00a6a31a 1545\r
1546\r
1547\r
613bcce9 1548 vector <Bool_t> xNegBand(sizeN);\r
1549 vector <Bool_t> xPosBand(sizeP);\r
1550 vector <Bool_t> gammaBand(sizeG);\r
00a6a31a 1551\r
1552\r
613bcce9 1553 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;\r
1554 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;\r
1555 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;\r
00a6a31a 1556\r
1557\r
613bcce9 1558 for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
00a6a31a 1559 \r
613bcce9 1560 Double_t aP[3]; positiveElectrons[iPos]->GetConstrainedPxPyPz(aP); \r
00a6a31a 1561\r
613bcce9 1562 TVector3 ePosV(aP[0],aP[1],aP[2]);\r
00a6a31a 1563\r
613bcce9 1564 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){\r
00a6a31a 1565 \r
613bcce9 1566 Double_t aN[3]; negativeElectrons[iNeg]->GetConstrainedPxPyPz(aN); \r
1567 TVector3 eNegV(aN[0],aN[1],aN[2]);\r
00a6a31a 1568\r
1569 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma\r
1570 xPosBand[iPos]=kFALSE;\r
1571 xNegBand[iNeg]=kFALSE;\r
1572 }\r
1573\r
613bcce9 1574 for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
1575 AliKFParticle* gammaCandidate = &gammas[iGam];\r
1576 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());\r
1577 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)\r
00a6a31a 1578 gammaBand[iGam]=kFALSE;\r
1579 }\r
1580 }\r
1581 }\r
1582\r
1583\r
1584\r
1585\r
613bcce9 1586 for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
00a6a31a 1587 if(xPosBand[iPos]){\r
613bcce9 1588 fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);\r
00a6a31a 1589 }\r
1590 }\r
613bcce9 1591 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){\r
00a6a31a 1592 if(xNegBand[iNeg]){\r
613bcce9 1593 fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);\r
00a6a31a 1594 }\r
1595 }\r
613bcce9 1596 for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
00a6a31a 1597 if(gammaBand[iGam]){\r
613bcce9 1598 fKFReconstructedGammasCut.push_back(gammas[iGam]);\r
00a6a31a 1599 }\r
1600 }\r
1601}\r
1602\r
1603\r
1604void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)\r
1605{\r
613bcce9 1606 // see header file for documentation\r
00a6a31a 1607 pid = -1;\r
1608 weight = -1;\r
1609\r
1610 double wpart[5];\r
1611 double wpartbayes[5];\r
1612\r
1613 //get probability of the diffenrent particle types\r
1614 track->GetESDpid(wpart);\r
1615\r
1616 // Tentative particle type "concentrations"\r
1617 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};\r
1618\r
1619 //Bayes' formula\r
1620 double rcc = 0.;\r
1621 for (int i = 0; i < 5; i++)\r
1622 {\r
1623 rcc+=(c[i] * wpart[i]);\r
1624 }\r
1625\r
1626\r
1627\r
1628 for (int i=0; i<5; i++) {\r
1629 if( rcc!=0){\r
1630 wpartbayes[i] = c[i] * wpart[i] / rcc;\r
1631 }\r
1632 }\r
1633\r
1634\r
1635\r
1636 Float_t max=0.;\r
1637 int ipid=-1;\r
1638 //find most probable particle in ESD pid\r
1639 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton\r
1640 for (int i = 0; i < 5; i++)\r
1641 {\r
1642 if (wpartbayes[i] > max)\r
1643 {\r
1644 ipid = i;\r
1645 max = wpartbayes[i];\r
1646 }\r
1647 }\r
1648\r
1649 pid = ipid;\r
1650 weight = max;\r
1651}\r
1652double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)\r
1653{\r
1654 // Calculates the number of sigma to the vertex.\r
1655\r
1656 Float_t b[2];\r
1657 Float_t bRes[2];\r
1658 Float_t bCov[3];\r
1659 t->GetImpactParameters(b,bCov);\r
1660 if (bCov[0]<=0 || bCov[2]<=0) {\r
1661 AliDebug(1, "Estimated b resolution lower or equal zero!");\r
1662 bCov[0]=0; bCov[2]=0;\r
1663 }\r
1664 bRes[0] = TMath::Sqrt(bCov[0]);\r
1665 bRes[1] = TMath::Sqrt(bCov[2]);\r
1666\r
1667 // -----------------------------------\r
1668 // How to get to a n-sigma cut?\r
1669 //\r
1670 // The accumulated statistics from 0 to d is\r
1671 //\r
1672 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)\r
1673 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)\r
1674 //\r
1675 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)\r
1676 // Can this be expressed in a different way?\r
1677\r
1678 if (bRes[0] == 0 || bRes[1] ==0)\r
1679 return -1;\r
1680\r
1681 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));\r
1682\r
1683 // stupid rounding problem screws up everything:\r
1684 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(\r
1685 if (TMath::Exp(-d * d / 2) < 1e-10)\r
1686 return 1000;\r
1687\r
1688\r
1689 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);\r
1690 return d;\r
1691}\r
613bcce9 1692vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){\r
00a6a31a 1693\r
613bcce9 1694 vector <TLorentzVector> tlVtrack(0);\r
00a6a31a 1695\r
613bcce9 1696 for(UInt_t itrack=0; itrack < esdTrack.size(); itrack++){\r
1697 double P[3]; esdTrack[itrack]->GetConstrainedPxPyPz(P);\r
1698 TLorentzVector currentTrack;\r
1699 currentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);\r
1700 tlVtrack.push_back(currentTrack);\r
00a6a31a 1701 }\r
1702\r
613bcce9 1703 return tlVtrack;\r
00a6a31a 1704}\r
1705\r