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