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