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