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