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