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