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