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