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