Added a new class for the AOD information
[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
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;
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
659 if(GammaEtaCut && GammaRCut){
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
680 if(GammaEtaCut && GammaRCut){
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
704 if(GammaEtaCut && GammaRCut){
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
731 if(GammaEtaCut && GammaRCut){
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
1158 TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
1159 TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
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());
1175 fHistograms->FillHistogram("ESD_Background_Pt", MomentumVectorbackgroundCandidate.Pt());
1176 fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());
1177 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
1178 fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());
1179 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
1180 fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
1181 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());
1182 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
1183 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());
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];
1206 TVector3 MomentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1207
1208 if( MomentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
1209 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
1210
1211 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
1212 CalculateJetCone(gammaIndex,fLeadingChargedIndex);
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}
1235void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex, Int_t fLeadingChargedIndex){
1236 // CalculateJetCone
1237
1238 Double_t cone;
1239 Double_t coneSize=0.3;
1240 Double_t ptJet=0;
1241
1242 AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
1243 TVector3 MomentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
1244
1245 AliESDtrack* leadingCharged = fChargedParticles[fLeadingChargedIndex];
1246 Double_t momLeadingCharged[3];
1247 leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
1248
1249 TVector3 MomentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
1250
1251 Double_t phi1=MomentumVectorLeadingCharged.Phi();
1252 Double_t eta1=MomentumVectorLeadingCharged.Eta();
1253 Double_t phi3=MomentumVectorCurrentGamma.Phi();
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);
1261 TVector3 MomentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1262 Double_t phi2=MomentumVectorChargedParticle.Phi();
1263 Double_t eta2=MomentumVectorChargedParticle.Eta();
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
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();
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
1296 if (MomentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
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];
1311 TVector3 MomentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
1312
1313 Double_t phi1=MomentumVectorgammaHighestPt.Phi();
1314 Double_t eta1=MomentumVectorgammaHighestPt.Eta();
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);
1322 TVector3 MomentumVectorChargedParticle(mom[0],mom[1],mom[2]);
1323 Double_t phi2=MomentumVectorChargedParticle.Phi();
1324 Double_t eta2=MomentumVectorChargedParticle.Eta();
1325 Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
1326
1327 if(MomentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
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 }
1341 if (MomentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
1342 fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
1343 }
1344
1345 if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
1346 if (MomentumVectorChargedParticle.Pt()> ptLeadingCharged && MomentumVectorChargedParticle.Pt()>0.1*MomentumVectorgammaHighestPt.Pt()){
1347 ptLeadingCharged=MomentumVectorChargedParticle.Pt();
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;
1362 Double_t fGammaPtHighest = -100.;
1363
1364 for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1365 AliKFParticle * gammaHighestPtCandidate = &fKFReconstructedGammas[firstGammaIndex];
1366 TVector3 MomentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
1367 if (MomentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
1368 fGammaPtHighest=MomentumVectorgammaHighestPtCandidate.Pt();
1369 //gammaHighestPt = gammaHighestPtCandidate;
1370 indexHighestPtGamma=firstGammaIndex;
1371 }
1372 }
1373
1374 return indexHighestPtGamma;
1375
1376}
1377
1378
1379void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
1380{
1381 // Terminate analysis
1382 //
1383 AliDebug(1,"Do nothing in Terminate");
1384}
1385
1386void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
1387{
1388 //AOD
1389 fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
1390 fAODBranch->SetName(fAODBranchName);
1391 AddAODBranch("TClonesArray", &fAODBranch);
1392
1393 // Create the output container
1394 if(fOutputContainer != NULL){
1395 delete fOutputContainer;
1396 fOutputContainer = NULL;
1397 }
1398 if(fOutputContainer == NULL){
1399 fOutputContainer = new TList();
1400 }
1401
1402 //Adding the histograms to the output container
1403 fHistograms->GetOutputContainer(fOutputContainer);
1404
1405
1406 if(fWriteNtuple){
1407 if(fGammaNtuple == NULL){
1408 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);
1409 }
1410 if(fNeutralMesonNtuple == NULL){
1411 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
1412 }
1413 TList * ntupleTList = new TList();
1414 ntupleTList->SetName("Ntuple");
1415 ntupleTList->Add((TNtuple*)fGammaNtuple);
1416 fOutputContainer->Add(ntupleTList);
1417 }
1418
1419 fOutputContainer->SetName(GetName());
1420}
1421
1422Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
1423 //helper function
1424 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
1425 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
1426 return v3D0.Angle(v3D1);
1427}
1428
1429void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
1430 // see header file for documentation
1431
1432 vector<Int_t> indexOfGammaParticle;
1433
1434 fStack = fV0Reader->GetMCStack();
1435
1436 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1437 return; // aborts if the primary vertex does not have contributors.
1438 }
1439
1440 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
1441 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
1442 if(particle->GetPdgCode()==22){ //Gamma
1443 if(particle->GetNDaughters() >= 2){
1444 TParticle* electron=NULL;
1445 TParticle* positron=NULL;
1446 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
1447 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
1448 if(tmpDaughter->GetUniqueID() == 5){
1449 if(tmpDaughter->GetPdgCode() == 11){
1450 electron = tmpDaughter;
1451 }
1452 else if(tmpDaughter->GetPdgCode() == -11){
1453 positron = tmpDaughter;
1454 }
1455 }
1456 }
1457 if(electron!=NULL && positron!=0){
1458 if(electron->R()<160){
1459 indexOfGammaParticle.push_back(iTracks);
1460 }
1461 }
1462 }
1463 }
1464 }
1465
1466 Int_t nFoundGammas=0;
1467 Int_t nNotFoundGammas=0;
1468
1469 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1470 for(Int_t i=0;i<numberOfV0s;i++){
1471 fV0Reader->GetV0(i);
1472
1473 if(fV0Reader->HasSameMCMother() == kFALSE){
1474 continue;
1475 }
1476
1477 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1478 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1479
1480 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1481 continue;
1482 }
1483 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1484 continue;
1485 }
1486
1487 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1488 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
1489 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
1490 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
1491 nFoundGammas++;
1492 }
1493 else{
1494 nNotFoundGammas++;
1495 }
1496 }
1497 }
1498 }
1499 // cout<<"Found: "<<nFoundGammas<<" of: "<<indexOfGammaParticle.size()<<endl;
1500}
1501
1502
1503void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
1504 // see header file for documantation
1505
1506 fESDEvent = fV0Reader->GetESDEvent();
1507
1508
1509 vector <AliESDtrack*> vESDeNegTemp(0);
1510 vector <AliESDtrack*> vESDePosTemp(0);
1511 vector <AliESDtrack*> vESDxNegTemp(0);
1512 vector <AliESDtrack*> vESDxPosTemp(0);
1513 vector <AliESDtrack*> vESDeNegNoJPsi(0);
1514 vector <AliESDtrack*> vESDePosNoJPsi(0);
1515
1516
1517
1518 fHistograms->FillTable("Table_Electrons",0);//Count number of Events
1519
1520 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
1521 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
1522
1523 if(!curTrack){
1524 //print warning here
1525 continue;
1526 }
1527
1528 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
1529 double r[3];curTrack->GetConstrainedXYZ(r);
1530
1531 TVector3 rXYZ(r);
1532
1533 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
1534
1535 Bool_t flagKink = kTRUE;
1536 Bool_t flagTPCrefit = kTRUE;
1537 Bool_t flagTRDrefit = kTRUE;
1538 Bool_t flagITSrefit = kTRUE;
1539 Bool_t flagTRDout = kTRUE;
1540 Bool_t flagVertex = kTRUE;
1541
1542
1543 //Cuts ---------------------------------------------------------------
1544
1545 if(curTrack->GetKinkIndex(0) > 0){
1546 fHistograms->FillHistogram("Table_Electrons",5);//Count kink
1547 flagKink = kFALSE;
1548 }
1549
1550 ULong_t trkStatus = curTrack->GetStatus();
1551
1552 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
1553
1554 if(!tpcRefit){
1555 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
1556 flagTPCrefit = kFALSE;
1557 }
1558
1559 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
1560 if(!itsRefit){
1561 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
1562 flagITSrefit = kFALSE;
1563 }
1564
1565 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
1566
1567 if(!trdRefit){
1568 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
1569 flagTRDrefit = kFALSE;
1570 }
1571
1572 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
1573
1574 if(!trdOut) {
1575 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
1576 flagTRDout = kFALSE;
1577 }
1578
1579 double nsigmaToVxt = GetSigmaToVertex(curTrack);
1580
1581 if(nsigmaToVxt > 3){
1582 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
1583 flagVertex = kFALSE;
1584 }
1585
1586 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
1587 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
1588
1589
1590 Stat_t pid, weight;
1591 GetPID(curTrack, pid, weight);
1592
1593 if(pid!=0){
1594 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
1595 }
1596
1597 if(pid == 0){
1598 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
1599 }
1600
1601
1602
1603
1604 Int_t labelMC = TMath::Abs(curTrack->GetLabel());
1605 TParticle* curParticle = fStack->Particle(labelMC);
1606
1607
1608
1609
1610 TLorentzVector curElec;
1611 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
1612
1613
1614
1615
1616 if(curTrack->GetSign() > 0){
1617
1618 vESDxPosTemp.push_back(curTrack);
1619
1620 if( pid == 0){
1621
1622 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1623 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
1624 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1625 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1626 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
1627 vESDePosTemp.push_back(curTrack);
1628
1629
1630
1631 }
1632
1633 }
1634 else {
1635 vESDxNegTemp.push_back(curTrack);
1636
1637 if( pid == 0){
1638
1639 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
1640 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
1641 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
1642 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
1643 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
1644 vESDeNegTemp.push_back(curTrack);
1645
1646
1647
1648
1649 }
1650 }
1651
1652 }
1653
1654
1655 Bool_t ePosJPsi = kFALSE;
1656 Bool_t eNegJPsi = kFALSE;
1657 Bool_t ePosPi0 = kFALSE;
1658 Bool_t eNegPi0 = kFALSE;
1659
1660 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
1661
1662 for(UInt_t iNeg=0; iNeg < vESDeNegTemp.size(); iNeg++){
1663 if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)
1664 if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){
1665 Int_t labelMother = fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);
1666 TParticle* partMother = fStack ->Particle(labelMother);
1667 if (partMother->GetPdgCode() == 111){
1668 ieNegPi0 = iNeg;
1669 eNegPi0 = kTRUE;
1670 }
1671 if(partMother->GetPdgCode() == 443){ //Mother JPsi
1672 fHistograms->FillTable("Table_Electrons",14);
1673 ieNegJPsi = iNeg;
1674 eNegJPsi = kTRUE;
1675 }
1676 else{
1677 vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
1678 // cout<<"ESD No Positivo JPsi "<<endl;
1679 }
1680
1681 }
1682 }
1683
1684 for(UInt_t iPos=0; iPos < vESDePosTemp.size(); iPos++){
1685 if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)
1686 if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){
1687 Int_t labelMother = fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0);
1688 TParticle* partMother = fStack ->Particle(labelMother);
1689 if (partMother->GetPdgCode() == 111){
1690 iePosPi0 = iPos;
1691 ePosPi0 = kTRUE;
1692 }
1693 if(partMother->GetPdgCode() == 443){ //Mother JPsi
1694 fHistograms->FillTable("Table_Electrons",15);
1695 iePosJPsi = iPos;
1696 ePosJPsi = kTRUE;
1697 }
1698 else{
1699 vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
1700 // cout<<"ESD No Negativo JPsi "<<endl;
1701 }
1702
1703 }
1704 }
1705
1706 if( eNegJPsi && ePosJPsi ){
1707 TVector3 tempeNegV,tempePosV;
1708 tempeNegV.SetXYZ(vESDeNegTemp[ieNegJPsi]->Px(),vESDeNegTemp[ieNegJPsi]->Py(),vESDeNegTemp[ieNegJPsi]->Pz());
1709 tempePosV.SetXYZ(vESDePosTemp[iePosJPsi]->Px(),vESDePosTemp[iePosJPsi]->Py(),vESDePosTemp[iePosJPsi]->Pz());
1710 fHistograms->FillTable("Table_Electrons",16);
1711 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
1712 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegJPsi]->GetLabel())),
1713 fStack->Particle(TMath::Abs(vESDePosTemp[iePosJPsi]->GetLabel()))));
1714 }
1715
1716 if( eNegPi0 && ePosPi0 ){
1717 TVector3 tempeNegV,tempePosV;
1718 tempeNegV.SetXYZ(vESDeNegTemp[ieNegPi0]->Px(),vESDeNegTemp[ieNegPi0]->Py(),vESDeNegTemp[ieNegPi0]->Pz());
1719 tempePosV.SetXYZ(vESDePosTemp[iePosPi0]->Px(),vESDePosTemp[iePosPi0]->Py(),vESDePosTemp[iePosPi0]->Pz());
1720 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
1721 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegPi0]->GetLabel())),
1722 fStack->Particle(TMath::Abs(vESDePosTemp[iePosPi0]->GetLabel()))));
1723 }
1724
1725
1726 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
1727
1728 CleanWithAngleCuts(vESDeNegTemp,vESDePosTemp,fKFReconstructedGammas);
1729
1730 vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
1731 vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
1732
1733
1734 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
1735
1736
1737
1738
1739 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
1740
1741
1742 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
1743 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
1744
1745
1746
1747 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",
1748 fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);
1749
1750 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
1751 fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);
1752
1753 //BackGround
1754
1755 //Like Sign e+e-
1756 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
1757 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
1758 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
1759 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
1760
1761 // Like Sign e+e- no JPsi
1762 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
1763 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
1764
1765 //Mixed Event
1766
1767 if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){
1768 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
1769 fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);
1770 fPreviousEventTLVNegElectron = vCurrentTLVeNeg;
1771 fPreviousEventTLVPosElectron = vCurrentTLVePos;
1772
1773 }
1774
1775 /*
1776 //Photons P
1777 Double_t vtx[3];
1778 vtx[0]=0;vtx[1]=0;vtx[2]=0;
1779 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
1780
1781 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
1782
1783
1784
1785 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
1786 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
1787
1788 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
1789
1790 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
1791
1792 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
1793
1794
1795 }
1796
1797
1798 */
1799
1800
1801}
1802
1803void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
1804 //see header file for documentation
1805 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
1806 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
1807 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
1808 }
1809 }
1810}
1811void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){
1812 //see header file for documentation
1813 for( UInt_t n=0; n < eNeg.size(); n++){
1814
1815 TLorentzVector en = eNeg.at(n);
1816 for (UInt_t p=0; p < ePos.size(); p++){
1817 TLorentzVector ep = ePos.at(p);
1818 TLorentzVector np = ep + en;
1819 fHistograms->FillHistogram(histoName.Data(),np.M());
1820 }
1821 }
1822
1823}
1824
1825void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,
1826 vector <TLorentzVector> tlVeNeg,vector<TLorentzVector> tlVePos)
1827{
1828 //see header file for documentation
1829
1830 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++ ){
1831
1832 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
1833
1834 TLorentzVector xy = tlVePos[iPos] + tlVeNeg[iNeg];
1835
1836 for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){
1837
1838 AliKFParticle * gammaCandidate = &fKFGammas[iGam];
1839 TLorentzVector g;
1840
1841 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
1842 TLorentzVector xyg = xy + g;
1843 fHistograms->FillHistogram(histoMass.Data(),xyg.M());
1844 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
1845 }
1846 }
1847 }
1848
1849}
1850void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)
1851{
1852 // see header file for documentation
1853 for(UInt_t i=0; i < e.size(); i++)
1854 {
1855 for (UInt_t j=i+1; j < e.size(); j++)
1856 {
1857 TLorentzVector ee = e[i] + e[j];
1858
1859 fHistograms->FillHistogram(hBg.Data(),ee.M());
1860 }
1861 }
1862}
1863
1864
1865void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> negativeElectrons,
1866 vector <AliESDtrack*> positiveElectrons, vector <AliKFParticle> gammas){
1867 // see header file for documentation
1868
1869 UInt_t sizeN = negativeElectrons.size();
1870 UInt_t sizeP = positiveElectrons.size();
1871 UInt_t sizeG = gammas.size();
1872
1873
1874
1875 vector <Bool_t> xNegBand(sizeN);
1876 vector <Bool_t> xPosBand(sizeP);
1877 vector <Bool_t> gammaBand(sizeG);
1878
1879
1880 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
1881 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
1882 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
1883
1884
1885 for(UInt_t iPos=0; iPos < sizeP; iPos++){
1886
1887 Double_t aP[3]; positiveElectrons[iPos]->GetConstrainedPxPyPz(aP);
1888
1889 TVector3 ePosV(aP[0],aP[1],aP[2]);
1890
1891 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
1892
1893 Double_t aN[3]; negativeElectrons[iNeg]->GetConstrainedPxPyPz(aN);
1894 TVector3 eNegV(aN[0],aN[1],aN[2]);
1895
1896 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
1897 xPosBand[iPos]=kFALSE;
1898 xNegBand[iNeg]=kFALSE;
1899 }
1900
1901 for(UInt_t iGam=0; iGam < sizeG; iGam++){
1902 AliKFParticle* gammaCandidate = &gammas[iGam];
1903 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
1904 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
1905 gammaBand[iGam]=kFALSE;
1906 }
1907 }
1908 }
1909
1910
1911
1912
1913 for(UInt_t iPos=0; iPos < sizeP; iPos++){
1914 if(xPosBand[iPos]){
1915 fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
1916 }
1917 }
1918 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
1919 if(xNegBand[iNeg]){
1920 fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
1921 }
1922 }
1923 for(UInt_t iGam=0; iGam < sizeG; iGam++){
1924 if(gammaBand[iGam]){
1925 fKFReconstructedGammasCut.push_back(gammas[iGam]);
1926 }
1927 }
1928}
1929
1930
1931void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
1932{
1933 // see header file for documentation
1934 pid = -1;
1935 weight = -1;
1936
1937 double wpart[5];
1938 double wpartbayes[5];
1939
1940 //get probability of the diffenrent particle types
1941 track->GetESDpid(wpart);
1942
1943 // Tentative particle type "concentrations"
1944 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
1945
1946 //Bayes' formula
1947 double rcc = 0.;
1948 for (int i = 0; i < 5; i++)
1949 {
1950 rcc+=(c[i] * wpart[i]);
1951 }
1952
1953
1954
1955 for (int i=0; i<5; i++) {
1956 if( rcc!=0){
1957 wpartbayes[i] = c[i] * wpart[i] / rcc;
1958 }
1959 }
1960
1961
1962
1963 Float_t max=0.;
1964 int ipid=-1;
1965 //find most probable particle in ESD pid
1966 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
1967 for (int i = 0; i < 5; i++)
1968 {
1969 if (wpartbayes[i] > max)
1970 {
1971 ipid = i;
1972 max = wpartbayes[i];
1973 }
1974 }
1975
1976 pid = ipid;
1977 weight = max;
1978}
1979double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
1980{
1981 // Calculates the number of sigma to the vertex.
1982
1983 Float_t b[2];
1984 Float_t bRes[2];
1985 Float_t bCov[3];
1986 t->GetImpactParameters(b,bCov);
1987 if (bCov[0]<=0 || bCov[2]<=0) {
1988 AliDebug(1, "Estimated b resolution lower or equal zero!");
1989 bCov[0]=0; bCov[2]=0;
1990 }
1991 bRes[0] = TMath::Sqrt(bCov[0]);
1992 bRes[1] = TMath::Sqrt(bCov[2]);
1993
1994 // -----------------------------------
1995 // How to get to a n-sigma cut?
1996 //
1997 // The accumulated statistics from 0 to d is
1998 //
1999 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
2000 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
2001 //
2002 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
2003 // Can this be expressed in a different way?
2004
2005 if (bRes[0] == 0 || bRes[1] ==0)
2006 return -1;
2007
2008 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
2009
2010 // stupid rounding problem screws up everything:
2011 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
2012 if (TMath::Exp(-d * d / 2) < 1e-10)
2013 return 1000;
2014
2015
2016 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
2017 return d;
2018}
2019vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
2020
2021 vector <TLorentzVector> tlVtrack(0);
2022
2023 for(UInt_t itrack=0; itrack < esdTrack.size(); itrack++){
2024 double P[3]; esdTrack[itrack]->GetConstrainedPxPyPz(P);
2025 TLorentzVector currentTrack;
2026 currentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);
2027 tlVtrack.push_back(currentTrack);
2028 }
2029
2030 return tlVtrack;
2031}
2032