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