1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
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 **************************************************************************/
16 ////////////////////////////////////////////////
17 //---------------------------------------------
18 // Class used to do analysis on conversion pairs
19 //---------------------------------------------
20 ////////////////////////////////////////////////
26 #include "AliAnalysisTaskGammaConversion.h"
29 #include "AliESDtrackCuts.h"
31 //#include "AliCFManager.h" // for CF
32 //#include "AliCFContainer.h" // for CF
33 #include "AliESDInputHandler.h"
34 #include "AliAnalysisManager.h"
35 #include "AliGammaConversionAODObject.h"
36 #include "AliGammaConversionBGHandler.h"
37 #include "AliESDCaloCluster.h" // for combining PHOS and GammaConv
38 #include "AliKFVertex.h"
46 class AliMCEventHandler;
47 class AliESDInputHandler;
48 class AliAnalysisManager;
55 ClassImp(AliAnalysisTaskGammaConversion)
58 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
62 fMCTruth(NULL), // for CF
63 fGCMCEvent(NULL), // for CF
65 fOutputContainer(NULL),
66 fCFManager(0x0), // for CF
68 fTriggerCINT1B(kFALSE),
70 fDoNeutralMeson(kFALSE),
71 fDoOmegaMeson(kFALSE),
74 fRecalculateV0ForGamma(kFALSE),
75 fKFReconstructedGammasTClone(NULL),
76 fKFReconstructedPi0sTClone(NULL),
77 fKFRecalculatedGammasTClone(NULL),
78 fCurrentEventPosElectronTClone(NULL),
79 fCurrentEventNegElectronTClone(NULL),
80 fKFReconstructedGammasCutTClone(NULL),
81 fPreviousEventTLVNegElectronTClone(NULL),
82 fPreviousEventTLVPosElectronTClone(NULL),
87 fElectronRecalculatedv1(),
88 fElectronRecalculatedv2(),
96 fMinOpeningAngleGhostCut(0.),
98 fCalculateBackground(kFALSE),
101 fNeutralMesonNtuple(NULL),
102 fTotalNumberOfAddedNtupleEntries(0),
103 fChargedParticles(NULL),
104 fChargedParticlesId(),
106 fMinPtForGammaJet(1.),
107 fMinIsoConeSize(0.2),
109 fMinPtGamChargedCorr(0.5),
111 fLeadingChargedIndex(-1),
116 fAODBranchName("GammaConv"),
117 fDoNeutralMesonV0MCCheck(kFALSE),
118 fKFReconstructedGammasV0Index()
120 // Default constructor
122 /* Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
123 // Common I/O in slot 0
124 DefineInput (0, TChain::Class());
125 DefineOutput(0, TTree::Class());
127 // Your private output
128 DefineOutput(1, TList::Class());
130 // Define standard ESD track cuts for Gamma-hadron correlation
135 AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
136 AliAnalysisTaskSE(name),
139 fMCTruth(NULL), // for CF
140 fGCMCEvent(NULL), // for CF
142 fOutputContainer(0x0),
143 fCFManager(0x0), // for CF
145 fTriggerCINT1B(kFALSE),
147 fDoNeutralMeson(kFALSE),
148 fDoOmegaMeson(kFALSE),
151 fRecalculateV0ForGamma(kFALSE),
152 fKFReconstructedGammasTClone(NULL),
153 fKFReconstructedPi0sTClone(NULL),
154 fKFRecalculatedGammasTClone(NULL),
155 fCurrentEventPosElectronTClone(NULL),
156 fCurrentEventNegElectronTClone(NULL),
157 fKFReconstructedGammasCutTClone(NULL),
158 fPreviousEventTLVNegElectronTClone(NULL),
159 fPreviousEventTLVPosElectronTClone(NULL),
164 fElectronRecalculatedv1(),
165 fElectronRecalculatedv2(),
173 fMinOpeningAngleGhostCut(0.),
175 fCalculateBackground(kFALSE),
176 fWriteNtuple(kFALSE),
178 fNeutralMesonNtuple(NULL),
179 fTotalNumberOfAddedNtupleEntries(0),
180 fChargedParticles(NULL),
181 fChargedParticlesId(),
183 fMinPtForGammaJet(1.),
184 fMinIsoConeSize(0.2),
186 fMinPtGamChargedCorr(0.5),
188 fLeadingChargedIndex(-1),
193 fAODBranchName("GammaConv"),
194 fDoNeutralMesonV0MCCheck(kFALSE),
195 fKFReconstructedGammasV0Index()
197 // Common I/O in slot 0
198 DefineInput (0, TChain::Class());
199 DefineOutput(0, TTree::Class());
201 // Your private output
202 DefineOutput(1, TList::Class());
203 DefineOutput(2, AliCFContainer::Class()); // for CF
206 // Define standard ESD track cuts for Gamma-hadron correlation
211 AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
213 // Remove all pointers
215 if(fOutputContainer){
216 fOutputContainer->Clear() ;
217 delete fOutputContainer ;
232 delete fEsdTrackCuts;
242 void AliAnalysisTaskGammaConversion::Init()
245 // AliLog::SetGlobalLogLevel(AliLog::kError);
247 void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
250 if (fEsdTrackCuts!=NULL){
251 delete fEsdTrackCuts;
253 fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
254 //standard cuts from:
255 //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
257 // Cuts used up to 3rd of March
259 // fEsdTrackCuts->SetMinNClustersTPC(50);
260 // fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
261 // fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
262 // fEsdTrackCuts->SetRequireITSRefit(kTRUE);
263 // fEsdTrackCuts->SetMaxNsigmaToVertex(3);
264 // fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
266 //------- To be tested-----------
268 Int_t minNClustersTPC = 70;
269 Double_t maxChi2PerClusterTPC = 4.0;
270 Double_t maxDCAtoVertexXY = 2.4; // cm
271 Double_t maxDCAtoVertexZ = 3.2; // cm
272 fEsdTrackCuts->SetRequireSigmaToVertex(kFALSE);
273 fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
274 fEsdTrackCuts->SetRequireITSRefit(kTRUE);
275 // fEsdTrackCuts->SetRequireTPCStandAlone(kTRUE);
276 fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
277 fEsdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
278 fEsdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
279 fEsdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
280 fEsdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
281 fEsdTrackCuts->SetDCAToVertex2D(kTRUE);
282 fEsdTrackCuts->SetEtaRange(-0.8, 0.8);
283 fEsdTrackCuts->SetPtRange(0.15);
285 fEsdTrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kFirst);
289 //----- From Jacek 10.03.03 ------------------/
290 // minNClustersTPC = 70;
291 // maxChi2PerClusterTPC = 4.0;
292 // maxDCAtoVertexXY = 2.4; // cm
293 // maxDCAtoVertexZ = 3.2; // cm
295 // esdTrackCuts->SetRequireSigmaToVertex(kFALSE);
296 // esdTrackCuts->SetRequireTPCRefit(kFALSE);
297 // esdTrackCuts->SetRequireTPCStandAlone(kTRUE);
298 // esdTrackCuts->SetAcceptKinkDaughters(kFALSE);
299 // esdTrackCuts->SetMinNClustersTPC(minNClustersTPC);
300 // esdTrackCuts->SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC);
301 // esdTrackCuts->SetMaxDCAToVertexXY(maxDCAtoVertexXY);
302 // esdTrackCuts->SetMaxDCAToVertexZ(maxDCAtoVertexZ);
303 // esdTrackCuts->SetDCAToVertex2D(kTRUE);
307 // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
308 // fV0Reader->SetESDtrackCuts(fEsdTrackCuts);
311 void AliAnalysisTaskGammaConversion::UserExec(Option_t */*option*/)
313 // Execute analysis for current event
315 // Load the esdpid from the esdhandler if exists (tender was applied) otherwise set the Bethe Bloch parameters
317 AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
318 AliESDInputHandler *esdHandler=0x0;
319 if ( (esdHandler=dynamic_cast<AliESDInputHandler*>(man->GetInputEventHandler())) && esdHandler->GetESDpid() ){
320 AliV0Reader::SetESDpid(esdHandler->GetESDpid());
322 //load esd pid bethe bloch parameters depending on the existance of the MC handler
323 // yes: MC parameters
324 // no: data parameters
325 if (!AliV0Reader::GetESDpid()){
327 AliV0Reader::InitESDpid();
329 AliV0Reader::InitESDpid(1);
335 if(fV0Reader == NULL){
336 // Write warning here cuts and so on are default if this ever happens
339 fV0Reader->SetInputAndMCEvent(InputEvent(), MCEvent());
341 fV0Reader->Initialize();
342 fDoMCTruth = fV0Reader->GetDoMCTruth();
346 //ConnectInputData("");
348 //Each event needs an empty branch
349 // fAODBranch->Clear();
350 fAODBranch->Delete();
352 if(fKFReconstructedGammasTClone == NULL){
353 fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
355 if(fCurrentEventPosElectronTClone == NULL){
356 fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
358 if(fCurrentEventNegElectronTClone == NULL){
359 fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
361 if(fKFReconstructedGammasCutTClone == NULL){
362 fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
364 if(fPreviousEventTLVNegElectronTClone == NULL){
365 fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
367 if(fPreviousEventTLVPosElectronTClone == NULL){
368 fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
370 if(fChargedParticles == NULL){
371 fChargedParticles = new TClonesArray("AliESDtrack",0);
374 if(fKFReconstructedPi0sTClone == NULL){
375 fKFReconstructedPi0sTClone = new TClonesArray("AliKFParticle",0);
378 if(fKFRecalculatedGammasTClone == NULL){
379 fKFRecalculatedGammasTClone = new TClonesArray("AliKFParticle",0);
384 fKFReconstructedGammasTClone->Delete();
385 fCurrentEventPosElectronTClone->Delete();
386 fCurrentEventNegElectronTClone->Delete();
387 fKFReconstructedGammasCutTClone->Delete();
388 fPreviousEventTLVNegElectronTClone->Delete();
389 fPreviousEventTLVPosElectronTClone->Delete();
390 fKFReconstructedPi0sTClone->Delete();
391 fKFRecalculatedGammasTClone->Delete();
400 fElectronRecalculatedv1.clear();
401 fElectronRecalculatedv2.clear();
404 fChargedParticles->Delete();
406 fChargedParticlesId.clear();
408 fKFReconstructedGammasV0Index.clear();
410 //Clear the data in the v0Reader
411 // fV0Reader->UpdateEventByEventData();
413 //Take Only events with proper trigger
416 if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
420 if(!fV0Reader->CheckForPrimaryVertexZ() ){
424 // Process the MC information
429 //Process the v0 information with no cuts
432 // Process the v0 information
437 FillAODWithConversionGammas() ;
439 // Process reconstructed gammas
440 if(fDoNeutralMeson == kTRUE){
441 ProcessGammasForNeutralMesonAnalysis();
445 if(fDoMCTruth == kTRUE){
448 //Process reconstructed gammas electrons for Chi_c Analysis
449 if(fDoChic == kTRUE){
450 ProcessGammaElectronsForChicAnalysis();
452 // Process reconstructed gammas for gamma Jet/hadron correlations
454 ProcessGammasForGammaJetAnalysis();
457 //calculate background if flag is set
458 if(fCalculateBackground){
459 CalculateBackground();
462 if(fDoNeutralMeson == kTRUE){
463 // ProcessConvPHOSGammasForNeutralMesonAnalysis();
464 if(fDoOmegaMeson == kTRUE){
465 ProcessGammasForOmegaMesonAnalysis();
469 //Clear the data in the v0Reader
470 fV0Reader->UpdateEventByEventData();
471 if(fRecalculateV0ForGamma==kTRUE){
472 RecalculateV0ForGamma();
474 PostData(1, fOutputContainer);
475 PostData(2, fCFManager->GetParticleContainer()); // for CF
479 // void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
480 // // see header file for documentation
481 // // printf(" ConnectInputData %s\n", GetName());
483 // AliAnalysisTaskSE::ConnectInputData(option);
485 // if(fV0Reader == NULL){
486 // // Write warning here cuts and so on are default if this ever happens
488 // fV0Reader->Initialize();
489 // fDoMCTruth = fV0Reader->GetDoMCTruth();
494 void AliAnalysisTaskGammaConversion::ProcessMCData(){
495 // see header file for documentation
496 //InputEvent(), MCEvent());
498 fStack = fV0Reader->GetMCStack();
499 fMCTruth = fV0Reader->GetMCTruth(); // for CF
500 fGCMCEvent = fV0Reader->GetMCEvent(); // for CF
502 fStack= MCEvent()->Stack();
503 fGCMCEvent=MCEvent();
506 Double_t containerInput[3];
508 if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
509 fCFManager->SetEventInfo(fGCMCEvent);
514 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
515 return; // aborts if the primary vertex does not have contributors.
518 for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
519 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
526 ///////////////////////Begin Chic Analysis/////////////////////////////
528 if(particle->GetPdgCode() == 443){//Is JPsi
529 if(particle->GetNDaughters()==2){
530 if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
531 TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
533 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
534 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
535 if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
536 fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
538 if( TMath::Abs(daug0->Eta()) < 0.9){
539 if(daug0->GetPdgCode() == -11)
540 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
542 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
545 if(TMath::Abs(daug1->Eta()) < 0.9){
546 if(daug1->GetPdgCode() == -11)
547 fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
549 fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
554 // const int CHI_C0 = 10441;
555 // const int CHI_C1 = 20443;
556 // const int CHI_C2 = 445
557 if(particle->GetPdgCode() == 22){//gamma from JPsi
558 if(particle->GetMother(0) > -1){
559 if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
560 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
561 fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
562 if(TMath::Abs(particle->Eta()) < 1.2)
563 fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
567 if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
568 if( particle->GetNDaughters() == 2){
569 TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
570 TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
572 if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
573 if( daug0->GetPdgCode() == 443){
574 TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
575 TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
576 if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
577 fHistograms->FillTable("Table_Electrons",18);
580 else if (daug1->GetPdgCode() == 443){
581 TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
582 TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
583 if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
584 fHistograms->FillTable("Table_Electrons",18);
591 /////////////////////End Chic Analysis////////////////////////////
594 if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
596 if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
598 Double_t tmpPhi=particle->Phi();
600 if(particle->Phi()> TMath::Pi()){
601 tmpPhi = particle->Phi()-(2*TMath::Pi());
605 if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
609 rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
613 if (particle->GetPdgCode() == 22){
616 if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
617 continue; // no photon as mothers!
620 if(particle->GetMother(0) >= fStack->GetNprimary()){
621 continue; // the gamma has a mother, and it is not a primary particle
624 fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
625 fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
626 fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
627 fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
628 fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
632 containerInput[0] = particle->Pt();
633 containerInput[1] = particle->Eta();
634 if(particle->GetMother(0) >=0){
635 containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
638 containerInput[2]=-1;
641 fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); // generated gamma
643 if(particle->GetMother(0) < 0){ // direct gamma
644 fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
645 fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
646 fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
647 fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
648 fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
651 // looking for conversion (electron + positron from pairbuilding (= 5) )
652 TParticle* ePos = NULL;
653 TParticle* eNeg = NULL;
655 if(particle->GetNDaughters() >= 2){
656 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
657 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
658 if(tmpDaughter->GetUniqueID() == 5){
659 if(tmpDaughter->GetPdgCode() == 11){
662 else if(tmpDaughter->GetPdgCode() == -11){
670 if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
675 Double_t ePosPhi = ePos->Phi();
676 if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
678 Double_t eNegPhi = eNeg->Phi();
679 if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
682 if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
683 continue; // no reconstruction below the Pt cut
686 if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
690 if(ePos->R()>fV0Reader->GetMaxRCut()){
691 continue; // cuts on distance from collision point
694 if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
695 continue; // outside material
699 if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
700 continue; // line cut to exclude regions where we do not reconstruct
706 fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable); // reconstructable gamma
708 fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
709 fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
710 fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
711 fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
712 fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
713 fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
715 fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
716 fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
717 fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
718 fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
720 fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
721 fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
722 fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
723 fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
727 Int_t rBin = fHistograms->GetRBin(ePos->R());
728 Int_t zBin = fHistograms->GetZBin(ePos->Vz());
729 Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
732 TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
734 TString nameMCMappingPhiR="";
735 nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
736 // fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
738 TString nameMCMappingPhi="";
739 nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
740 // fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
741 //fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
743 TString nameMCMappingR="";
744 nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
745 // fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
746 //fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
748 TString nameMCMappingPhiInR="";
749 nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
750 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
751 fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
753 TString nameMCMappingZInR="";
754 nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
755 fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
758 TString nameMCMappingPhiInZ="";
759 nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
760 // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
761 fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
765 TString nameMCMappingFMDPhiInZ="";
766 nameMCMappingFMDPhiInZ.Form("MC_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
767 fHistograms->FillHistogram(nameMCMappingFMDPhiInZ, vtxPos.Phi());
770 TString nameMCMappingRInZ="";
771 nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
772 fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
774 if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
775 TString nameMCMappingMidPtPhiInR="";
776 nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
777 fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
779 TString nameMCMappingMidPtZInR="";
780 nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
781 fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
784 TString nameMCMappingMidPtPhiInZ="";
785 nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
786 fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
790 TString nameMCMappingMidPtFMDPhiInZ="";
791 nameMCMappingMidPtFMDPhiInZ.Form("MC_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
792 fHistograms->FillHistogram(nameMCMappingMidPtFMDPhiInZ, vtxPos.Phi());
796 TString nameMCMappingMidPtRInZ="";
797 nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
798 fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
804 fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
805 fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
806 fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
807 fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
808 fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
809 fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
812 if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
813 fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
814 fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
815 fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
816 fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
817 fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
819 } // end direct gamma
820 else{ // mother exits
821 /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
822 fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
823 fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
825 fMCGammaChic.push_back(particle);
828 } // end if mother exits
829 } // end if particle is a photon
833 // process motherparticles (2 gammas as daughters)
834 // the motherparticle had already to pass the R and the eta cut, but no line cut.
835 // the line cut is just valid for the conversions!
837 if(particle->GetNDaughters() == 2){
839 TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
840 TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
842 if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
844 // Check the acceptance for both gammas
845 Bool_t gammaEtaCut = kTRUE;
846 if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
847 Bool_t gammaRCut = kTRUE;
848 if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
852 // check for conversions now -> have to pass eta, R and line cut!
853 Bool_t daughter0Electron = kFALSE;
854 Bool_t daughter0Positron = kFALSE;
855 Bool_t daughter1Electron = kFALSE;
856 Bool_t daughter1Positron = kFALSE;
858 if(daughter0->GetNDaughters() >= 2){ // first gamma
859 for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
860 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
861 if(tmpDaughter->GetUniqueID() == 5){
862 if(tmpDaughter->GetPdgCode() == 11){
863 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
864 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
865 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
866 daughter0Electron = kTRUE;
872 else if(tmpDaughter->GetPdgCode() == -11){
873 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
874 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
875 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
876 daughter0Positron = kTRUE;
886 if(daughter1->GetNDaughters() >= 2){ // second gamma
887 for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
888 TParticle *tmpDaughter = fStack->Particle(TrackIndex);
889 if(tmpDaughter->GetUniqueID() == 5){
890 if(tmpDaughter->GetPdgCode() == 11){
891 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
892 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
893 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
894 daughter1Electron = kTRUE;
900 else if(tmpDaughter->GetPdgCode() == -11){
901 if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
902 if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
903 if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
904 daughter1Positron = kTRUE;
916 if(particle->GetPdgCode()==111){ //Pi0
917 if( iTracks >= fStack->GetNprimary()){
918 fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
919 fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
920 fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
921 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
922 fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
923 fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
924 fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
925 fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
926 fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
928 if(gammaEtaCut && gammaRCut){
929 //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
930 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
931 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
932 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
933 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
934 fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
939 fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
940 fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
941 fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
942 fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
943 fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
944 fHistograms->FillHistogram("MC_Pi0_R", particle->R());
945 fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
946 fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
947 fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
948 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_Fiducial", particle->Pt());
950 if(gammaEtaCut && gammaRCut){
951 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
952 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
953 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
954 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_withinAcceptance_Fiducial", particle->Pt());
956 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
957 fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
958 fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
959 fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
960 fHistograms->FillHistogram("MC_Pi0_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
961 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
962 fHistograms->FillHistogram("MC_Pi0_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
965 if((daughter0->Energy()+daughter1->Energy())!= 0.){
966 alfa= TMath::Abs((daughter0->Energy()-daughter1->Energy())/(daughter0->Energy()+daughter1->Energy()));
968 fHistograms->FillHistogram("MC_Pi0_alpha",alfa);
969 if(TMath::Abs(particle->Eta())<0.9)fHistograms->FillHistogram("MC_Pi0_Pt_ConvGamma_withinAcceptance_Fiducial", particle->Pt());
976 if(particle->GetPdgCode()==221){ //Eta
977 fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
978 fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
979 fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
980 fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
981 fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
982 fHistograms->FillHistogram("MC_Eta_R", particle->R());
983 fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
984 fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
985 fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
987 if(gammaEtaCut && gammaRCut){
988 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
989 fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
990 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
991 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
992 fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
993 fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
994 fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
995 fHistograms->FillHistogram("MC_Eta_ConvGamma_OpeningAngle_Pt", particle->Pt(),GetMCOpeningAngle(daughter0,daughter1));
996 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter0->Pt());
997 fHistograms->FillHistogram("MC_Eta_ConvGamma_PtGamma_Pt", particle->Pt(),daughter1->Pt());
1005 // all motherparticles with 2 gammas as daughters
1006 fHistograms->FillHistogram("MC_Mother_R", particle->R());
1007 fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
1008 fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
1009 fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
1010 fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
1011 fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
1012 fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
1013 fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
1014 fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
1015 fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
1016 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
1018 if(gammaEtaCut && gammaRCut){
1019 // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
1020 fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
1021 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
1022 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
1023 if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
1024 fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
1025 fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
1026 fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
1031 } // end passed R and eta cut
1033 } // end if(particle->GetNDaughters() == 2)
1035 }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
1037 } // end ProcessMCData
1041 void AliAnalysisTaskGammaConversion::FillNtuple(){
1042 //Fills the ntuple with the different values
1044 if(fGammaNtuple == NULL){
1047 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1048 for(Int_t i=0;i<numberOfV0s;i++){
1049 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};
1050 AliESDv0 * cV0 = fV0Reader->GetV0(i);
1053 fV0Reader->GetPIDProbability(negPID,posPID);
1054 values[0]=cV0->GetOnFlyStatus();
1055 values[1]=fV0Reader->CheckForPrimaryVertex();
1058 values[4]=fV0Reader->GetX();
1059 values[5]=fV0Reader->GetY();
1060 values[6]=fV0Reader->GetZ();
1061 values[7]=fV0Reader->GetXYRadius();
1062 values[8]=fV0Reader->GetMotherCandidateNDF();
1063 values[9]=fV0Reader->GetMotherCandidateChi2();
1064 values[10]=fV0Reader->GetMotherCandidateEnergy();
1065 values[11]=fV0Reader->GetMotherCandidateEta();
1066 values[12]=fV0Reader->GetMotherCandidatePt();
1067 values[13]=fV0Reader->GetMotherCandidateMass();
1068 values[14]=fV0Reader->GetMotherCandidateWidth();
1069 // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
1070 values[16]=fV0Reader->GetOpeningAngle();
1071 values[17]=fV0Reader->GetNegativeTrackEnergy();
1072 values[18]=fV0Reader->GetNegativeTrackPt();
1073 values[19]=fV0Reader->GetNegativeTrackEta();
1074 values[20]=fV0Reader->GetNegativeTrackPhi();
1075 values[21]=fV0Reader->GetPositiveTrackEnergy();
1076 values[22]=fV0Reader->GetPositiveTrackPt();
1077 values[23]=fV0Reader->GetPositiveTrackEta();
1078 values[24]=fV0Reader->GetPositiveTrackPhi();
1079 values[25]=fV0Reader->HasSameMCMother();
1080 if(values[25] != 0){
1081 values[26]=fV0Reader->GetMotherMCParticlePDGCode();
1082 values[15]=fV0Reader->GetMotherMCParticle()->Pt();
1084 fTotalNumberOfAddedNtupleEntries++;
1085 fGammaNtuple->Fill(values);
1087 fV0Reader->ResetV0IndexNumber();
1091 void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
1092 // Process all the V0's without applying any cuts to it
1094 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
1095 for(Int_t i=0;i<numberOfV0s;i++){
1096 /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
1098 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
1102 // if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
1103 if( !fV0Reader->CheckV0FinderStatus(i)){
1108 if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
1109 !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
1113 if( fV0Reader->GetNegativeESDTrack()->GetSign()== fV0Reader->GetPositiveESDTrack()->GetSign()){
1117 if( fV0Reader->GetNegativeESDTrack()->GetKinkIndex(0) > 0 ||
1118 fV0Reader->GetPositiveESDTrack()->GetKinkIndex(0) > 0) {
1121 if(TMath::Abs(fV0Reader->GetMotherCandidateEta())> fV0Reader->GetEtaCut()){
1124 if(TMath::Abs(fV0Reader->GetPositiveTrackEta())> fV0Reader->GetEtaCut()){
1127 if(TMath::Abs(fV0Reader->GetNegativeTrackEta())> fV0Reader->GetEtaCut()){
1130 if((TMath::Abs(fV0Reader->GetZ())*fV0Reader->GetLineCutZRSlope())-fV0Reader->GetLineCutZValue() > fV0Reader->GetXYRadius() ){ // cuts out regions where we do not reconstruct
1135 if(fV0Reader->HasSameMCMother() == kFALSE){
1139 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1140 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1142 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1145 if(negativeMC->GetPdgCode() == positiveMC->GetPdgCode()){
1149 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){ // id 5 is conversion
1153 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1155 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1156 fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1157 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1158 fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1159 fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1160 fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1161 fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1162 fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1163 fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1164 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1166 fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1167 fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1169 fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1170 fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
1171 fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1172 fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1173 fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1174 fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1175 fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1176 fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1178 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1179 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1180 fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1181 fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1183 //store MCTruth properties
1184 fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1185 fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1186 fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1190 fV0Reader->ResetV0IndexNumber();
1193 void AliAnalysisTaskGammaConversion::ProcessV0s(){
1194 // see header file for documentation
1197 if(fWriteNtuple == kTRUE){
1201 Int_t nSurvivingV0s=0;
1202 while(fV0Reader->NextV0()){
1206 //-------------------------- filling v0 information -------------------------------------
1207 fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
1208 fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1209 fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1210 fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1212 fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
1213 fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
1214 fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
1215 fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
1216 fHistograms->FillHistogram("ESD_E_nTPCClusters", fV0Reader->GetNegativeTracknTPCClusters());
1217 fHistograms->FillHistogram("ESD_E_nITSClusters", fV0Reader->GetNegativeTracknITSClusters());
1219 fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
1220 fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
1221 fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
1222 fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
1223 fHistograms->FillHistogram("ESD_P_nTPCClusters", fV0Reader->GetPositiveTracknTPCClusters());
1224 fHistograms->FillHistogram("ESD_P_nITSClusters", fV0Reader->GetPositiveTracknITSClusters());
1226 fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1227 fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1228 fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1229 fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1230 fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1231 fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1232 fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1233 fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1234 fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1235 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1237 fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1238 fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1240 fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1241 fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1242 fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1243 fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1245 fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1246 fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1247 fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1248 fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1252 fV0Reader->GetPIDProbability(negPID,posPID);
1253 fHistograms->FillHistogram("ESD_ConvGamma_E_EProbP",fV0Reader->GetNegativeTrackP(),negPID);
1254 fHistograms->FillHistogram("ESD_ConvGamma_P_EProbP",fV0Reader->GetPositiveTrackP(),posPID);
1256 Double_t negPIDmupi=0;
1257 Double_t posPIDmupi=0;
1258 fV0Reader->GetPIDProbabilityMuonPion(negPIDmupi,posPIDmupi);
1259 fHistograms->FillHistogram("ESD_ConvGamma_E_mupiProbP",fV0Reader->GetNegativeTrackP(),negPIDmupi);
1260 fHistograms->FillHistogram("ESD_ConvGamma_P_mupiProbP",fV0Reader->GetPositiveTrackP(),posPIDmupi);
1262 Double_t armenterosQtAlfa[2];
1263 fV0Reader->GetArmenterosQtAlfa(fV0Reader-> GetNegativeKFParticle(),
1264 fV0Reader-> GetPositiveKFParticle(),
1265 fV0Reader->GetMotherCandidateKFCombination(),
1268 fHistograms->FillHistogram("ESD_ConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1272 Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
1273 Int_t zBin = fHistograms->GetZBin(fV0Reader->GetZ());
1274 Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
1277 TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
1279 // Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
1281 TString nameESDMappingPhiR="";
1282 nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
1283 //fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
1285 TString nameESDMappingPhi="";
1286 nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
1287 //fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
1289 TString nameESDMappingR="";
1290 nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
1291 //fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
1293 TString nameESDMappingPhiInR="";
1294 nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
1295 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1296 fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
1298 TString nameESDMappingZInR="";
1299 nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
1300 fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
1302 TString nameESDMappingPhiInZ="";
1303 nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
1304 // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
1305 fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
1307 if(fV0Reader->GetXYRadius()<rFMD){
1308 TString nameESDMappingFMDPhiInZ="";
1309 nameESDMappingFMDPhiInZ.Form("ESD_Conversion_Mapping_FMD_Phi_in_Z_%02d",zBin);
1310 fHistograms->FillHistogram(nameESDMappingFMDPhiInZ, vtxConv.Phi());
1314 TString nameESDMappingRInZ="";
1315 nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
1316 fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
1318 if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
1319 TString nameESDMappingMidPtPhiInR="";
1320 nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
1321 fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
1323 TString nameESDMappingMidPtZInR="";
1324 nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
1325 fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
1327 TString nameESDMappingMidPtPhiInZ="";
1328 nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
1329 fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
1330 if(fV0Reader->GetXYRadius()<rFMD){
1331 TString nameESDMappingMidPtFMDPhiInZ="";
1332 nameESDMappingMidPtFMDPhiInZ.Form("ESD_Conversion_Mapping_MidPt_FMD_Phi_in_Z_%02d",zBin);
1333 fHistograms->FillHistogram(nameESDMappingMidPtFMDPhiInZ, vtxConv.Phi());
1337 TString nameESDMappingMidPtRInZ="";
1338 nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
1339 fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
1346 new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
1347 fKFReconstructedGammasV0Index.push_back(fV0Reader->GetCurrentV0IndexNumber()-1);
1348 // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
1349 fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
1350 fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
1352 //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
1355 if(fV0Reader->HasSameMCMother() == kFALSE){
1358 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1359 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1361 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
1365 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
1368 if(negativeMC->GetUniqueID() == 4 && positiveMC->GetUniqueID() ==4){// fill r distribution for Dalitz decays
1369 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 111){ //pi0
1370 fHistograms->FillHistogram("ESD_TrueDalitzContamination_R", fV0Reader->GetXYRadius());
1374 if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
1378 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1381 Double_t containerInput[3];
1382 containerInput[0] = fV0Reader->GetMotherCandidatePt();
1383 containerInput[1] = fV0Reader->GetMotherCandidateEta();
1384 containerInput[2] = fV0Reader->GetMotherCandidateMass();
1385 fCFManager->GetParticleContainer()->Fill(containerInput,kStepTrueGamma); // for CF
1388 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
1389 fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
1390 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
1391 fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
1392 fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
1393 fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
1394 fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
1395 fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
1396 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
1397 fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
1398 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
1399 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
1400 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1401 fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
1403 fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
1404 fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
1407 fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
1408 fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
1409 fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
1410 fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
1412 fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
1413 fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
1414 fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
1415 fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
1417 fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
1418 fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
1419 fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
1420 fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
1421 fHistograms->FillHistogram("ESD_TrueConvGamma_alfa_qt",armenterosQtAlfa[1],armenterosQtAlfa[0]);
1425 //store MCTruth properties
1426 fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
1427 fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
1428 fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
1431 Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
1432 Double_t esdpt = fV0Reader->GetMotherCandidatePt();
1433 Double_t resdPt = 0.;
1435 resdPt = ((esdpt - mcpt)/mcpt)*100.;
1438 cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
1441 fHistograms->FillHistogram("Resolution_Gamma_dPt_Pt", mcpt, resdPt);
1442 fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
1443 fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
1444 fHistograms->FillHistogram("Resolution_Gamma_dPt_Phi", fV0Reader->GetMotherCandidatePhi(), resdPt);
1446 Double_t resdZ = 0.;
1447 if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
1448 resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100.;
1450 Double_t resdZAbs = 0.;
1451 resdZAbs = (fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz());
1453 fHistograms->FillHistogram("Resolution_dZAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdZAbs);
1454 fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
1455 fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
1456 fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
1458 // new for dPt_Pt-histograms for Electron and Positron
1459 Double_t mcEpt = fV0Reader->GetNegativeMCParticle()->Pt();
1460 Double_t resEdPt = 0.;
1462 resEdPt = ((fV0Reader->GetNegativeTrackPt()-mcEpt)/mcEpt)*100.;
1464 UInt_t statusN = fV0Reader->GetNegativeESDTrack()->GetStatus();
1465 // AliESDtrack * negTrk = fV0Reader->GetNegativeESDTrack();
1466 UInt_t kTRDoutN = (statusN & AliESDtrack::kTRDout);
1468 Int_t ITSclsE= fV0Reader->GetNegativeTracknITSClusters();
1469 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1471 case 0: // 0 ITS clusters
1472 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS0", mcEpt, resEdPt);
1474 case 1: // 1 ITS cluster
1475 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS1", mcEpt, resEdPt);
1477 case 2: // 2 ITS clusters
1478 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS2", mcEpt, resEdPt);
1480 case 3: // 3 ITS clusters
1481 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS3", mcEpt, resEdPt);
1483 case 4: // 4 ITS clusters
1484 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS4", mcEpt, resEdPt);
1486 case 5: // 5 ITS clusters
1487 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS5", mcEpt, resEdPt);
1489 case 6: // 6 ITS clusters
1490 fHistograms->FillHistogram("Resolution_E_dPt_Pt_ITS6", mcEpt, resEdPt);
1493 //Filling histograms with respect to Electron resolution
1494 fHistograms->FillHistogram("Resolution_E_dPt_Pt", mcEpt, resEdPt);
1495 fHistograms->FillHistogram("Resolution_E_dPt_Phi", fV0Reader->GetNegativeTrackPhi(), resEdPt);
1497 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1498 fHistograms->FillHistogram("Resolution_E_nTRDtracklets_MCPt", mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDntracklets());
1499 fHistograms->FillHistogram("Resolution_E_nTRDclusters_ESDPt",fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1500 fHistograms->FillHistogram("Resolution_E_nTRDclusters_MCPt",mcEpt, fV0Reader->GetNegativeESDTrack()->GetTRDncls());
1501 fHistograms->FillHistogram("Resolution_E_TRDsignal_ESDPt", fV0Reader->GetNegativeTrackPt(), fV0Reader->GetNegativeESDTrack()->GetTRDsignal());
1504 Double_t mcPpt = fV0Reader->GetPositiveMCParticle()->Pt();
1505 Double_t resPdPt = 0;
1507 resPdPt = ((fV0Reader->GetPositiveTrackPt()-mcPpt)/mcPpt)*100.;
1510 UInt_t statusP = fV0Reader->GetPositiveESDTrack()->GetStatus();
1511 // AliESDtrack * posTr= fV0Reader->GetPositiveESDTrack();
1512 UInt_t kTRDoutP = (statusP & AliESDtrack::kTRDout);
1514 Int_t ITSclsP = fV0Reader->GetPositiveTracknITSClusters();
1515 // filling Resolution_Pt_dPt with respect to the Number of ITS clusters for Positrons
1517 case 0: // 0 ITS clusters
1518 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS0", mcPpt, resPdPt);
1520 case 1: // 1 ITS cluster
1521 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS1", mcPpt, resPdPt);
1523 case 2: // 2 ITS clusters
1524 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS2", mcPpt, resPdPt);
1526 case 3: // 3 ITS clusters
1527 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS3", mcPpt, resPdPt);
1529 case 4: // 4 ITS clusters
1530 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS4", mcPpt, resPdPt);
1532 case 5: // 5 ITS clusters
1533 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS5", mcPpt, resPdPt);
1535 case 6: // 6 ITS clusters
1536 fHistograms->FillHistogram("Resolution_P_dPt_Pt_ITS6", mcPpt, resPdPt);
1539 //Filling histograms with respect to Positron resolution
1540 fHistograms->FillHistogram("Resolution_P_dPt_Pt", mcPpt, resPdPt);
1541 fHistograms->FillHistogram("Resolution_P_dPt_Phi", fV0Reader->GetPositiveTrackPhi(), resPdPt);
1543 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1544 fHistograms->FillHistogram("Resolution_P_nTRDtracklets_MCPt", mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDntracklets());
1545 fHistograms->FillHistogram("Resolution_P_nTRDclusters_ESDPt",fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1546 fHistograms->FillHistogram("Resolution_P_nTRDclusters_MCPt",mcPpt, fV0Reader->GetPositiveESDTrack()->GetTRDncls());
1547 fHistograms->FillHistogram("Resolution_P_TRDsignal_ESDPt", fV0Reader->GetPositiveTrackPt(), fV0Reader->GetPositiveESDTrack()->GetTRDsignal());
1551 Double_t resdR = 0.;
1552 if(fV0Reader->GetNegativeMCParticle()->R() != 0){
1553 resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100.;
1555 Double_t resdRAbs = 0.;
1556 resdRAbs = (fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R());
1558 fHistograms->FillHistogram("Resolution_dRAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdRAbs);
1559 fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
1560 fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
1561 fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
1562 fHistograms->FillHistogram("Resolution_R_dPt", fV0Reader->GetNegativeMCParticle()->R(), resdPt);
1564 Double_t resdPhiAbs=0.;
1566 resdPhiAbs= (fV0Reader->GetMotherCandidatePhi()-fV0Reader->GetNegativeMCParticle()->Phi());
1567 fHistograms->FillHistogram("Resolution_dPhiAbs_VS_R", fV0Reader->GetNegativeMCParticle()->R(), resdPhiAbs);
1569 }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
1571 }//while(fV0Reader->NextV0)
1572 fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
1573 fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
1574 fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
1576 fV0Reader->ResetV0IndexNumber();
1579 void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
1580 // Fill AOD with reconstructed Gamma
1582 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
1583 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
1584 //Create AOD particle object from AliKFParticle
1586 /* AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1587 //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
1588 //but this means that I have to work a little bit more in my side.
1589 //AODPWG4Particle objects are simpler and lighter, I think
1590 AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
1591 gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
1592 gamma.SetCaloLabel(-1,-1); //How to get the MC label of the 2 electrons that form the gamma?
1593 gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
1594 gamma.SetPdg(AliCaloPID::kPhotonConv); //photon id
1595 gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
1597 //Add it to the aod list
1598 Int_t i = fAODBranch->GetEntriesFast();
1599 new((*fAODBranch)[i]) AliAODPWG4Particle(gamma);
1601 // AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
1602 AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
1603 AliGammaConversionAODObject aodObject;
1604 aodObject.SetPx(gammakf->GetPx());
1605 aodObject.SetPy(gammakf->GetPy());
1606 aodObject.SetPz(gammakf->GetPz());
1607 aodObject.SetLabel1(fElectronv1[gammaIndex]);
1608 aodObject.SetLabel2(fElectronv2[gammaIndex]);
1609 Int_t i = fAODBranch->GetEntriesFast();
1610 new((*fAODBranch)[i]) AliGammaConversionAODObject(aodObject);
1615 void AliAnalysisTaskGammaConversion::ProcessGammasForOmegaMesonAnalysis(){
1616 // omega meson analysis pi0+gamma decay
1617 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
1618 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
1619 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1621 AliKFParticle * omegaCandidateGammaDaughter = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1622 if(fGammav1[firstPi0Index]==firstGammaIndex || fGammav2[firstPi0Index]==firstGammaIndex){
1626 AliKFParticle omegaCandidate(*omegaCandidatePi0Daughter,*omegaCandidateGammaDaughter);
1627 Double_t massOmegaCandidate = 0.;
1628 Double_t widthOmegaCandidate = 0.;
1630 omegaCandidate.GetMass(massOmegaCandidate,widthOmegaCandidate);
1632 fHistograms->FillHistogram("ESD_Omega_InvMass_vs_Pt",massOmegaCandidate ,omegaCandidate.GetPt());
1633 fHistograms->FillHistogram("ESD_Omega_InvMass",massOmegaCandidate);
1635 //delete omegaCandidate;
1637 }// end of omega reconstruction in pi0+gamma channel
1639 if(fDoJet == kTRUE){
1640 AliKFParticle* negPiKF=NULL;
1641 AliKFParticle* posPiKF=NULL;
1643 // look at the pi+pi+pi0 channel
1644 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
1645 AliESDtrack* posTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
1646 if (posTrack->GetSign()<0) continue;
1647 if (posPiKF) delete posPiKF; posPiKF=NULL;
1648 posPiKF = new AliKFParticle( *(posTrack) ,211);
1650 for(Int_t jCh=0;jCh<fChargedParticles->GetEntriesFast();jCh++){
1651 AliESDtrack* negTrack = (AliESDtrack*)(fChargedParticles->At(jCh));
1652 if( negTrack->GetSign()>0) continue;
1653 if (negPiKF) delete negPiKF; negPiKF=NULL;
1654 negPiKF = new AliKFParticle( *(negTrack) ,-211);
1655 AliKFParticle omegaCandidatePipPinPi0(*omegaCandidatePi0Daughter,*posPiKF,*negPiKF);
1656 Double_t massOmegaCandidatePipPinPi0 = 0.;
1657 Double_t widthOmegaCandidatePipPinPi0 = 0.;
1659 omegaCandidatePipPinPi0.GetMass(massOmegaCandidatePipPinPi0,widthOmegaCandidatePipPinPi0);
1661 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass_vs_Pt",massOmegaCandidatePipPinPi0 ,omegaCandidatePipPinPi0.GetPt());
1662 fHistograms->FillHistogram("ESD_OmegaPipPinPi0_InvMass",massOmegaCandidatePipPinPi0);
1664 // delete omegaCandidatePipPinPi0;
1667 } // checking ig gammajet because in that case the chargedparticle list is created
1673 if(fCalculateBackground){
1674 // Background calculation for the omega
1675 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
1676 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
1677 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1678 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
1679 for(Int_t firstPi0Index=0;firstPi0Index<fKFReconstructedPi0sTClone->GetEntriesFast();firstPi0Index++){
1680 AliKFParticle * omegaCandidatePi0Daughter = (AliKFParticle *)fKFReconstructedPi0sTClone->At(firstPi0Index);
1681 AliKFParticle * omegaBckCandidate = new AliKFParticle(*omegaCandidatePi0Daughter,previousGoodV0);
1682 Double_t massOmegaBckCandidate = 0.;
1683 Double_t widthOmegaBckCandidate = 0.;
1685 omegaBckCandidate->GetMass(massOmegaBckCandidate,widthOmegaBckCandidate);
1686 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass_vs_Pt",massOmegaBckCandidate ,omegaBckCandidate->GetPt());
1687 fHistograms->FillHistogram("ESD_Omega_Bck_InvMass",massOmegaBckCandidate);
1689 delete omegaBckCandidate;
1693 } // end of checking if background calculation is available
1696 void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
1697 // see header file for documentation
1699 // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
1700 // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
1702 if(fKFReconstructedGammasTClone->GetEntriesFast()>fV0Reader->GetNumberOfV0s()){
1703 cout<<"Warning, number of entries in the tclone is bigger than number of v0s"<<endl;
1706 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1707 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
1709 // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
1710 // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
1712 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1713 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
1715 if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1718 if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
1722 AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
1724 Double_t massTwoGammaCandidate = 0.;
1725 Double_t widthTwoGammaCandidate = 0.;
1726 Double_t chi2TwoGammaCandidate =10000.;
1727 twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
1728 if(twoGammaCandidate->GetNDF()>0){
1729 chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
1731 fHistograms->FillHistogram("ESD_Mother_Chi2",chi2TwoGammaCandidate);
1732 if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
1734 TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
1735 TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
1737 Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
1739 if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){
1743 rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
1746 if( (twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()) != 0){
1747 alfa=TMath::Abs((twoGammaDecayCandidateDaughter0->GetE()-twoGammaDecayCandidateDaughter1->GetE())
1748 /(twoGammaDecayCandidateDaughter0->GetE()+twoGammaDecayCandidateDaughter1->GetE()));
1751 if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
1752 delete twoGammaCandidate;
1753 continue; // minimum opening angle to avoid using ghosttracks
1756 fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
1757 fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
1758 fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
1759 fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
1760 fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
1761 fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
1762 fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
1763 fHistograms->FillHistogram("ESD_Mother_alfa", alfa);
1764 fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
1765 fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
1766 fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
1767 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1768 if(alfa<fV0Reader->GetAlphaCutMeson()){
1769 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_alpha",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1771 fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
1773 /* Kenneth, just for testing*/
1774 AliGammaConversionBGHandler * bgHandlerTest = fV0Reader->GetBGHandler();
1776 Int_t zbin= bgHandlerTest->GetZBinIndex(fV0Reader->GetVertexZ());
1777 Int_t mbintest= bgHandlerTest->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1779 fHistograms->FillHistogram(Form("%d%dESD_Mother_InvMass",zbin,mbintest),massTwoGammaCandidate);
1781 /* end Kenneth, just for testing*/
1783 if(fCalculateBackground){
1784 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
1785 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
1786 fHistograms->FillHistogram(Form("%dESD_Mother_InvMass_vs_Pt",mbin),massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1788 // if(fDoNeutralMesonV0MCCheck){
1790 //Kenneth: Checking the eta of the gamma to check the difference between 0.9 and 1.2
1791 Int_t indexKF1 = fKFReconstructedGammasV0Index.at(firstGammaIndex);
1792 if(indexKF1<fV0Reader->GetNumberOfV0s()){
1793 fV0Reader->GetV0(indexKF1);//updates to the correct v0
1794 Double_t eta1 = fV0Reader->GetMotherCandidateEta();
1795 Bool_t isRealPi0=kFALSE;
1796 Bool_t isRealEta=kFALSE;
1797 Int_t gamma1MotherLabel=-1;
1798 if(fV0Reader->HasSameMCMother() == kTRUE){
1799 //cout<<"This v0 is a real v0!!!!"<<endl;
1800 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1801 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1802 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1803 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1804 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1805 gamma1MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1810 Int_t indexKF2 = fKFReconstructedGammasV0Index.at(secondGammaIndex);
1811 if(indexKF1 == indexKF2){
1812 cout<<"index of the two KF particles are the same.... should not happen"<<endl;
1814 if(indexKF2<fV0Reader->GetNumberOfV0s()){
1815 fV0Reader->GetV0(indexKF2);
1816 Double_t eta2 = fV0Reader->GetMotherCandidateEta();
1817 Int_t gamma2MotherLabel=-1;
1818 if(fV0Reader->HasSameMCMother() == kTRUE){
1819 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
1820 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
1821 if(TMath::Abs(negativeMC->GetPdgCode())==11 && TMath::Abs(positiveMC->GetPdgCode())==11){
1822 if(negativeMC->GetUniqueID() == 5 && positiveMC->GetUniqueID() ==5){
1823 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
1824 gamma2MotherLabel=fV0Reader->GetMotherMCParticle()->GetFirstMother();
1829 if(gamma1MotherLabel>=0 && gamma1MotherLabel==gamma2MotherLabel){
1830 if(fV0Reader->CheckIfPi0IsMother(gamma1MotherLabel)){
1833 if(fV0Reader->CheckIfEtaIsMother(gamma1MotherLabel)){
1839 if(TMath::Abs(eta1)>0.9 && TMath::Abs(eta2)>0.9){
1840 // fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1841 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1842 if(isRealPi0 || isRealEta){
1843 fHistograms->FillHistogram("ESD_TruePi0_InvMass_1212",massTwoGammaCandidate);
1844 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_1212",openingAngleTwoGammaCandidate);
1845 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1846 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1847 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1850 else if(TMath::Abs(eta1)>0.9 || TMath::Abs(eta2)>0.9){
1851 // fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1852 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1853 if(isRealPi0 || isRealEta){
1854 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0912",massTwoGammaCandidate);
1855 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0912",openingAngleTwoGammaCandidate);
1856 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1857 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1858 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1862 // fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1863 // fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1864 if(isRealPi0 || isRealEta){
1865 fHistograms->FillHistogram("ESD_TruePi0_InvMass_0909",massTwoGammaCandidate);
1866 fHistograms->FillHistogram("ESD_TruePi0_OpeningAngle_0909",openingAngleTwoGammaCandidate);
1867 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1868 fHistograms->FillHistogram("ESD_TruePi0_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1869 fHistograms->FillHistogram("ESD_TruePi0_InvMass",massTwoGammaCandidate);
1875 if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
1876 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
1877 fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
1880 if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1881 fHistograms->FillHistogram("ESD_Mother_InvMass_1212",massTwoGammaCandidate);
1882 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt1212",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1884 else if(TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())>0.9 || TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())>0.9){
1885 fHistograms->FillHistogram("ESD_Mother_InvMass_0912",massTwoGammaCandidate);
1886 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0912",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1889 fHistograms->FillHistogram("ESD_Mother_InvMass_0909",massTwoGammaCandidate);
1890 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt0909",massTwoGammaCandidate,momentumVectorTwoGammaCandidate.Pt());
1892 Double_t lowMassPi0=0.1;
1893 Double_t highMassPi0=0.15;
1894 if (massTwoGammaCandidate > lowMassPi0 && massTwoGammaCandidate < highMassPi0 ){
1895 new((*fKFReconstructedPi0sTClone)[fKFReconstructedPi0sTClone->GetEntriesFast()]) AliKFParticle(*twoGammaCandidate);
1896 fGammav1.push_back(firstGammaIndex);
1897 fGammav2.push_back(secondGammaIndex);
1902 delete twoGammaCandidate;
1907 void AliAnalysisTaskGammaConversion::ProcessConvPHOSGammasForNeutralMesonAnalysis(){
1909 // see header file for documentation
1910 // Analyse Pi0 with one photon from Phos and 1 photon from conversions
1915 vtx[0] = fV0Reader->GetPrimaryVertex()->GetX();
1916 vtx[1] = fV0Reader->GetPrimaryVertex()->GetY();
1917 vtx[2] = fV0Reader->GetPrimaryVertex()->GetZ();
1920 // Loop over all CaloClusters and consider only the PHOS ones:
1921 AliESDCaloCluster *clu;
1922 TLorentzVector pPHOS;
1923 TLorentzVector gammaPHOS;
1924 TLorentzVector gammaGammaConv;
1925 TLorentzVector pi0GammaConvPHOS;
1926 TLorentzVector gammaGammaConvBck;
1927 TLorentzVector pi0GammaConvPHOSBck;
1930 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
1931 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
1932 if ( !clu->IsPHOS() || clu->E()<0.1 ) continue;
1933 clu ->GetMomentum(pPHOS ,vtx);
1934 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1935 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1936 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz(),0.);
1937 gammaPHOS.SetXYZM(pPHOS.Px(),pPHOS.Py(),pPHOS.Pz(),0.);
1938 pi0GammaConvPHOS=gammaGammaConv+gammaPHOS;
1939 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS",pi0GammaConvPHOS.M());
1940 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvPHOS",pi0GammaConvPHOS.M(),pi0GammaConvPHOS.Pt());
1942 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),twoGammaDecayCandidateDaughter0->Py(),twoGammaDecayCandidateDaughter0->Pz());
1943 TVector3 v3D1(gammaPHOS.Px(),gammaPHOS.Py(),gammaPHOS.Pz());
1944 Double_t opanConvPHOS= v3D0.Angle(v3D1);
1945 if ( opanConvPHOS < 0.35){
1946 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanLow",pi0GammaConvPHOS.M());
1948 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvPHOS_OpanHigh",pi0GammaConvPHOS.M());
1953 // Now the LorentVector pPHOS is obtained and can be paired with the converted proton
1955 //==== End of the PHOS cluster selection ============
1956 TLorentzVector pEMCAL;
1957 TLorentzVector gammaEMCAL;
1958 TLorentzVector pi0GammaConvEMCAL;
1959 TLorentzVector pi0GammaConvEMCALBck;
1961 for (Int_t i=0; i<fV0Reader->GetESDEvent()->GetNumberOfCaloClusters(); i++) {
1962 clu = fV0Reader->GetESDEvent()->GetCaloCluster(i);
1963 if ( !clu->IsEMCAL() || clu->E()<0.1 ) continue;
1964 if (clu->GetNCells() <= 1) continue;
1965 if ( clu->GetTOF()*1e9 < 550 || clu->GetTOF()*1e9 > 750) continue;
1967 clu ->GetMomentum(pEMCAL ,vtx);
1968 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
1969 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
1970 gammaGammaConv.SetXYZM(twoGammaDecayCandidateDaughter0->Px(),
1971 twoGammaDecayCandidateDaughter0->Py(),
1972 twoGammaDecayCandidateDaughter0->Pz(),0.);
1973 gammaEMCAL.SetXYZM(pEMCAL.Px(),pEMCAL.Py(),pEMCAL.Pz(),0.);
1974 pi0GammaConvEMCAL=gammaGammaConv+gammaEMCAL;
1975 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL",pi0GammaConvEMCAL.M());
1976 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL",pi0GammaConvEMCAL.M(),pi0GammaConvEMCAL.Pt());
1977 TVector3 v3D0(twoGammaDecayCandidateDaughter0->Px(),
1978 twoGammaDecayCandidateDaughter0->Py(),
1979 twoGammaDecayCandidateDaughter0->Pz());
1980 TVector3 v3D1(gammaEMCAL.Px(),gammaEMCAL.Py(),gammaEMCAL.Pz());
1983 Double_t opanConvEMCAL= v3D0.Angle(v3D1);
1984 if ( opanConvEMCAL < 0.35){
1985 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanLow",pi0GammaConvEMCAL.M());
1987 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_OpanHigh",pi0GammaConvEMCAL.M());
1991 if(fCalculateBackground){
1992 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
1993 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
1994 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
1995 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
1996 gammaGammaConvBck.SetXYZM(previousGoodV0.Px(),
1997 previousGoodV0.Py(),
1998 previousGoodV0.Pz(),0.);
1999 pi0GammaConvEMCALBck=gammaGammaConvBck+gammaEMCAL;
2000 fHistograms->FillHistogram("ESD_Mother_InvMass_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M());
2001 fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_GammaConvEMCAL_Bck",pi0GammaConvEMCALBck.M(),
2002 pi0GammaConvEMCALBck.Pt());
2006 // Now the LorentVector pEMCAL is obtained and can be paired with the converted proton
2007 } // end of checking if background photons are available
2009 //==== End of the PHOS cluster selection ============
2013 void AliAnalysisTaskGammaConversion::CalculateBackground(){
2014 // see header file for documentation
2017 TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
2019 for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
2020 AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
2021 for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
2022 AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
2023 for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
2024 AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
2026 AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
2028 Double_t massBG =0.;
2029 Double_t widthBG = 0.;
2030 Double_t chi2BG =10000.;
2031 backgroundCandidate->GetMass(massBG,widthBG);
2032 if(backgroundCandidate->GetNDF()>0){
2033 chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
2034 if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
2036 TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
2037 TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
2039 Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
2042 if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
2043 else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
2047 if( (currentEventGoodV0.GetE()+previousGoodV0.GetE()) != 0){
2048 alfa=TMath::Abs((currentEventGoodV0.GetE()-previousGoodV0.GetE())
2049 /(currentEventGoodV0.GetE()+previousGoodV0.GetE()));
2053 if(openingAngleBG < fMinOpeningAngleGhostCut ){
2054 delete backgroundCandidate;
2055 continue; // minimum opening angle to avoid using ghosttracks
2059 fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
2060 fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
2061 fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
2062 fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
2063 fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
2064 fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
2065 fHistograms->FillHistogram("ESD_Background_Mass", massBG);
2066 fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2067 fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2068 fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
2069 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
2070 fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
2072 if(alfa<fV0Reader->GetAlphaCutMeson()){
2073 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_alpha",massBG,momentumVectorbackgroundCandidate.Pt());
2075 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2076 fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
2077 fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
2082 AliGammaConversionBGHandler * bgHandler = fV0Reader->GetBGHandler();
2084 Int_t zbin= bgHandler->GetZBinIndex(fV0Reader->GetVertexZ());
2085 Int_t mbin= bgHandler->GetMultiplicityBinIndex(fV0Reader->CountESDTracks());
2087 fHistograms->FillHistogram(Form("%d%dESD_Background_GammaDaughter_OpeningAngle",zbin,mbin), openingAngleBG);
2088 fHistograms->FillHistogram(Form("%d%dESD_Background_Energy",zbin,mbin), backgroundCandidate->GetE());
2089 fHistograms->FillHistogram(Form("%d%dESD_Background_Pt",zbin,mbin), momentumVectorbackgroundCandidate.Pt());
2090 fHistograms->FillHistogram(Form("%d%dESD_Background_Eta",zbin,mbin), momentumVectorbackgroundCandidate.Eta());
2091 fHistograms->FillHistogram(Form("%d%dESD_Background_Rapidity",zbin,mbin), rapidity);
2092 fHistograms->FillHistogram(Form("%d%dESD_Background_Phi",zbin,mbin), spaceVectorbackgroundCandidate.Phi());
2093 fHistograms->FillHistogram(Form("%d%dESD_Background_Mass",zbin,mbin), massBG);
2094 fHistograms->FillHistogram(Form("%d%dESD_Background_R",zbin,mbin), spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
2095 fHistograms->FillHistogram(Form("%d%dESD_Background_ZR",zbin,mbin), backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
2096 fHistograms->FillHistogram(Form("%d%dESD_Background_XY",zbin,mbin), backgroundCandidate->GetX(), backgroundCandidate->GetY());
2097 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2098 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass",zbin,mbin),massBG);
2100 if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
2101 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_vs_Pt_Fiducial",zbin,mbin),massBG,momentumVectorbackgroundCandidate.Pt());
2102 fHistograms->FillHistogram(Form("%d%dESD_Background_InvMass_Fiducial",zbin,mbin),massBG);
2106 delete backgroundCandidate;
2113 void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
2114 //ProcessGammasForGammaJetAnalysis
2116 Double_t distIsoMin;
2118 CreateListOfChargedParticles();
2121 // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
2122 for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
2123 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
2124 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
2126 if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
2127 distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
2129 if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
2130 CalculateJetCone(gammaIndex);
2136 //____________________________________________________________________
2137 Bool_t AliAnalysisTaskGammaConversion::IsGoodImpPar(AliESDtrack *const track)
2140 // check whether particle has good DCAr(Pt) impact
2141 // parameter. Only for TPC+ITS tracks (7*sigma cut)
2142 // Origin: Andrea Dainese
2145 Float_t d0z0[2],covd0z0[3];
2146 track->GetImpactParameters(d0z0,covd0z0);
2147 Float_t sigma= 0.0050+0.0060/TMath::Power(track->Pt(),0.9);
2148 Float_t d0max = 7.*sigma;
2149 if(TMath::Abs(d0z0[0]) < d0max) return kTRUE;
2155 void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
2156 // CreateListOfChargedParticles
2158 fESDEvent = fV0Reader->GetESDEvent();
2159 Int_t numberOfESDTracks=0;
2160 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2161 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2166 if(!IsGoodImpPar(curTrack)){
2170 if(fEsdTrackCuts->AcceptTrack(curTrack) ){
2171 new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
2172 // fChargedParticles.push_back(curTrack);
2173 fChargedParticlesId.push_back(iTracks);
2174 numberOfESDTracks++;
2177 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracks",numberOfESDTracks);
2179 if (fV0Reader->GetNumberOfContributorsVtx()>=1){
2180 fHistograms->FillHistogram("ESD_NumberOfGoodESDTracksVtx",numberOfESDTracks);
2183 void AliAnalysisTaskGammaConversion::RecalculateV0ForGamma(){
2185 Double_t massE=0.00051099892;
2186 TLorentzVector curElecPos;
2187 TLorentzVector curElecNeg;
2188 TLorentzVector curGamma;
2190 TLorentzVector curGammaAt;
2191 TLorentzVector curElecPosAt;
2192 TLorentzVector curElecNegAt;
2193 AliKFVertex primVtxGamma(*(fESDEvent->GetPrimaryVertex()));
2194 AliKFVertex primVtxImprovedGamma = primVtxGamma;
2196 const AliESDVertex *vtxT3D=fESDEvent->GetPrimaryVertex();
2198 Double_t xPrimaryVertex=vtxT3D->GetXv();
2199 Double_t yPrimaryVertex=vtxT3D->GetYv();
2200 Double_t zPrimaryVertex=vtxT3D->GetZv();
2201 // Float_t primvertex[3]={xPrimaryVertex,yPrimaryVertex,zPrimaryVertex};
2203 Float_t nsigmaTPCtrackPos;
2204 Float_t nsigmaTPCtrackNeg;
2205 Float_t nsigmaTPCtrackPosToPion;
2206 Float_t nsigmaTPCtrackNegToPion;
2207 AliKFParticle* negKF=NULL;
2208 AliKFParticle* posKF=NULL;
2210 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2211 AliESDtrack* posTrack = fESDEvent->GetTrack(iTracks);
2215 if (posKF) delete posKF; posKF=NULL;
2216 if(posTrack->GetSign()<0) continue;
2217 if(!(posTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2218 if(posTrack->GetKinkIndex(0)>0 ) continue;
2219 if(posTrack->GetNcls(1)<50)continue;
2221 // posTrack->GetConstrainedPxPyPz(momPos);
2222 posTrack->GetPxPyPz(momPos);
2223 AliESDtrack *ptrk=fESDEvent->GetTrack(iTracks);
2224 curElecPos.SetXYZM(momPos[0],momPos[1],momPos[2],massE);
2225 if(TMath::Abs(curElecPos.Eta())<0.9) continue;
2226 posKF = new AliKFParticle( *(posTrack),-11);
2228 nsigmaTPCtrackPos = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kElectron);
2229 nsigmaTPCtrackPosToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(posTrack,AliPID::kPion);
2231 if ( nsigmaTPCtrackPos>5.|| nsigmaTPCtrackPos<-2.){
2235 if(pow((momPos[0]*momPos[0]+momPos[1]*momPos[1]+momPos[2]*momPos[2]),0.5)>0.5 && nsigmaTPCtrackPosToPion<1){
2241 for(Int_t jTracks = 0; jTracks < fESDEvent->GetNumberOfTracks(); jTracks++){
2242 AliESDtrack* negTrack = fESDEvent->GetTrack(jTracks);
2246 if (negKF) delete negKF; negKF=NULL;
2247 if(negTrack->GetSign()>0) continue;
2248 if(!(negTrack->GetStatus() & AliESDtrack::kTPCrefit))continue;
2249 if(negTrack->GetKinkIndex(0)>0 ) continue;
2250 if(negTrack->GetNcls(1)<50)continue;
2252 // negTrack->GetConstrainedPxPyPz(momNeg);
2253 negTrack->GetPxPyPz(momNeg);
2255 nsigmaTPCtrackNeg = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kElectron);
2256 nsigmaTPCtrackNegToPion = fV0Reader->GetESDpid()->NumberOfSigmasTPC(negTrack,AliPID::kPion);
2257 if ( nsigmaTPCtrackNeg>5. || nsigmaTPCtrackNeg<-2.){
2260 if(pow((momNeg[0]*momNeg[0]+momNeg[1]*momNeg[1]+momNeg[2]*momNeg[2]),0.5)>0.5 && nsigmaTPCtrackNegToPion<1){
2263 AliESDtrack *ntrk=fESDEvent->GetTrack(jTracks);
2264 curElecNeg.SetXYZM(momNeg[0],momNeg[1],momNeg[2],massE);
2265 if(TMath::Abs(curElecNeg.Eta())<0.9) continue;
2266 negKF = new AliKFParticle( *(negTrack) ,11);
2268 Double_t b=fESDEvent->GetMagneticField();
2269 Double_t xn, xp, dca=ntrk->GetDCA(ptrk,b,xn,xp);
2270 AliExternalTrackParam nt(*ntrk), pt(*ptrk);
2271 nt.PropagateTo(xn,b); pt.PropagateTo(xp,b);
2274 //--- Like in ITSV0Finder
2275 AliExternalTrackParam ntAt0(*ntrk), ptAt0(*ptrk);
2276 Double_t xxP,yyP,alphaP;
2279 // if (!ptAt0.GetGlobalXYZat(ptAt0->GetX(),xxP,yyP,zzP)) continue;
2280 if (!ptAt0.GetXYZAt(ptAt0.GetX(),b,rP)) continue;
2283 alphaP = TMath::ATan2(yyP,xxP);
2286 ptAt0.Propagate(alphaP,0,b);
2287 Float_t ptfacP = (1.+100.*TMath::Abs(ptAt0.GetC(b)));
2289 // Double_t distP = ptAt0.GetY();
2290 Double_t normP = ptfacP*TMath::Sqrt(ptAt0.GetSigmaY2());
2291 Double_t normdist0P = TMath::Abs(ptAt0.GetY()/normP);
2292 Double_t normdist1P = TMath::Abs((ptAt0.GetZ()-zPrimaryVertex)/(ptfacP*TMath::Sqrt(ptAt0.GetSigmaZ2())));
2293 Double_t normdistP = TMath::Sqrt(normdist0P*normdist0P+normdist1P*normdist1P);
2296 Double_t xxN,yyN,alphaN;
2298 // if (!ntAt0.GetGlobalXYZat(ntAt0->GetX(),xxN,yyN,zzN)) continue;
2299 if (!ntAt0.GetXYZAt(ntAt0.GetX(),b,rN)) continue;
2303 alphaN = TMath::ATan2(yyN,xxN);
2305 ntAt0.Propagate(alphaN,0,b);
2307 Float_t ptfacN = (1.+100.*TMath::Abs(ntAt0.GetC(b)));
2308 // Double_t distN = ntAt0.GetY();
2309 Double_t normN = ptfacN*TMath::Sqrt(ntAt0.GetSigmaY2());
2310 Double_t normdist0N = TMath::Abs(ntAt0.GetY()/normN);
2311 Double_t normdist1N = TMath::Abs((ntAt0.GetZ()-zPrimaryVertex)/(ptfacN*TMath::Sqrt(ntAt0.GetSigmaZ2())));
2312 Double_t normdistN = TMath::Sqrt(normdist0N*normdist0N+normdist1N*normdist1N);
2314 //-----------------------------
2316 Double_t momNegAt[3];
2317 nt.GetPxPyPz(momNegAt);
2318 curElecNegAt.SetXYZM(momNegAt[0],momNegAt[1],momNegAt[2],massE);
2320 Double_t momPosAt[3];
2321 pt.GetPxPyPz(momPosAt);
2322 curElecPosAt.SetXYZM(momPosAt[0],momPosAt[1],momPosAt[2],massE);
2327 // Double_t dneg= negTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2328 // Double_t dpos= posTrack->GetD(xPrimaryVertex,yPrimaryVertex,b);
2329 AliESDv0 vertex(nt,jTracks,pt,iTracks);
2332 Float_t cpa=vertex.GetV0CosineOfPointingAngle(xPrimaryVertex,yPrimaryVertex,zPrimaryVertex);
2336 // cout<< "v0Rr::"<< v0Rr<<endl;
2337 // if (pvertex.GetRr()<0.5){
2340 // cout<<"vertex.GetChi2V0()"<<vertex.GetChi2V0()<<endl;
2341 if(cpa<0.9)continue;
2342 // if (vertex.GetChi2V0() > 30) continue;
2343 // cout<<"xp+xn::"<<xp<<" "<<xn<<endl;
2344 if ((xn+xp) < 0.4) continue;
2345 if (TMath::Abs(ntrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05)
2346 if (TMath::Abs(ptrk->GetD(xPrimaryVertex,yPrimaryVertex,b))<0.05) continue;
2348 //cout<<"pass"<<endl;
2350 AliKFParticle v0GammaC;
2353 v0GammaC.SetMassConstraint(0,0.001);
2354 primVtxImprovedGamma+=v0GammaC;
2355 v0GammaC.SetProductionVertex(primVtxImprovedGamma);
2358 curGamma=curElecNeg+curElecPos;
2359 curGammaAt=curElecNegAt+curElecPosAt;
2361 // invariant mass versus pt of K0short
2363 Double_t chi2V0GammaC=100000.;
2364 if( v0GammaC.GetNDF() != 0) {
2365 chi2V0GammaC = v0GammaC.GetChi2()/v0GammaC.GetNDF();
2367 cout<< "ERROR::v0K0C.GetNDF()" << endl;
2370 if(chi2V0GammaC<200 &&chi2V0GammaC>0 ){
2371 if(fHistograms != NULL){
2372 fHistograms->FillHistogram("ESD_RecalculateV0_InvMass",v0GammaC.GetMass());
2373 fHistograms->FillHistogram("ESD_RecalculateV0_Pt",v0GammaC.GetPt());
2374 fHistograms->FillHistogram("ESD_RecalculateV0_E_dEdxP",curElecNegAt.P(),negTrack->GetTPCsignal());
2375 fHistograms->FillHistogram("ESD_RecalculateV0_P_dEdxP",curElecPosAt.P(),posTrack->GetTPCsignal());
2376 fHistograms->FillHistogram("ESD_RecalculateV0_cpa",cpa);
2377 fHistograms->FillHistogram("ESD_RecalculateV0_dca",dca);
2378 fHistograms->FillHistogram("ESD_RecalculateV0_normdistP",normdistP);
2379 fHistograms->FillHistogram("ESD_RecalculateV0_normdistN",normdistN);
2381 new((*fKFRecalculatedGammasTClone)[fKFRecalculatedGammasTClone->GetEntriesFast()]) AliKFParticle(v0GammaC);
2382 fElectronRecalculatedv1.push_back(iTracks);
2383 fElectronRecalculatedv2.push_back(jTracks);
2389 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();firstGammaIndex++){
2390 for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFRecalculatedGammasTClone->GetEntriesFast();secondGammaIndex++){
2391 AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(firstGammaIndex);
2392 AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFRecalculatedGammasTClone->At(secondGammaIndex);
2394 if(fElectronRecalculatedv1[firstGammaIndex]==fElectronRecalculatedv1[secondGammaIndex]){
2397 if( fElectronRecalculatedv2[firstGammaIndex]==fElectronRecalculatedv2[secondGammaIndex]){
2401 AliKFParticle twoGammaCandidate(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
2402 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass",twoGammaCandidate.GetMass());
2403 fHistograms->FillHistogram("ESD_RecalculateGG_InvMass_vs_Pt",twoGammaCandidate.GetMass(),twoGammaCandidate.GetPt());
2407 void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
2411 Double_t coneSize=0.3;
2414 // AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
2415 AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
2417 TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
2419 AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
2421 Double_t momLeadingCharged[3];
2422 leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
2424 TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
2426 Double_t phi1=momentumVectorLeadingCharged.Phi();
2427 Double_t eta1=momentumVectorLeadingCharged.Eta();
2428 Double_t phi3=momentumVectorCurrentGamma.Phi();
2430 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2431 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2432 Int_t chId = fChargedParticlesId[iCh];
2433 if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
2435 curTrack->GetConstrainedPxPyPz(mom);
2436 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
2437 Double_t phi2=momentumVectorChargedParticle.Phi();
2438 Double_t eta2=momentumVectorChargedParticle.Eta();
2442 if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
2443 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
2445 if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
2446 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
2448 if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
2449 cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
2453 if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
2454 ptJet+= momentumVectorChargedParticle.Pt();
2455 Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
2456 Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
2457 fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
2458 fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
2462 Double_t dphiHdrGam=phi3-phi2;
2463 if ( dphiHdrGam < (-TMath::PiOver2())){
2464 dphiHdrGam+=(TMath::TwoPi());
2467 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
2468 dphiHdrGam-=(TMath::TwoPi());
2471 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
2472 fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
2479 Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
2480 // GetMinimumDistanceToCharge
2482 Double_t fIsoMin=100.;
2483 Double_t ptLeadingCharged=-1.;
2485 fLeadingChargedIndex=-1;
2487 AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
2488 TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
2490 Double_t phi1=momentumVectorgammaHighestPt.Phi();
2491 Double_t eta1=momentumVectorgammaHighestPt.Eta();
2493 for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
2494 AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
2495 Int_t chId = fChargedParticlesId[iCh];
2496 if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
2498 curTrack->GetConstrainedPxPyPz(mom);
2499 TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
2500 Double_t phi2=momentumVectorChargedParticle.Phi();
2501 Double_t eta2=momentumVectorChargedParticle.Eta();
2502 Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
2504 if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
2510 Double_t dphiHdrGam=phi1-phi2;
2511 if ( dphiHdrGam < (-TMath::PiOver2())){
2512 dphiHdrGam+=(TMath::TwoPi());
2515 if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
2516 dphiHdrGam-=(TMath::TwoPi());
2518 if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
2519 fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
2522 if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
2523 if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
2524 ptLeadingCharged=momentumVectorChargedParticle.Pt();
2525 fLeadingChargedIndex=iCh;
2530 fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
2535 Int_t AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
2536 //GetIndexHighestPtGamma
2538 Int_t indexHighestPtGamma=-1;
2540 fGammaPtHighest = -100.;
2542 for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
2543 AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
2544 TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
2545 if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
2546 fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
2547 //gammaHighestPt = gammaHighestPtCandidate;
2548 indexHighestPtGamma=firstGammaIndex;
2552 return indexHighestPtGamma;
2557 void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
2559 // Terminate analysis
2561 AliDebug(1,"Do nothing in Terminate");
2564 void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
2567 if(fAODBranch==NULL){
2568 fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
2570 fAODBranch->Delete();
2571 fAODBranch->SetName(fAODBranchName);
2572 AddAODBranch("TClonesArray", &fAODBranch);
2574 // Create the output container
2575 if(fOutputContainer != NULL){
2576 delete fOutputContainer;
2577 fOutputContainer = NULL;
2579 if(fOutputContainer == NULL){
2580 fOutputContainer = new TList();
2581 fOutputContainer->SetOwner(kTRUE);
2584 //Adding the histograms to the output container
2585 fHistograms->GetOutputContainer(fOutputContainer);
2589 if(fGammaNtuple == NULL){
2590 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);
2592 if(fNeutralMesonNtuple == NULL){
2593 fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
2595 TList * ntupleTList = new TList();
2596 ntupleTList->SetOwner(kTRUE);
2597 ntupleTList->SetName("Ntuple");
2598 ntupleTList->Add((TNtuple*)fGammaNtuple);
2599 fOutputContainer->Add(ntupleTList);
2602 fOutputContainer->SetName(GetName());
2605 Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
2607 TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
2608 TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
2609 return v3D0.Angle(v3D1);
2612 void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
2613 // see header file for documentation
2615 vector<Int_t> indexOfGammaParticle;
2617 fStack = fV0Reader->GetMCStack();
2619 if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
2620 return; // aborts if the primary vertex does not have contributors.
2623 for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
2624 TParticle* particle = (TParticle *)fStack->Particle(iTracks);
2625 if(particle->GetPdgCode()==22){ //Gamma
2626 if(particle->GetNDaughters() >= 2){
2627 TParticle* electron=NULL;
2628 TParticle* positron=NULL;
2629 for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
2630 TParticle *tmpDaughter = fStack->Particle(daughterIndex);
2631 if(tmpDaughter->GetUniqueID() == 5){
2632 if(tmpDaughter->GetPdgCode() == 11){
2633 electron = tmpDaughter;
2635 else if(tmpDaughter->GetPdgCode() == -11){
2636 positron = tmpDaughter;
2640 if(electron!=NULL && positron!=0){
2641 if(electron->R()<160){
2642 indexOfGammaParticle.push_back(iTracks);
2649 Int_t nFoundGammas=0;
2650 Int_t nNotFoundGammas=0;
2652 Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
2653 for(Int_t i=0;i<numberOfV0s;i++){
2654 fV0Reader->GetV0(i);
2656 if(fV0Reader->HasSameMCMother() == kFALSE){
2660 TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
2661 TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
2663 if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
2666 if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
2670 if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
2671 //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
2672 for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
2673 if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
2685 void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
2686 // see header file for documantation
2688 fESDEvent = fV0Reader->GetESDEvent();
2691 TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
2692 TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
2693 TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
2694 TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
2695 TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
2696 TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
2699 vector <AliESDtrack*> vESDeNegTemp(0);
2700 vector <AliESDtrack*> vESDePosTemp(0);
2701 vector <AliESDtrack*> vESDxNegTemp(0);
2702 vector <AliESDtrack*> vESDxPosTemp(0);
2703 vector <AliESDtrack*> vESDeNegNoJPsi(0);
2704 vector <AliESDtrack*> vESDePosNoJPsi(0);
2708 fHistograms->FillTable("Table_Electrons",0);//Count number of Events
2710 for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
2711 AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
2714 //print warning here
2718 double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
2719 double r[3];curTrack->GetConstrainedXYZ(r);
2723 fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
2725 Bool_t flagKink = kTRUE;
2726 Bool_t flagTPCrefit = kTRUE;
2727 Bool_t flagTRDrefit = kTRUE;
2728 Bool_t flagITSrefit = kTRUE;
2729 Bool_t flagTRDout = kTRUE;
2730 Bool_t flagVertex = kTRUE;
2733 //Cuts ---------------------------------------------------------------
2735 if(curTrack->GetKinkIndex(0) > 0){
2736 fHistograms->FillHistogram("Table_Electrons",5);//Count kink
2740 ULong_t trkStatus = curTrack->GetStatus();
2742 ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
2745 fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
2746 flagTPCrefit = kFALSE;
2749 ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
2751 fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
2752 flagITSrefit = kFALSE;
2755 ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
2758 fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
2759 flagTRDrefit = kFALSE;
2762 ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
2765 fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
2766 flagTRDout = kFALSE;
2769 double nsigmaToVxt = GetSigmaToVertex(curTrack);
2771 if(nsigmaToVxt > 3){
2772 fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
2773 flagVertex = kFALSE;
2776 if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
2777 fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
2781 GetPID(curTrack, pid, weight);
2784 fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
2788 fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
2796 TLorentzVector curElec;
2797 curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
2801 Int_t labelMC = TMath::Abs(curTrack->GetLabel());
2802 TParticle* curParticle = fStack->Particle(labelMC);
2803 if(curTrack->GetSign() > 0){
2805 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2806 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2809 fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2810 fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2816 if(curTrack->GetSign() > 0){
2818 // vESDxPosTemp.push_back(curTrack);
2819 new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
2823 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2824 fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
2825 // fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
2826 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2827 // fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
2828 // vESDePosTemp.push_back(curTrack);
2829 new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
2836 new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
2840 fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
2841 fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
2842 fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
2843 new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
2852 Bool_t ePosJPsi = kFALSE;
2853 Bool_t eNegJPsi = kFALSE;
2854 Bool_t ePosPi0 = kFALSE;
2855 Bool_t eNegPi0 = kFALSE;
2857 UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
2859 for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
2860 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
2861 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
2862 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
2863 TParticle* partMother = fStack ->Particle(labelMother);
2864 if (partMother->GetPdgCode() == 111){
2868 if(partMother->GetPdgCode() == 443){ //Mother JPsi
2869 fHistograms->FillTable("Table_Electrons",14);
2874 // vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
2875 new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
2876 // cout<<"ESD No Positivo JPsi "<<endl;
2882 for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
2883 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
2884 if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
2885 Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
2886 TParticle* partMother = fStack ->Particle(labelMother);
2887 if (partMother->GetPdgCode() == 111){
2891 if(partMother->GetPdgCode() == 443){ //Mother JPsi
2892 fHistograms->FillTable("Table_Electrons",15);
2897 // vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
2898 new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
2899 // cout<<"ESD No Negativo JPsi "<<endl;
2905 if( eNegJPsi && ePosJPsi ){
2906 TVector3 tempeNegV,tempePosV;
2907 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
2908 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
2909 fHistograms->FillTable("Table_Electrons",16);
2910 fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
2911 fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
2912 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
2915 if( eNegPi0 && ePosPi0 ){
2916 TVector3 tempeNegV,tempePosV;
2917 tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
2918 tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
2919 fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
2920 fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
2921 fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
2925 FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
2927 CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
2929 // vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
2930 // vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
2932 TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
2933 TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
2936 FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
2941 //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
2944 FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
2945 FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
2949 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2951 FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
2952 *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
2957 ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
2958 ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
2959 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
2960 ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
2962 // Like Sign e+e- no JPsi
2963 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
2964 ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
2968 if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
2969 FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
2970 *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
2971 *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
2972 *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
2979 vtx[0]=0;vtx[1]=0;vtx[2]=0;
2980 for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
2982 // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
2986 Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
2987 // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
2989 // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
2991 if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
2993 fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
3002 vESDeNegTemp->Delete();
3003 vESDePosTemp->Delete();
3004 vESDxNegTemp->Delete();
3005 vESDxPosTemp->Delete();
3006 vESDeNegNoJPsi->Delete();
3007 vESDePosNoJPsi->Delete();
3009 delete vESDeNegTemp;
3010 delete vESDePosTemp;
3011 delete vESDxNegTemp;
3012 delete vESDxPosTemp;
3013 delete vESDeNegNoJPsi;
3014 delete vESDePosNoJPsi;
3018 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
3019 //see header file for documentation
3020 for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
3021 for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
3022 fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
3027 void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
3028 //see header file for documentation
3029 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
3030 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3031 fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
3035 void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
3036 //see header file for documentation
3037 for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
3038 TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
3039 for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
3040 TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
3041 TLorentzVector np = ep + en;
3042 fHistograms->FillHistogram(histoName.Data(),np.M());
3047 void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
3048 TClonesArray const tlVeNeg,TClonesArray const tlVePos)
3050 //see header file for documentation
3052 for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
3054 for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
3056 TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
3058 for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
3060 // AliKFParticle * gammaCandidate = &fKFGammas[iGam];
3061 AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
3064 g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
3065 TLorentzVector xyg = xy + g;
3066 fHistograms->FillHistogram(histoMass.Data(),xyg.M());
3067 fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
3073 void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
3075 // see header file for documentation
3076 for(Int_t i=0; i < e.GetEntriesFast(); i++)
3078 for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
3080 TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
3082 fHistograms->FillHistogram(hBg.Data(),ee.M());
3088 void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
3089 TClonesArray const positiveElectrons,
3090 TClonesArray const gammas){
3091 // see header file for documentation
3093 UInt_t sizeN = negativeElectrons.GetEntriesFast();
3094 UInt_t sizeP = positiveElectrons.GetEntriesFast();
3095 UInt_t sizeG = gammas.GetEntriesFast();
3099 vector <Bool_t> xNegBand(sizeN);
3100 vector <Bool_t> xPosBand(sizeP);
3101 vector <Bool_t> gammaBand(sizeG);
3104 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
3105 for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
3106 for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
3109 for(UInt_t iPos=0; iPos < sizeP; iPos++){
3112 ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
3114 TVector3 ePosV(aP[0],aP[1],aP[2]);
3116 for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
3119 ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
3120 TVector3 eNegV(aN[0],aN[1],aN[2]);
3122 if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
3123 xPosBand[iPos]=kFALSE;
3124 xNegBand[iNeg]=kFALSE;
3127 for(UInt_t iGam=0; iGam < sizeG; iGam++){
3128 AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
3129 TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
3130 if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
3131 gammaBand[iGam]=kFALSE;
3139 for(UInt_t iPos=0; iPos < sizeP; iPos++){
3141 new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
3142 // fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
3145 for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
3147 new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
3148 // fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
3151 for(UInt_t iGam=0; iGam < sizeG; iGam++){
3152 if(gammaBand[iGam]){
3153 new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
3154 //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
3160 void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
3162 // see header file for documentation
3167 double wpartbayes[5];
3169 //get probability of the diffenrent particle types
3170 track->GetESDpid(wpart);
3172 // Tentative particle type "concentrations"
3173 double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
3177 for (int i = 0; i < 5; i++)
3179 rcc+=(c[i] * wpart[i]);
3184 for (int i=0; i<5; i++) {
3185 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)
3186 wpartbayes[i] = c[i] * wpart[i] / rcc;
3194 //find most probable particle in ESD pid
3195 //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
3196 for (int i = 0; i < 5; i++)
3198 if (wpartbayes[i] > max)
3201 max = wpartbayes[i];
3208 double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
3210 // Calculates the number of sigma to the vertex.
3215 t->GetImpactParameters(b,bCov);
3216 if (bCov[0]<=0 || bCov[2]<=0) {
3217 AliDebug(1, "Estimated b resolution lower or equal zero!");
3218 bCov[0]=0; bCov[2]=0;
3220 bRes[0] = TMath::Sqrt(bCov[0]);
3221 bRes[1] = TMath::Sqrt(bCov[2]);
3223 // -----------------------------------
3224 // How to get to a n-sigma cut?
3226 // The accumulated statistics from 0 to d is
3228 // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
3229 // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
3231 // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
3232 // Can this be expressed in a different way?
3234 if (bRes[0] == 0 || bRes[1] ==0)
3237 double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
3239 // stupid rounding problem screws up everything:
3240 // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
3241 if (TMath::Exp(-d * d / 2) < 1e-10)
3245 d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
3249 //vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
3250 TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
3251 //Return TLoresntz vector of track?
3252 // vector <TLorentzVector> tlVtrack(0);
3253 TClonesArray array("TLorentzVector",0);
3255 for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
3257 //esdTrack[itrack]->GetConstrainedPxPyPz(p);
3258 ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
3259 TLorentzVector currentTrack;
3260 currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
3261 new((array)[array.GetEntriesFast()]) TLorentzVector(currentTrack);
3262 // tlVtrack.push_back(currentTrack);