-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *\r
- * Version 1.1 *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-\r
-////////////////////////////////////////////////\r
-//--------------------------------------------- \r
-// Class used to do analysis on conversion pairs\r
-//---------------------------------------------\r
-////////////////////////////////////////////////\r
-\r
-// root\r
-#include <TChain.h>\r
-\r
-// analysis\r
-#include "AliAnalysisTaskGammaConversion.h"\r
-#include "AliStack.h"\r
-#include "AliLog.h"\r
-#include "TNtuple.h"\r
-\r
-class AliKFVertex;\r
-class AliAODHandler;\r
-class AliAODEvent;\r
-class ALiESDEvent;\r
-class AliMCEvent;\r
-class AliMCEventHandler;\r
-class AliESDInputHandler;\r
-class AliAnalysisManager;\r
-class Riostream;\r
-class TFile;\r
-class TInterpreter;\r
-class TSystem;\r
-class TROOT;\r
-\r
-ClassImp(AliAnalysisTaskGammaConversion)\r
-\r
-\r
-AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():\r
-AliAnalysisTaskSE(),\r
- fV0Reader(NULL),\r
- fStack(NULL),\r
- fESDEvent(NULL), \r
- fOutputContainer(NULL),\r
- fHistograms(NULL),\r
- fDoMCTruth(kFALSE),\r
- fMCAllGammas(),\r
- fMCPi0s(),\r
- fMCEtas(),\r
- fMCGammaChic(),\r
- fKFReconstructedGammas(),\r
- fIsTrueReconstructedGammas(),\r
- fElectronv1(),\r
- fElectronv2(),\r
- fCurrentEventPosElectron(),\r
- fPreviousEventPosElectron(),\r
- fCurrentEventNegElectron(),\r
- fPreviousEventNegElectron(),\r
- fKFReconstructedGammasCut(), \r
- fPreviousEventTLVNegElectron(),\r
- fPreviousEventTLVPosElectron(), \r
- fElectronMass(-1),\r
- fGammaMass(-1),\r
- fPi0Mass(-1),\r
- fEtaMass(-1),\r
- fGammaWidth(-1),\r
- fPi0Width(-1),\r
- fEtaWidth(-1),\r
- fMinOpeningAngleGhostCut(0.),\r
- fCalculateBackground(kFALSE),\r
- fWriteNtuple(kFALSE),\r
- fGammaNtuple(NULL),\r
- fNeutralMesonNtuple(NULL),\r
- fTotalNumberOfAddedNtupleEntries(0)\r
-{\r
- // Default constructor\r
- // Common I/O in slot 0\r
- DefineInput (0, TChain::Class());\r
- DefineOutput(0, TTree::Class());\r
- \r
- // Your private output\r
- DefineOutput(1, TList::Class());\r
-}\r
-\r
-AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):\r
- AliAnalysisTaskSE(name),\r
- fV0Reader(NULL),\r
- fStack(NULL),\r
- fESDEvent(NULL), \r
- fOutputContainer(0x0),\r
- fHistograms(NULL),\r
- fDoMCTruth(kFALSE),\r
- fMCAllGammas(),\r
- fMCPi0s(),\r
- fMCEtas(),\r
- fMCGammaChic(),\r
- fKFReconstructedGammas(),\r
- fIsTrueReconstructedGammas(),\r
- fElectronv1(),\r
- fElectronv2(),\r
- fCurrentEventPosElectron(),\r
- fPreviousEventPosElectron(),\r
- fCurrentEventNegElectron(),\r
- fPreviousEventNegElectron(),\r
- fKFReconstructedGammasCut(), \r
- fPreviousEventTLVNegElectron(),\r
- fPreviousEventTLVPosElectron(),\r
- fElectronMass(-1),\r
- fGammaMass(-1),\r
- fPi0Mass(-1),\r
- fEtaMass(-1),\r
- fGammaWidth(-1),\r
- fPi0Width(-1),\r
- fEtaWidth(-1),\r
- fMinOpeningAngleGhostCut(0.),\r
- fCalculateBackground(kFALSE),\r
- fWriteNtuple(kFALSE),\r
- fGammaNtuple(NULL),\r
- fNeutralMesonNtuple(NULL),\r
- fTotalNumberOfAddedNtupleEntries(0)\r
-{\r
- // Common I/O in slot 0\r
- DefineInput (0, TChain::Class());\r
- DefineOutput(0, TTree::Class());\r
- \r
- // Your private output\r
- DefineOutput(1, TList::Class());\r
-}\r
-\r
-AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion() \r
-{\r
- // Remove all pointers\r
- \r
- if(fOutputContainer){\r
- fOutputContainer->Clear() ; \r
- delete fOutputContainer ;\r
- }\r
- if(fHistograms){\r
- delete fHistograms;\r
- }\r
- if(fV0Reader){\r
- delete fV0Reader;\r
- }\r
-}\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::Init()\r
-{\r
- // Initialization\r
- AliLog::SetGlobalLogLevel(AliLog::kError);\r
-}\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)\r
-{\r
- // Execute analysis for current event\r
- \r
- ConnectInputData("");\r
- \r
- //clear vectors\r
- fMCAllGammas.clear();\r
- fMCPi0s.clear();\r
- fMCEtas.clear();\r
- fMCGammaChic.clear();\r
- \r
- fKFReconstructedGammas.clear();\r
- fIsTrueReconstructedGammas.clear();\r
- fElectronv1.clear();\r
- fElectronv2.clear();\r
- fCurrentEventPosElectron.clear();\r
- fCurrentEventNegElectron.clear(); \r
- fKFReconstructedGammasCut.clear(); \r
- \r
- //Clear the data in the v0Reader\r
- fV0Reader->UpdateEventByEventData();\r
-\r
- \r
- // Process the MC information\r
- if(fDoMCTruth){\r
- ProcessMCData();\r
- }\r
- \r
- //Process the v0 information with no cuts\r
- ProcessV0sNoCut();\r
- \r
- // Process the v0 information\r
- ProcessV0s();\r
- \r
- //calculate background if flag is set\r
- if(fCalculateBackground){\r
- CalculateBackground();\r
- }\r
- \r
- // Process reconstructed gammas\r
- ProcessGammasForNeutralMesonAnalysis();\r
-\r
- CheckV0Efficiency();\r
- \r
- //Process reconstructed gammas electrons for Chi_c Analysis\r
- ProcessGammaElectronsForChicAnalysis();\r
-\r
- PostData(1, fOutputContainer);\r
- \r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t */*option*/){\r
- // see header file for documentation\r
- \r
- if(fV0Reader == NULL){\r
- // Write warning here cuts and so on are default if this ever happens\r
- }\r
- fV0Reader->Initialize();\r
-}\r
-\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::ProcessMCData(){\r
- // see header file for documentation\r
- \r
- fStack = fV0Reader->GetMCStack();\r
-\r
- if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
- return; // aborts if the primary vertex does not have contributors.\r
- }\r
-\r
- for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {\r
- TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
- \r
- if (!particle) {\r
- //print warning here\r
- continue;\r
- }\r
-\r
- ///////////////////////Begin Chic Analysis/////////////////////////////\r
-\r
-\r
- if(particle->GetPdgCode() == 443){//Is JPsi\r
-\r
- if(particle->GetNDaughters()==2){\r
- if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&\r
- TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){\r
- TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
- TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
- if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)\r
- fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance\r
-\r
- if( TMath::Abs(daug0->Eta()) < 0.9){\r
- if(daug0->GetPdgCode() == -11)\r
- fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
- else\r
- fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
-\r
- }\r
- if(TMath::Abs(daug1->Eta()) < 0.9){\r
- if(daug1->GetPdgCode() == -11)\r
- fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance\r
- else\r
- fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance\r
- }\r
- }\r
- }\r
- }\r
- // const int CHI_C0 = 10441;\r
- // const int CHI_C1 = 20443;\r
- // const int CHI_C2 = 445\r
- if(particle->GetPdgCode() == 22){//gamma from JPsi\r
- if(particle->GetMother(0) > -1){\r
- if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||\r
- fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||\r
- fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){\r
- if(TMath::Abs(particle->Eta()) < 1.2)\r
- fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance\r
- }\r
- }\r
- }\r
- if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){\r
- if( particle->GetNDaughters() == 2){\r
- TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());\r
- TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());\r
-\r
- if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){\r
- if( daug0->GetPdgCode() == 443){\r
- TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());\r
- TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());\r
- if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
- fHistograms->FillTable("Table_Electrons",18);\r
-\r
- }//if\r
- else if (daug1->GetPdgCode() == 443){\r
- TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());\r
- TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());\r
- if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )\r
- fHistograms->FillTable("Table_Electrons",18);\r
- }//else if\r
- }//gamma o Jpsi\r
- }//GetNDaughters\r
- }\r
-\r
-\r
- /////////////////////End Chic Analysis////////////////////////////\r
-\r
-\r
- if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;\r
- \r
- if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point\r
- \r
- Double_t tmpPhi=particle->Phi();\r
- \r
- if(particle->Phi()> TMath::Pi()){\r
- tmpPhi = particle->Phi()-(2*TMath::Pi());\r
- }\r
- \r
- Double_t rapidity;\r
- if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){\r
- rapidity=0;\r
- }\r
- else{\r
- rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));\r
- } \r
- \r
- //process the gammas\r
- if (particle->GetPdgCode() == 22){\r
- \r
- if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){\r
- continue; // no photon as mothers!\r
- }\r
-\r
- if(particle->GetMother(0) >= fStack->GetNprimary()){\r
- continue; // the gamma has a mother, and it is not a primary particle\r
- }\r
-\r
- fMCAllGammas.push_back(particle);\r
- \r
- fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);\r
- \r
- \r
- if(particle->GetMother(0) < 0){ // direct gamma\r
- fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());\r
- fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity); \r
- }\r
- \r
- \r
- // looking for conversion (electron + positron from pairbuilding (= 5) )\r
- TParticle* ePos = NULL;\r
- TParticle* eNeg = NULL;\r
- \r
- if(particle->GetNDaughters() >= 2){\r
- for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
- TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
- if(tmpDaughter->GetUniqueID() == 5){\r
- if(tmpDaughter->GetPdgCode() == 11){\r
- eNeg = tmpDaughter;\r
- }\r
- else if(tmpDaughter->GetPdgCode() == -11){\r
- ePos = tmpDaughter;\r
- }\r
- }\r
- }\r
- }\r
- \r
- \r
- if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production\r
- continue;\r
- }\r
- \r
- \r
- Double_t ePosPhi = ePos->Phi();\r
- if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());\r
- \r
- Double_t eNegPhi = eNeg->Phi();\r
- if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());\r
- \r
- \r
- if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){\r
- continue; // no reconstruction below the Pt cut\r
- }\r
- \r
- if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){\r
- continue;\r
- } \r
- \r
- if(ePos->R()>fV0Reader->GetMaxRCut()){\r
- continue; // cuts on distance from collision point\r
- }\r
- \r
- \r
- if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){\r
- continue; // line cut to exclude regions where we do not reconstruct\r
- } \r
- \r
- fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);\r
- fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());\r
- \r
- fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());\r
- fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());\r
- fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());\r
- fHistograms->FillHistogram("MC_E_Phi", eNegPhi);\r
- \r
- fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());\r
- fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());\r
- fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());\r
- fHistograms->FillHistogram("MC_P_Phi", ePosPhi);\r
- \r
- \r
- \r
- //cout << "filled histos for converted gamma, ePos, eNeg" << endl;\r
- \r
- // begin Mapping \r
- Int_t rBin = fHistograms->GetRBin(ePos->R());\r
- Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());\r
- \r
- TString nameMCMappingPhiR="";\r
- nameMCMappingPhiR.Form("MC_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
- fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());\r
- \r
- TString nameMCMappingPhi="";\r
- nameMCMappingPhi.Form("MC_Conversion_Mapping-Phi%02d",phiBin);\r
- fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());\r
- \r
- TString nameMCMappingR="";\r
- nameMCMappingR.Form("MC_Conversion_Mapping-R%02d",rBin);\r
- fHistograms->FillHistogram(nameMCMappingR, particle->Eta());\r
- \r
- TString nameMCMappingPhiInR="";\r
- nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_R-%02d",rBin);\r
- fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);\r
- //end mapping\r
- \r
- fHistograms->FillHistogram("MC_Conversion_R",ePos->R());\r
- fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());\r
- fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());\r
- fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));\r
- \r
- //cout << "mapping is done" << endl;\r
- \r
- \r
- if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted\r
- fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());\r
- fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);\r
- \r
- } // end direct gamma\r
- else{ // mother exits \r
- if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0 \r
- fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S\r
- fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2\r
- ){ \r
- fMCGammaChic.push_back(particle);\r
- }\r
- } // end if mother exits\r
- } // end if particle is a photon\r
- \r
- if(particle->GetNDaughters() == 2){\r
- \r
- TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());\r
- TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());\r
- \r
- if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters\r
- \r
- \r
- \r
- // check for conversions now -> have to pass eta and line cut!\r
- Bool_t daughter0Electron = kFALSE;\r
- Bool_t daughter0Positron = kFALSE;\r
- Bool_t daughter1Electron = kFALSE;\r
- Bool_t daughter1Positron = kFALSE;\r
- \r
- \r
- \r
- if(daughter0->GetNDaughters() >= 2){\r
- for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){\r
- TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
- if(tmpDaughter->GetUniqueID() == 5){\r
- if(tmpDaughter->GetPdgCode() == 11){\r
- if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
- if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
- daughter0Electron = kTRUE;\r
- }\r
- \r
- }\r
- }\r
- else if(tmpDaughter->GetPdgCode() == -11){\r
- if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
- if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
- daughter0Positron = kTRUE;\r
- }\r
- \r
- }\r
- \r
- }\r
- }\r
- }\r
- }\r
- \r
- \r
- \r
- if(daughter1->GetNDaughters() >= 2){\r
- for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){\r
- TParticle *tmpDaughter = fStack->Particle(TrackIndex);\r
- if(tmpDaughter->GetUniqueID() == 5){\r
- if(tmpDaughter->GetPdgCode() == 11){\r
- if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
- if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
- daughter1Electron = kTRUE;\r
- }\r
- \r
- }\r
- }\r
- else if(tmpDaughter->GetPdgCode() == -11){\r
- if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){\r
- if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){\r
- daughter1Positron = kTRUE;\r
- }\r
- \r
- }\r
- \r
- }\r
- }\r
- }\r
- }\r
- \r
- \r
- \r
- \r
- if(particle->GetPdgCode()==111){ //Pi0\r
- if( iTracks >= fStack->GetNprimary()){\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling\r
- \r
- if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
- if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
- }\r
- }\r
- }\r
- else{\r
- fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta()); \r
- fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);\r
- fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);\r
- fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_Pi0_R", particle->R());\r
- fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());\r
- fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
- fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
- \r
- if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
- fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
- if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
- fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
- }\r
- }\r
- }\r
- }\r
- \r
- if(particle->GetPdgCode()==221){ //Eta\r
- fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);\r
- fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);\r
- fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_Eta_R", particle->R());\r
- fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());\r
- fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
- fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling\r
- \r
- if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
- fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
- if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
- fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
- }\r
- \r
- }\r
- \r
- }\r
- \r
- // all motherparticles with 2 gammas as daughters\r
- fHistograms->FillHistogram("MC_Mother_R", particle->R());\r
- fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());\r
- fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());\r
- fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());\r
- fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));\r
- fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());\r
- fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());\r
- fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());\r
- fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);\r
- fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);\r
- fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt()); \r
- if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){\r
- fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);\r
- fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt()); \r
- if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){\r
- fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());\r
- fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);\r
- fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt()); \r
-\r
- }\r
- \r
- \r
- }\r
- \r
- //cout << "mother histos are filled" << endl;\r
- \r
- } // end if(particle->GetNDaughters() == 2)\r
- \r
- }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)\r
- \r
- //cout << "right before the end of processMCdata" << endl;\r
- \r
-} // end ProcessMCData\r
-\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::FillNtuple(){\r
- //Fills the ntuple with the different values\r
-\r
- if(fGammaNtuple == NULL){\r
- return;\r
- }\r
- Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
- for(Int_t i=0;i<numberOfV0s;i++){\r
- 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};\r
- AliESDv0 * cV0 = fV0Reader->GetV0(i);\r
- Double_t negPID=0;\r
- Double_t posPID=0;\r
- fV0Reader->GetPIDProbability(negPID,posPID);\r
- values[0]=cV0->GetOnFlyStatus();\r
- values[1]=fV0Reader->CheckForPrimaryVertex();\r
- values[2]=negPID;\r
- values[3]=posPID;\r
- values[4]=fV0Reader->GetX();\r
- values[5]=fV0Reader->GetY();\r
- values[6]=fV0Reader->GetZ();\r
- values[7]=fV0Reader->GetXYRadius();\r
- values[8]=fV0Reader->GetMotherCandidateNDF();\r
- values[9]=fV0Reader->GetMotherCandidateChi2();\r
- values[10]=fV0Reader->GetMotherCandidateEnergy();\r
- values[11]=fV0Reader->GetMotherCandidateEta();\r
- values[12]=fV0Reader->GetMotherCandidatePt();\r
- values[13]=fV0Reader->GetMotherCandidateMass();\r
- values[14]=fV0Reader->GetMotherCandidateWidth();\r
- // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array\r
- values[16]=fV0Reader->GetOpeningAngle();\r
- values[17]=fV0Reader->GetNegativeTrackEnergy();\r
- values[18]=fV0Reader->GetNegativeTrackPt();\r
- values[19]=fV0Reader->GetNegativeTrackEta();\r
- values[20]=fV0Reader->GetNegativeTrackPhi();\r
- values[21]=fV0Reader->GetPositiveTrackEnergy();\r
- values[22]=fV0Reader->GetPositiveTrackPt();\r
- values[23]=fV0Reader->GetPositiveTrackEta();\r
- values[24]=fV0Reader->GetPositiveTrackPhi();\r
- values[25]=fV0Reader->HasSameMCMother();\r
- if(values[25] != 0){\r
- values[26]=fV0Reader->GetMotherMCParticlePDGCode();\r
- values[15]=fV0Reader->GetMotherMCParticle()->Pt();\r
- }\r
- fTotalNumberOfAddedNtupleEntries++;\r
- fGammaNtuple->Fill(values);\r
- }\r
- fV0Reader->ResetV0IndexNumber();\r
- \r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){\r
- // Process all the V0's without applying any cuts to it\r
-\r
- Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
- for(Int_t i=0;i<numberOfV0s;i++){\r
- /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);\r
-\r
- if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
- return;\r
- }\r
- \r
- if(fDoMCTruth){\r
- \r
- if(fV0Reader->HasSameMCMother() == kFALSE){\r
- continue;\r
- }\r
- \r
- TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
- TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
-\r
- if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
- continue;\r
- }\r
- if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
- continue;\r
- }\r
- \r
- if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
- \r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta()); \r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
- fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
- \r
- fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
- fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
- \r
- /*\r
- ESD_NoCutConvGamma_Pt\r
- ESD_NoCutConvGamma_Energy\r
- ESD_NoCutConvGamma_Eta\r
- ESD_NoCutConvGamma_Phi\r
- ESD_NoCutConvGamma_Mass\r
- ESD_NoCutConvGamma_Width\r
- ESD_NoCutConvGamma_Chi2\r
- ESD_NoCutConvGamma_NDF\r
- ESD_NoCutConvGamma_PtvsEta\r
- ESD_NoCutConversion_XY\r
- ESD_NoCutConversion_R\r
- ESD_NoCutConversion_ZR\r
- ESD_NoCutConversion_OpeningAngle\r
- */\r
- }\r
- }\r
- }\r
- fV0Reader->ResetV0IndexNumber();\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::ProcessV0s(){\r
- // see header file for documentation\r
- \r
- if(fWriteNtuple == kTRUE){\r
- FillNtuple();\r
- }\r
- \r
- Int_t nSurvivingV0s=0;\r
- while(fV0Reader->NextV0()){\r
- nSurvivingV0s++;\r
- \r
- \r
- //-------------------------- filling v0 information -------------------------------------\r
- fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
- fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle()); \r
- \r
- fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());\r
- fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());\r
- fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());\r
- fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());\r
- \r
- fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());\r
- fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());\r
- fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());\r
- fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());\r
- \r
- fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
- fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
- fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
- \r
- \r
- // begin mapping\r
- Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());\r
- Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());\r
- Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();\r
- \r
- TString nameESDMappingPhiR="";\r
- nameESDMappingPhiR.Form("ESD_Conversion_Mapping-Phi%02d-R%02d",phiBin,rBin);\r
- fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);\r
- \r
- TString nameESDMappingPhi="";\r
- nameESDMappingPhi.Form("ESD_Conversion_Mapping-Phi%02d",phiBin);\r
- fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);\r
- \r
- TString nameESDMappingR="";\r
- nameESDMappingR.Form("ESD_Conversion_Mapping-R%02d",rBin);\r
- fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta); \r
- \r
- TString nameESDMappingPhiInR="";\r
- nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_R-%02d",rBin);\r
- fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());\r
- // end mapping\r
- \r
- fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());\r
- fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());\r
- fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());\r
-\r
- \r
- //----------------------------------- checking for "real" conversions (MC match) --------------------------------------\r
- if(fDoMCTruth){\r
- \r
- if(fV0Reader->HasSameMCMother() == kFALSE){\r
- fIsTrueReconstructedGammas.push_back(kFALSE);\r
- continue;\r
- }\r
- \r
- \r
- TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
- TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
-\r
- if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
- fIsTrueReconstructedGammas.push_back(kFALSE);\r
- continue;\r
- }\r
- if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
- fIsTrueReconstructedGammas.push_back(kFALSE);\r
- continue;\r
- }\r
- \r
- if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
- fIsTrueReconstructedGammas.push_back(kTRUE);\r
- \r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta()); \r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
- fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());\r
- \r
- fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
- fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
-\r
- \r
- /*\r
- fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());\r
- fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());\r
- */\r
-\r
-\r
-\r
- //resolution\r
- Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();\r
- Double_t esdpt = fV0Reader->GetMotherCandidatePt();\r
- Double_t resdPt = 0;\r
- if(mcpt != 0){\r
- resdPt = ((esdpt - mcpt)/mcpt)*100;\r
- }\r
- \r
- fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);\r
- fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);\r
- fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);\r
- \r
- Double_t resdZ = 0;\r
- if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){\r
- resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;\r
- }\r
- \r
- fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);\r
- fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());\r
- fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());\r
- \r
- Double_t resdR = 0;\r
- if(fV0Reader->GetNegativeMCParticle()->R() != 0){\r
- resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;\r
- }\r
-\r
- fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);\r
- fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());\r
- fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());\r
- fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);\r
- }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)\r
- else{\r
- fIsTrueReconstructedGammas.push_back(kFALSE);\r
- }\r
- }//if(fDoMCTruth)\r
- }//while(fV0Reader->NextV0)\r
- fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);\r
- fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());\r
- \r
- //cout << "nearly at the end of doMCTruth" << endl;\r
- \r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){\r
- // see header file for documentation\r
- \r
- for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){\r
- for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){\r
- \r
- AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];\r
- AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];\r
- \r
- if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
- continue;\r
- }\r
- if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){\r
- continue;\r
- }\r
-\r
- /*\r
- if(fIsTrueReconstructedGammas[firstGammaIndex] == kFALSE || fIsTrueReconstructedGammas[secondGammaIndex] == kFALSE){\r
- continue;\r
- }\r
- */\r
-\r
- AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);\r
- \r
- Double_t massTwoGammaCandidate = 0.;\r
- Double_t widthTwoGammaCandidate = 0.;\r
- Double_t chi2TwoGammaCandidate =10000.; \r
- twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);\r
- if(twoGammaCandidate->GetNDF()>0){\r
- chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();\r
- \r
- if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){ \r
- \r
- TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());\r
- TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());\r
- \r
- Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1); \r
- Double_t rapidity;\r
- if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){\r
- rapidity=0;\r
- }\r
- else{\r
- rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));\r
- }\r
- \r
- if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut) continue; // minimum opening angle to avoid using ghosttracks\r
-\r
- fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);\r
- fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());\r
- fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());\r
- fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity); \r
- fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());\r
- fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);\r
- fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!\r
- fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());\r
- fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);\r
- }\r
- }\r
- delete twoGammaCandidate;\r
- \r
- //cout << "nearly at the end of processgamma for neutral meson ..." << endl;\r
- \r
- \r
- }\r
- }\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::CalculateBackground(){\r
- // see header file for documentation\r
- \r
- vector<AliKFParticle> vectorCurrentEventGoodV0s = fV0Reader->GetCurrentEventGoodV0s();\r
- vector<AliKFParticle> vectorPreviousEventGoodV0s = fV0Reader->GetPreviousEventGoodV0s();\r
- for(UInt_t iCurrent=0;iCurrent<vectorCurrentEventGoodV0s.size();iCurrent++){\r
- AliKFParticle * currentEventGoodV0 = &vectorCurrentEventGoodV0s.at(iCurrent);\r
- for(UInt_t iPrevious=0;iPrevious<vectorPreviousEventGoodV0s.size();iPrevious++){\r
- AliKFParticle * previousGoodV0 = &vectorPreviousEventGoodV0s.at(iPrevious);\r
-\r
- AliKFParticle *backgroundCandidate = new AliKFParticle(*currentEventGoodV0,*previousGoodV0);\r
- \r
- Double_t massBG =0.;\r
- Double_t widthBG = 0.;\r
- Double_t chi2BG =10000.; \r
- backgroundCandidate->GetMass(massBG,widthBG);\r
- if(backgroundCandidate->GetNDF()>0){\r
- chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();\r
- if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){\r
- \r
- TVector3 MomentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());\r
- TVector3 SpaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());\r
- \r
- Double_t openingAngleBG = currentEventGoodV0->GetAngle(*previousGoodV0);\r
-\r
- Double_t rapidity;\r
- if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;\r
- else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));\r
-\r
- \r
- \r
- \r
- if(openingAngleBG < fMinOpeningAngleGhostCut ) continue; // minimum opening angle to avoid using ghosttracks\r
- \r
- \r
- fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);\r
- fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());\r
- fHistograms->FillHistogram("ESD_Background_Pt", MomentumVectorbackgroundCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Background_Eta", MomentumVectorbackgroundCandidate.Eta());\r
- fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);\r
- fHistograms->FillHistogram("ESD_Background_Phi", SpaceVectorbackgroundCandidate.Phi());\r
- fHistograms->FillHistogram("ESD_Background_Mass", massBG);\r
- fHistograms->FillHistogram("ESD_Background_R", SpaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!\r
- fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), SpaceVectorbackgroundCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());\r
- fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,MomentumVectorbackgroundCandidate.Pt());\r
- fHistograms->FillHistogram("ESD_Background_InvMass",massBG);\r
- }\r
- }\r
- delete backgroundCandidate; \r
- //cout << "nearly at the end of background" << endl;\r
- \r
- }\r
- }\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)\r
-{\r
- // Terminate analysis\r
- //\r
- AliDebug(1,"Do nothing in Terminate");\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()\r
-{\r
- // Create the output container\r
- if(fOutputContainer != NULL){\r
- delete fOutputContainer;\r
- fOutputContainer = NULL;\r
- }\r
- if(fOutputContainer == NULL){\r
- fOutputContainer = new TList();\r
- }\r
- \r
- //Adding the histograms to the output container\r
- fHistograms->GetOutputContainer(fOutputContainer);\r
- \r
- \r
- if(fWriteNtuple){\r
- if(fGammaNtuple == NULL){\r
- 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);\r
- }\r
- if(fNeutralMesonNtuple == NULL){\r
- fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");\r
- }\r
- TList * ntupleTList = new TList();\r
- ntupleTList->SetName("Ntuple");\r
- ntupleTList->Add((TNtuple*)fGammaNtuple);\r
- fOutputContainer->Add(ntupleTList);\r
- }\r
- \r
- fOutputContainer->SetName(GetName());\r
-}\r
-\r
-Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{\r
- //helper function\r
- TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());\r
- TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());\r
- return v3D0.Angle(v3D1);\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){\r
- // see header file for documentation\r
-\r
- vector<Int_t> indexOfGammaParticle;\r
-\r
- fStack = fV0Reader->GetMCStack();\r
-\r
- if(fV0Reader->CheckForPrimaryVertex() == kFALSE){\r
- return; // aborts if the primary vertex does not have contributors.\r
- }\r
-\r
- for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {\r
- TParticle* particle = (TParticle *)fStack->Particle(iTracks);\r
- if(particle->GetPdgCode()==22){ //Gamma\r
- if(particle->GetNDaughters() >= 2){\r
- TParticle* electron=NULL;\r
- TParticle* positron=NULL; \r
- for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){\r
- TParticle *tmpDaughter = fStack->Particle(daughterIndex);\r
- if(tmpDaughter->GetUniqueID() == 5){\r
- if(tmpDaughter->GetPdgCode() == 11){\r
- electron = tmpDaughter;\r
- }\r
- else if(tmpDaughter->GetPdgCode() == -11){\r
- positron = tmpDaughter;\r
- }\r
- }\r
- }\r
- if(electron!=NULL && positron!=0){\r
- if(electron->R()<160){\r
- indexOfGammaParticle.push_back(iTracks);\r
- }\r
- }\r
- }\r
- }\r
- }\r
-\r
- Int_t nFoundGammas=0;\r
- Int_t nNotFoundGammas=0;\r
-\r
- Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();\r
- for(Int_t i=0;i<numberOfV0s;i++){\r
- fV0Reader->GetV0(i);\r
- \r
- if(fV0Reader->HasSameMCMother() == kFALSE){\r
- continue;\r
- }\r
- \r
- TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();\r
- TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();\r
- \r
- if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){\r
- continue;\r
- }\r
- if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){\r
- continue;\r
- }\r
- \r
- if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){\r
- //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();\r
- for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){\r
- if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){\r
- nFoundGammas++;\r
- }\r
- else{\r
- nNotFoundGammas++;\r
- }\r
- }\r
- }\r
- }\r
- // cout<<"Found: "<<nFoundGammas<<" of: "<<indexOfGammaParticle.size()<<endl;\r
-}\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){\r
- // see header file for documantation\r
-\r
- fESDEvent = fV0Reader->GetESDEvent();\r
-\r
-\r
- vector <AliESDtrack*> vESDeNegTemp(0);\r
- vector <AliESDtrack*> vESDePosTemp(0);\r
- vector <AliESDtrack*> vESDxNegTemp(0);\r
- vector <AliESDtrack*> vESDxPosTemp(0);\r
- vector <AliESDtrack*> vESDeNegNoJPsi(0);\r
- vector <AliESDtrack*> vESDePosNoJPsi(0); \r
-\r
-\r
-\r
- fHistograms->FillTable("Table_Electrons",0);//Count number of Events\r
-\r
- for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){\r
- AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);\r
-\r
- if(!curTrack){\r
- //print warning here\r
- continue;\r
- }\r
-\r
- double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;\r
- double r[3];curTrack->GetConstrainedXYZ(r);\r
-\r
- TVector3 rXYZ(r);\r
-\r
- fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks\r
-\r
- Bool_t flagKink = kTRUE;\r
- Bool_t flagTPCrefit = kTRUE;\r
- Bool_t flagTRDrefit = kTRUE;\r
- Bool_t flagITSrefit = kTRUE;\r
- Bool_t flagTRDout = kTRUE;\r
- Bool_t flagVertex = kTRUE;\r
-\r
-\r
- //Cuts ---------------------------------------------------------------\r
-\r
- if(curTrack->GetKinkIndex(0) > 0){\r
- fHistograms->FillHistogram("Table_Electrons",5);//Count kink\r
- flagKink = kFALSE;\r
- }\r
-\r
- ULong_t trkStatus = curTrack->GetStatus();\r
-\r
- ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);\r
-\r
- if(!tpcRefit){\r
- fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit\r
- flagTPCrefit = kFALSE;\r
- }\r
-\r
- ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);\r
- if(!itsRefit){\r
- fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit\r
- flagITSrefit = kFALSE;\r
- }\r
-\r
- ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);\r
-\r
- if(!trdRefit){\r
- fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit\r
- flagTRDrefit = kFALSE;\r
- }\r
-\r
- ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);\r
-\r
- if(!trdOut) {\r
- fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout\r
- flagTRDout = kFALSE;\r
- }\r
-\r
- double nsigmaToVxt = GetSigmaToVertex(curTrack);\r
-\r
- if(nsigmaToVxt > 3){\r
- fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3\r
- flagVertex = kFALSE;\r
- }\r
-\r
- if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;\r
- fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts\r
-\r
-\r
- Stat_t pid, weight;\r
- GetPID(curTrack, pid, weight);\r
-\r
- if(pid!=0){\r
- fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0\r
- }\r
-\r
- if(pid == 0){\r
- fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0\r
- }\r
-\r
-\r
-\r
-\r
- Int_t labelMC = TMath::Abs(curTrack->GetLabel());\r
- TParticle* curParticle = fStack->Particle(labelMC);\r
-\r
-\r
-\r
-\r
- TLorentzVector curElec;\r
- curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);\r
-\r
-\r
-\r
-\r
- if(curTrack->GetSign() > 0){\r
-\r
- vESDxPosTemp.push_back(curTrack);\r
-\r
- if( pid == 0){\r
-\r
- fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
- fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());\r
- fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
- fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
- fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
- vESDePosTemp.push_back(curTrack);\r
-\r
-\r
-\r
- }\r
-\r
- }\r
- else {\r
- vESDxNegTemp.push_back(curTrack);\r
-\r
- if( pid == 0){\r
-\r
- fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());\r
- fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());\r
- fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());\r
- fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());\r
- fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());\r
- vESDeNegTemp.push_back(curTrack);\r
-\r
-\r
-\r
-\r
- }\r
- }\r
-\r
- }\r
-\r
-\r
- Bool_t ePosJPsi = kFALSE;\r
- Bool_t eNegJPsi = kFALSE; \r
- Bool_t ePosPi0 = kFALSE;\r
- Bool_t eNegPi0 = kFALSE;\r
- \r
- UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;\r
- \r
- for(UInt_t iNeg=0; iNeg < vESDeNegTemp.size(); iNeg++){\r
- if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetPdgCode() == 11)\r
- if(fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0) > -1){\r
- Int_t labelMother = fStack->Particle(TMath::Abs(vESDeNegTemp[iNeg]->GetLabel()))->GetMother(0);\r
- TParticle* partMother = fStack ->Particle(labelMother);\r
- if (partMother->GetPdgCode() == 111){\r
- ieNegPi0 = iNeg;\r
- eNegPi0 = kTRUE;\r
- }\r
- if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
- fHistograms->FillTable("Table_Electrons",14);\r
- ieNegJPsi = iNeg;\r
- eNegJPsi = kTRUE;\r
- }\r
- else{ \r
- vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);\r
- // cout<<"ESD No Positivo JPsi "<<endl;\r
- }\r
-\r
- }\r
- } \r
-\r
- for(UInt_t iPos=0; iPos < vESDePosTemp.size(); iPos++){\r
- if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetPdgCode() == -11)\r
- if(fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0) > -1){\r
- Int_t labelMother = fStack->Particle(TMath::Abs(vESDePosTemp[iPos]->GetLabel()))->GetMother(0);\r
- TParticle* partMother = fStack ->Particle(labelMother);\r
- if (partMother->GetPdgCode() == 111){\r
- iePosPi0 = iPos;\r
- ePosPi0 = kTRUE;\r
- }\r
- if(partMother->GetPdgCode() == 443){ //Mother JPsi\r
- fHistograms->FillTable("Table_Electrons",15);\r
- iePosJPsi = iPos;\r
- ePosJPsi = kTRUE;\r
- }\r
- else{\r
- vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);\r
- // cout<<"ESD No Negativo JPsi "<<endl;\r
- }\r
-\r
- }\r
- }\r
- \r
- if( eNegJPsi && ePosJPsi ){\r
- TVector3 tempeNegV,tempePosV;\r
- tempeNegV.SetXYZ(vESDeNegTemp[ieNegJPsi]->Px(),vESDeNegTemp[ieNegJPsi]->Py(),vESDeNegTemp[ieNegJPsi]->Pz()); \r
- tempePosV.SetXYZ(vESDePosTemp[iePosJPsi]->Px(),vESDePosTemp[iePosJPsi]->Py(),vESDePosTemp[iePosJPsi]->Pz());\r
- fHistograms->FillTable("Table_Electrons",16);\r
- fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV)); \r
- fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegJPsi]->GetLabel())),\r
- fStack->Particle(TMath::Abs(vESDePosTemp[iePosJPsi]->GetLabel())))); \r
- }\r
- \r
- if( eNegPi0 && ePosPi0 ){\r
- TVector3 tempeNegV,tempePosV;\r
- tempeNegV.SetXYZ(vESDeNegTemp[ieNegPi0]->Px(),vESDeNegTemp[ieNegPi0]->Py(),vESDeNegTemp[ieNegPi0]->Pz());\r
- tempePosV.SetXYZ(vESDePosTemp[iePosPi0]->Px(),vESDePosTemp[iePosPi0]->Py(),vESDePosTemp[iePosPi0]->Pz());\r
- fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));\r
- fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(vESDeNegTemp[ieNegPi0]->GetLabel())),\r
- fStack->Particle(TMath::Abs(vESDePosTemp[iePosPi0]->GetLabel())))); \r
- }\r
- \r
-\r
- FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));\r
-\r
- CleanWithAngleCuts(vESDeNegTemp,vESDePosTemp,fKFReconstructedGammas);\r
- \r
- vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);\r
- vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);\r
-\r
-\r
- FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);\r
-\r
- \r
-\r
-\r
- //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);\r
-\r
-\r
- FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);\r
- FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));\r
-\r
- \r
-\r
- FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",\r
- fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
-\r
- FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",\r
- fKFReconstructedGammasCut,vCurrentTLVeNeg,vCurrentTLVePos);\r
-\r
- //BackGround\r
-\r
- //Like Sign e+e-\r
- ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);\r
- ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);\r
- ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);\r
- ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);\r
-\r
- // Like Sign e+e- no JPsi\r
- ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));\r
- ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));\r
-\r
- //Mixed Event\r
-\r
- if( fCurrentEventPosElectron.size() > 0 && fCurrentEventNegElectron.size() > 0 && fKFReconstructedGammasCut.size() > 0 ){\r
- FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",\r
- fKFReconstructedGammasCut,fPreviousEventTLVNegElectron,fPreviousEventTLVPosElectron);\r
- fPreviousEventTLVNegElectron = vCurrentTLVeNeg;\r
- fPreviousEventTLVPosElectron = vCurrentTLVePos;\r
-\r
- }\r
-\r
- /*\r
- //Photons P\r
- Double_t vtx[3];\r
- vtx[0]=0;vtx[1]=0;vtx[2]=0;\r
- for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){\r
-\r
- // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;\r
-\r
-\r
-\r
- Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();\r
- // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;\r
-\r
- // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;\r
-\r
- if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )\r
-\r
- fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());\r
-\r
-\r
- }\r
-\r
-\r
- */\r
-\r
-\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){\r
- //see header file for documentation\r
- for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){\r
- for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
- fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));\r
- }\r
- }\r
-}\r
-void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, vector <TLorentzVector> eNeg, vector <TLorentzVector> ePos){\r
- //see header file for documentation\r
- for( UInt_t n=0; n < eNeg.size(); n++){\r
-\r
- TLorentzVector en = eNeg.at(n);\r
- for (UInt_t p=0; p < ePos.size(); p++){\r
- TLorentzVector ep = ePos.at(p);\r
- TLorentzVector np = ep + en;\r
- fHistograms->FillHistogram(histoName.Data(),np.M());\r
- }\r
- }\r
-\r
-}\r
-\r
-void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,vector <AliKFParticle> fKFGammas,\r
- vector <TLorentzVector> tlVeNeg,vector<TLorentzVector> tlVePos)\r
-{\r
- //see header file for documentation\r
-\r
- for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++ ){\r
-\r
- for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){\r
-\r
- TLorentzVector xy = tlVePos[iPos] + tlVeNeg[iNeg];\r
-\r
- for (UInt_t iGam=0; iGam < fKFGammas.size(); iGam++){\r
-\r
- AliKFParticle * gammaCandidate = &fKFGammas[iGam];\r
- TLorentzVector g;\r
-\r
- g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);\r
- TLorentzVector xyg = xy + g;\r
- fHistograms->FillHistogram(histoMass.Data(),xyg.M());\r
- fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));\r
- }\r
- }\r
- }\r
-\r
-}\r
-void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, vector <TLorentzVector> e)\r
-{\r
- // see header file for documentation\r
- for(UInt_t i=0; i < e.size(); i++)\r
- {\r
- for (UInt_t j=i+1; j < e.size(); j++)\r
- {\r
- TLorentzVector ee = e[i] + e[j];\r
-\r
- fHistograms->FillHistogram(hBg.Data(),ee.M());\r
- }\r
- }\r
-}\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(vector <AliESDtrack*> negativeElectrons,\r
- vector <AliESDtrack*> positiveElectrons, vector <AliKFParticle> gammas){\r
- // see header file for documentation\r
-\r
- UInt_t sizeN = negativeElectrons.size();\r
- UInt_t sizeP = positiveElectrons.size();\r
- UInt_t sizeG = gammas.size();\r
-\r
-\r
-\r
- vector <Bool_t> xNegBand(sizeN);\r
- vector <Bool_t> xPosBand(sizeP);\r
- vector <Bool_t> gammaBand(sizeG);\r
-\r
-\r
- for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;\r
- for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;\r
- for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;\r
-\r
-\r
- for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
- \r
- Double_t aP[3]; positiveElectrons[iPos]->GetConstrainedPxPyPz(aP); \r
-\r
- TVector3 ePosV(aP[0],aP[1],aP[2]);\r
-\r
- for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){\r
- \r
- Double_t aN[3]; negativeElectrons[iNeg]->GetConstrainedPxPyPz(aN); \r
- TVector3 eNegV(aN[0],aN[1],aN[2]);\r
-\r
- if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma\r
- xPosBand[iPos]=kFALSE;\r
- xNegBand[iNeg]=kFALSE;\r
- }\r
-\r
- for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
- AliKFParticle* gammaCandidate = &gammas[iGam];\r
- TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());\r
- if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)\r
- gammaBand[iGam]=kFALSE;\r
- }\r
- }\r
- }\r
-\r
-\r
-\r
-\r
- for(UInt_t iPos=0; iPos < sizeP; iPos++){\r
- if(xPosBand[iPos]){\r
- fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);\r
- }\r
- }\r
- for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){\r
- if(xNegBand[iNeg]){\r
- fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);\r
- }\r
- }\r
- for(UInt_t iGam=0; iGam < sizeG; iGam++){\r
- if(gammaBand[iGam]){\r
- fKFReconstructedGammasCut.push_back(gammas[iGam]);\r
- }\r
- }\r
-}\r
-\r
-\r
-void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)\r
-{\r
- // see header file for documentation\r
- pid = -1;\r
- weight = -1;\r
-\r
- double wpart[5];\r
- double wpartbayes[5];\r
-\r
- //get probability of the diffenrent particle types\r
- track->GetESDpid(wpart);\r
-\r
- // Tentative particle type "concentrations"\r
- double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};\r
-\r
- //Bayes' formula\r
- double rcc = 0.;\r
- for (int i = 0; i < 5; i++)\r
- {\r
- rcc+=(c[i] * wpart[i]);\r
- }\r
-\r
-\r
-\r
- for (int i=0; i<5; i++) {\r
- if( rcc!=0){\r
- wpartbayes[i] = c[i] * wpart[i] / rcc;\r
- }\r
- }\r
-\r
-\r
-\r
- Float_t max=0.;\r
- int ipid=-1;\r
- //find most probable particle in ESD pid\r
- //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton\r
- for (int i = 0; i < 5; i++)\r
- {\r
- if (wpartbayes[i] > max)\r
- {\r
- ipid = i;\r
- max = wpartbayes[i];\r
- }\r
- }\r
-\r
- pid = ipid;\r
- weight = max;\r
-}\r
-double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)\r
-{\r
- // Calculates the number of sigma to the vertex.\r
-\r
- Float_t b[2];\r
- Float_t bRes[2];\r
- Float_t bCov[3];\r
- t->GetImpactParameters(b,bCov);\r
- if (bCov[0]<=0 || bCov[2]<=0) {\r
- AliDebug(1, "Estimated b resolution lower or equal zero!");\r
- bCov[0]=0; bCov[2]=0;\r
- }\r
- bRes[0] = TMath::Sqrt(bCov[0]);\r
- bRes[1] = TMath::Sqrt(bCov[2]);\r
-\r
- // -----------------------------------\r
- // How to get to a n-sigma cut?\r
- //\r
- // The accumulated statistics from 0 to d is\r
- //\r
- // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)\r
- // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)\r
- //\r
- // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)\r
- // Can this be expressed in a different way?\r
-\r
- if (bRes[0] == 0 || bRes[1] ==0)\r
- return -1;\r
-\r
- double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));\r
-\r
- // stupid rounding problem screws up everything:\r
- // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(\r
- if (TMath::Exp(-d * d / 2) < 1e-10)\r
- return 1000;\r
-\r
-\r
- d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);\r
- return d;\r
-}\r
-vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){\r
-\r
- vector <TLorentzVector> tlVtrack(0);\r
-\r
- for(UInt_t itrack=0; itrack < esdTrack.size(); itrack++){\r
- double P[3]; esdTrack[itrack]->GetConstrainedPxPyPz(P);\r
- TLorentzVector currentTrack;\r
- currentTrack.SetXYZM(P[0],P[1],P[2],fElectronMass);\r
- tlVtrack.push_back(currentTrack);\r
- }\r
-\r
- return tlVtrack;\r
-}\r
-\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: Ana Marin, Kathrin Koch, Kenneth Aamodt *
+ * Version 1.1 *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+////////////////////////////////////////////////
+//---------------------------------------------
+// Class used to do analysis on conversion pairs
+//---------------------------------------------
+////////////////////////////////////////////////
+
+// root
+#include <TChain.h>
+
+// analysis
+#include "AliAnalysisTaskGammaConversion.h"
+#include "AliStack.h"
+#include "AliLog.h"
+#include "AliESDtrackCuts.h"
+#include "TNtuple.h"
+//#include "AliCFManager.h" // for CF
+//#include "AliCFContainer.h" // for CF
+#include "AliGammaConversionAODObject.h"
+#include "AliGammaConversionBGHandler.h"
+
+class AliCFContainer;
+class AliCFManager;
+class AliKFVertex;
+class AliAODHandler;
+class AliAODEvent;
+class ALiESDEvent;
+class AliMCEvent;
+class AliMCEventHandler;
+class AliESDInputHandler;
+class AliAnalysisManager;
+class Riostream;
+class TFile;
+class TInterpreter;
+class TSystem;
+class TROOT;
+
+ClassImp(AliAnalysisTaskGammaConversion)
+
+
+AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion():
+AliAnalysisTaskSE(),
+ fV0Reader(NULL),
+ fStack(NULL),
+ fMCTruth(NULL), // for CF
+ fGCMCEvent(NULL), // for CF
+ fESDEvent(NULL),
+ fOutputContainer(NULL),
+ fCFManager(0x0), // for CF
+ fHistograms(NULL),
+ fTriggerCINT1B(kFALSE),
+ fDoMCTruth(kFALSE),
+ fDoNeutralMeson(kFALSE),
+ fDoJet(kFALSE),
+ fDoChic(kFALSE),
+ fKFReconstructedGammasTClone(NULL),
+ fCurrentEventPosElectronTClone(NULL),
+ fCurrentEventNegElectronTClone(NULL),
+ fKFReconstructedGammasCutTClone(NULL),
+ fPreviousEventTLVNegElectronTClone(NULL),
+ fPreviousEventTLVPosElectronTClone(NULL),
+// fKFReconstructedGammas(),
+ fElectronv1(),
+ fElectronv2(),
+// fCurrentEventPosElectron(),
+// fCurrentEventNegElectron(),
+// fKFReconstructedGammasCut(),
+// fPreviousEventTLVNegElectron(),
+// fPreviousEventTLVPosElectron(),
+ fElectronMass(-1),
+ fGammaMass(-1),
+ fPi0Mass(-1),
+ fEtaMass(-1),
+ fGammaWidth(-1),
+ fPi0Width(-1),
+ fEtaWidth(-1),
+ fMinOpeningAngleGhostCut(0.),
+ fEsdTrackCuts(NULL),
+ fCalculateBackground(kFALSE),
+ fWriteNtuple(kFALSE),
+ fGammaNtuple(NULL),
+ fNeutralMesonNtuple(NULL),
+ fTotalNumberOfAddedNtupleEntries(0),
+ fChargedParticles(NULL),
+ fChargedParticlesId(),
+ fGammaPtHighest(0.),
+ fMinPtForGammaJet(1.),
+ fMinIsoConeSize(0.2),
+ fMinPtIsoCone(0.7),
+ fMinPtGamChargedCorr(0.5),
+ fMinPtJetCone(0.5),
+ fLeadingChargedIndex(-1),
+ fLowPtMapping(1.),
+ fHighPtMapping(3.),
+ fDoCF(kFALSE),
+ fAODBranch(NULL),
+ fAODBranchName("GammaConv")//,
+ // fAODObjects(NULL)
+{
+ // Default constructor
+
+ /* Kenneth: the default constructor should not have any define input/output or the call to SetESDtrackCuts
+ // Common I/O in slot 0
+ DefineInput (0, TChain::Class());
+ DefineOutput(0, TTree::Class());
+
+ // Your private output
+ DefineOutput(1, TList::Class());
+
+ // Define standard ESD track cuts for Gamma-hadron correlation
+ SetESDtrackCuts();
+ */
+}
+
+AliAnalysisTaskGammaConversion::AliAnalysisTaskGammaConversion(const char* name):
+ AliAnalysisTaskSE(name),
+ fV0Reader(NULL),
+ fStack(NULL),
+ fMCTruth(NULL), // for CF
+ fGCMCEvent(NULL), // for CF
+ fESDEvent(NULL),
+ fOutputContainer(0x0),
+ fCFManager(0x0), // for CF
+ fHistograms(NULL),
+ fTriggerCINT1B(kFALSE),
+ fDoMCTruth(kFALSE),
+ fDoNeutralMeson(kFALSE),
+ fDoJet(kFALSE),
+ fDoChic(kFALSE),
+ fKFReconstructedGammasTClone(NULL),
+ fCurrentEventPosElectronTClone(NULL),
+ fCurrentEventNegElectronTClone(NULL),
+ fKFReconstructedGammasCutTClone(NULL),
+ fPreviousEventTLVNegElectronTClone(NULL),
+ fPreviousEventTLVPosElectronTClone(NULL),
+ // fKFReconstructedGammas(),
+ fElectronv1(),
+ fElectronv2(),
+ // fCurrentEventPosElectron(),
+ // fCurrentEventNegElectron(),
+ // fKFReconstructedGammasCut(),
+ // fPreviousEventTLVNegElectron(),
+ // fPreviousEventTLVPosElectron(),
+ fElectronMass(-1),
+ fGammaMass(-1),
+ fPi0Mass(-1),
+ fEtaMass(-1),
+ fGammaWidth(-1),
+ fPi0Width(-1),
+ fEtaWidth(-1),
+ fMinOpeningAngleGhostCut(0.),
+ fEsdTrackCuts(NULL),
+ fCalculateBackground(kFALSE),
+ fWriteNtuple(kFALSE),
+ fGammaNtuple(NULL),
+ fNeutralMesonNtuple(NULL),
+ fTotalNumberOfAddedNtupleEntries(0),
+ fChargedParticles(NULL),
+ fChargedParticlesId(),
+ fGammaPtHighest(0.),
+ fMinPtForGammaJet(1.),
+ fMinIsoConeSize(0.2),
+ fMinPtIsoCone(0.7),
+ fMinPtGamChargedCorr(0.5),
+ fMinPtJetCone(0.5),
+ fLeadingChargedIndex(-1),
+ fLowPtMapping(1.),
+ fHighPtMapping(3.),
+ fDoCF(kFALSE),
+ fAODBranch(NULL),
+ fAODBranchName("GammaConv")//,
+ // fAODObjects(NULL)
+{
+ // Common I/O in slot 0
+ DefineInput (0, TChain::Class());
+ DefineOutput(0, TTree::Class());
+
+ // Your private output
+ DefineOutput(1, TList::Class());
+ DefineOutput(2, AliCFContainer::Class()); // for CF
+
+
+ // Define standard ESD track cuts for Gamma-hadron correlation
+ SetESDtrackCuts();
+}
+
+AliAnalysisTaskGammaConversion::~AliAnalysisTaskGammaConversion()
+{
+ // Remove all pointers
+
+ if(fOutputContainer){
+ fOutputContainer->Clear() ;
+ delete fOutputContainer ;
+ }
+ if(fHistograms){
+ delete fHistograms;
+ }
+ if(fV0Reader){
+ delete fV0Reader;
+ }
+
+ // for CF
+ if(fCFManager){
+ delete fCFManager;
+ }
+
+ if(fEsdTrackCuts){
+ delete fEsdTrackCuts;
+ }
+
+ if (fAODBranch) {
+ fAODBranch->Clear();
+ delete fAODBranch ;
+ }
+}
+
+
+void AliAnalysisTaskGammaConversion::Init()
+{
+ // Initialization
+ // AliLog::SetGlobalLogLevel(AliLog::kError);
+}
+void AliAnalysisTaskGammaConversion::SetESDtrackCuts()
+{
+ // SetESDtrackCuts
+ if (fEsdTrackCuts!=NULL){
+ delete fEsdTrackCuts;
+ }
+ fEsdTrackCuts = new AliESDtrackCuts("AliESDtrackCuts");
+ //standard cuts from:
+ //http://aliceinfo.cern.ch/alicvs/viewvc/PWG0/dNdEta/CreateCuts.C?revision=1.4&view=markup
+ //fEsdTrackCuts->SetMinNClustersTPC(50);
+ //fEsdTrackCuts->SetMaxChi2PerClusterTPC(3.5);
+ //fEsdTrackCuts->SetMaxCovDiagonalElements(2,2,0.5,0.5,2);
+ fEsdTrackCuts->SetRequireTPCRefit(kTRUE);
+ fEsdTrackCuts->SetRequireITSRefit(kTRUE);
+ fEsdTrackCuts->SetMaxNsigmaToVertex(3);
+ fEsdTrackCuts->SetRequireSigmaToVertex(kTRUE);
+ // fEsdTrackCuts->SetAcceptKinkDaughters(kFALSE);
+
+}
+
+void AliAnalysisTaskGammaConversion::Exec(Option_t */*option*/)
+{
+ // Execute analysis for current event
+ ConnectInputData("");
+
+ //Each event needs an empty branch
+ // fAODBranch->Clear();
+ fAODBranch->Delete();
+
+ if(fKFReconstructedGammasTClone == NULL){
+ fKFReconstructedGammasTClone = new TClonesArray("AliKFParticle",0);
+ }
+ if(fCurrentEventPosElectronTClone == NULL){
+ fCurrentEventPosElectronTClone = new TClonesArray("AliESDtrack",0);
+ }
+ if(fCurrentEventNegElectronTClone == NULL){
+ fCurrentEventNegElectronTClone = new TClonesArray("AliESDtrack",0);
+ }
+ if(fKFReconstructedGammasCutTClone == NULL){
+ fKFReconstructedGammasCutTClone = new TClonesArray("AliKFParticle",0);
+ }
+ if(fPreviousEventTLVNegElectronTClone == NULL){
+ fPreviousEventTLVNegElectronTClone = new TClonesArray("TLorentzVector",0);
+ }
+ if(fPreviousEventTLVPosElectronTClone == NULL){
+ fPreviousEventTLVPosElectronTClone = new TClonesArray("TLorentzVector",0);
+ }
+ if(fChargedParticles == NULL){
+ fChargedParticles = new TClonesArray("AliESDtrack",0);
+ }
+
+ //clear TClones
+ fKFReconstructedGammasTClone->Delete();
+ fCurrentEventPosElectronTClone->Delete();
+ fCurrentEventNegElectronTClone->Delete();
+ fKFReconstructedGammasCutTClone->Delete();
+ fPreviousEventTLVNegElectronTClone->Delete();
+ fPreviousEventTLVPosElectronTClone->Delete();
+
+ //clear vectors
+ fElectronv1.clear();
+ fElectronv2.clear();
+
+ fChargedParticles->Delete();
+
+ fChargedParticlesId.clear();
+
+ //Clear the data in the v0Reader
+ // fV0Reader->UpdateEventByEventData();
+
+ //Take Only events with proper trigger
+ /*
+ if(fTriggerCINT1B){
+ if(!fV0Reader->GetESDEvent()->IsTriggerClassFired("CINT1B-ABCE-NOPF-ALL")) return;
+ }
+ */
+
+ // Process the MC information
+ if(fDoMCTruth){
+ ProcessMCData();
+ }
+
+ //Process the v0 information with no cuts
+ ProcessV0sNoCut();
+
+ // Process the v0 information
+ ProcessV0s();
+
+
+ //Fill Gamma AOD
+ FillAODWithConversionGammas() ;
+
+ // Process reconstructed gammas
+ if(fDoNeutralMeson == kTRUE){
+ ProcessGammasForNeutralMesonAnalysis();
+ }
+
+ if(fDoMCTruth == kTRUE){
+ CheckV0Efficiency();
+ }
+ //Process reconstructed gammas electrons for Chi_c Analysis
+ if(fDoChic == kTRUE){
+ ProcessGammaElectronsForChicAnalysis();
+ }
+ // Process reconstructed gammas for gamma Jet/hadron correlations
+ if(fDoJet == kTRUE){
+ ProcessGammasForGammaJetAnalysis();
+ }
+
+ //calculate background if flag is set
+ if(fCalculateBackground){
+ CalculateBackground();
+ }
+
+ //Clear the data in the v0Reader
+ fV0Reader->UpdateEventByEventData();
+
+ PostData(1, fOutputContainer);
+ PostData(2, fCFManager->GetParticleContainer()); // for CF
+
+}
+
+void AliAnalysisTaskGammaConversion::ConnectInputData(Option_t *option){
+ // see header file for documentation
+
+ AliAnalysisTaskSE::ConnectInputData(option);
+
+ if(fV0Reader == NULL){
+ // Write warning here cuts and so on are default if this ever happens
+ }
+ fV0Reader->Initialize();
+ fDoMCTruth = fV0Reader->GetDoMCTruth();
+}
+
+
+
+void AliAnalysisTaskGammaConversion::ProcessMCData(){
+ // see header file for documentation
+
+ fStack = fV0Reader->GetMCStack();
+ fMCTruth = fV0Reader->GetMCTruth(); // for CF
+ fGCMCEvent = fV0Reader->GetMCEvent(); // for CF
+
+
+ // for CF
+ Double_t containerInput[3];
+ if(fDoCF){
+ if(!fGCMCEvent) cout << "NO MC INFO FOUND" << endl;
+ fCFManager->SetEventInfo(fGCMCEvent);
+ }
+ // end for CF
+
+
+ if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
+ return; // aborts if the primary vertex does not have contributors.
+ }
+
+ for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++) {
+ TParticle* particle = (TParticle *)fStack->Particle(iTracks);
+
+ if (!particle) {
+ //print warning here
+ continue;
+ }
+
+ ///////////////////////Begin Chic Analysis/////////////////////////////
+
+ if(particle->GetPdgCode() == 443){//Is JPsi
+ if(particle->GetNDaughters()==2){
+ if(TMath::Abs(fStack->Particle(particle->GetFirstDaughter())->GetPdgCode()) == 11 &&
+ TMath::Abs(fStack->Particle(particle->GetLastDaughter())->GetPdgCode()) == 11){
+ TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
+ TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
+ if(TMath::Abs(daug0->Eta()) < 0.9 && TMath::Abs(daug1->Eta()) < 0.9)
+ fHistograms->FillTable("Table_Electrons",3);//e+ e- from J/Psi inside acceptance
+
+ if( TMath::Abs(daug0->Eta()) < 0.9){
+ if(daug0->GetPdgCode() == -11)
+ fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
+ else
+ fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
+
+ }
+ if(TMath::Abs(daug1->Eta()) < 0.9){
+ if(daug1->GetPdgCode() == -11)
+ fHistograms->FillTable("Table_Electrons",1);//e+ from J/Psi inside acceptance
+ else
+ fHistograms->FillTable("Table_Electrons",2);//e- from J/Psi inside acceptance
+ }
+ }
+ }
+ }
+ // const int CHI_C0 = 10441;
+ // const int CHI_C1 = 20443;
+ // const int CHI_C2 = 445
+ if(particle->GetPdgCode() == 22){//gamma from JPsi
+ if(particle->GetMother(0) > -1){
+ if(fStack->Particle(particle->GetMother(0))->GetPdgCode() == 10441 ||
+ fStack->Particle(particle->GetMother(0))->GetPdgCode() == 20443 ||
+ fStack->Particle(particle->GetMother(0))->GetPdgCode() == 445){
+ if(TMath::Abs(particle->Eta()) < 1.2)
+ fHistograms->FillTable("Table_Electrons",17);// gamma from chic inside accptance
+ }
+ }
+ }
+ if(particle->GetPdgCode() == 10441 || particle->GetPdgCode() == 20443 || particle->GetPdgCode() == 445){
+ if( particle->GetNDaughters() == 2){
+ TParticle* daug0 = fStack->Particle(particle->GetFirstDaughter());
+ TParticle* daug1 = fStack->Particle(particle->GetLastDaughter());
+
+ if( (daug0->GetPdgCode() == 443 || daug0->GetPdgCode() == 22) && (daug1->GetPdgCode() == 443 || daug1->GetPdgCode() == 22) ){
+ if( daug0->GetPdgCode() == 443){
+ TParticle* daugE0 = fStack->Particle(daug0->GetFirstDaughter());
+ TParticle* daugE1 = fStack->Particle(daug0->GetLastDaughter());
+ if( TMath::Abs(daug1->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
+ fHistograms->FillTable("Table_Electrons",18);
+
+ }//if
+ else if (daug1->GetPdgCode() == 443){
+ TParticle* daugE0 = fStack->Particle(daug1->GetFirstDaughter());
+ TParticle* daugE1 = fStack->Particle(daug1->GetLastDaughter());
+ if( TMath::Abs(daug0->Eta()) < 1.2 && TMath::Abs(daugE0->Eta()) < 0.9 && TMath::Abs(daugE1->Eta()) < 0.9 )
+ fHistograms->FillTable("Table_Electrons",18);
+ }//else if
+ }//gamma o Jpsi
+ }//GetNDaughters
+ }
+
+
+ /////////////////////End Chic Analysis////////////////////////////
+
+
+ if(TMath::Abs(particle->Eta())> fV0Reader->GetEtaCut() ) continue;
+
+ if(particle->R()>fV0Reader->GetMaxRCut()) continue; // cuts on distance from collision point
+
+ Double_t tmpPhi=particle->Phi();
+
+ if(particle->Phi()> TMath::Pi()){
+ tmpPhi = particle->Phi()-(2*TMath::Pi());
+ }
+
+ Double_t rapidity;
+ if(particle->Energy() - particle->Pz() == 0 || particle->Energy() + particle->Pz() == 0){
+ rapidity=0;
+ }
+ else{
+ rapidity = 0.5*(TMath::Log((particle->Energy()+particle->Pz()) / (particle->Energy()-particle->Pz())));
+ }
+
+ //process the gammas
+ if (particle->GetPdgCode() == 22){
+
+
+ if(particle->GetMother(0) >-1 && fStack->Particle(particle->GetMother(0))->GetPdgCode() == 22){
+ continue; // no photon as mothers!
+ }
+
+ if(particle->GetMother(0) >= fStack->GetNprimary()){
+ continue; // the gamma has a mother, and it is not a primary particle
+ }
+
+ fHistograms->FillHistogram("MC_allGamma_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_allGamma_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_allGamma_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_allGamma_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_allGamma_Rapid", rapidity);
+
+ // for CF
+ if(fDoCF){
+ containerInput[0] = particle->Pt();
+ containerInput[1] = particle->Eta();
+ if(particle->GetMother(0) >=0){
+ containerInput[2] = fStack->Particle(particle->GetMother(0))->GetMass();
+ }
+ else{
+ containerInput[2]=-1;
+ }
+
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepGenerated); // generated gamma
+ }
+ if(particle->GetMother(0) < 0){ // direct gamma
+ fHistograms->FillHistogram("MC_allDirectGamma_Energy",particle->Energy());
+ fHistograms->FillHistogram("MC_allDirectGamma_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_allDirectGamma_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_allDirectGamma_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_allDirectGamma_Rapid", rapidity);
+ }
+
+ // looking for conversion (electron + positron from pairbuilding (= 5) )
+ TParticle* ePos = NULL;
+ TParticle* eNeg = NULL;
+
+ if(particle->GetNDaughters() >= 2){
+ for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
+ TParticle *tmpDaughter = fStack->Particle(daughterIndex);
+ if(tmpDaughter->GetUniqueID() == 5){
+ if(tmpDaughter->GetPdgCode() == 11){
+ eNeg = tmpDaughter;
+ }
+ else if(tmpDaughter->GetPdgCode() == -11){
+ ePos = tmpDaughter;
+ }
+ }
+ }
+ }
+
+
+ if(ePos == NULL || eNeg == NULL){ // means we do not have two daughters from pair production
+ continue;
+ }
+
+
+ Double_t ePosPhi = ePos->Phi();
+ if(ePos->Phi()> TMath::Pi()) ePosPhi = ePos->Phi()-(2*TMath::Pi());
+
+ Double_t eNegPhi = eNeg->Phi();
+ if(eNeg->Phi()> TMath::Pi()) eNegPhi = eNeg->Phi()-(2*TMath::Pi());
+
+
+ if(ePos->Pt()<fV0Reader->GetPtCut() || eNeg->Pt()<fV0Reader->GetPtCut()){
+ continue; // no reconstruction below the Pt cut
+ }
+
+ if(TMath::Abs(ePos->Eta())> fV0Reader->GetEtaCut() || TMath::Abs(eNeg->Eta())> fV0Reader->GetEtaCut()){
+ continue;
+ }
+
+ if(ePos->R()>fV0Reader->GetMaxRCut()){
+ continue; // cuts on distance from collision point
+ }
+
+ if(TMath::Abs(ePos->Vz()) > fV0Reader->GetMaxZCut()){
+ continue; // outside material
+ }
+
+
+ if((TMath::Abs(ePos->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() > ePos->R()){
+ continue; // line cut to exclude regions where we do not reconstruct
+ }
+
+
+ // for CF
+ if(fDoCF){
+ fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructable); // reconstructable gamma
+ }
+ fHistograms->FillHistogram("MC_ConvGamma_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_ConvGamma_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_ConvGamma_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_ConvGamma_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_ConvGamma_Rapid", rapidity);
+ fHistograms->FillHistogram("MC_ConvGamma_Pt_Eta", particle->Pt(),particle->Eta());
+
+ fHistograms->FillHistogram("MC_E_Energy", eNeg->Energy());
+ fHistograms->FillHistogram("MC_E_Pt", eNeg->Pt());
+ fHistograms->FillHistogram("MC_E_Eta", eNeg->Eta());
+ fHistograms->FillHistogram("MC_E_Phi", eNegPhi);
+
+ fHistograms->FillHistogram("MC_P_Energy", ePos->Energy());
+ fHistograms->FillHistogram("MC_P_Pt", ePos->Pt());
+ fHistograms->FillHistogram("MC_P_Eta", ePos->Eta());
+ fHistograms->FillHistogram("MC_P_Phi", ePosPhi);
+
+
+ // begin Mapping
+ Int_t rBin = fHistograms->GetRBin(ePos->R());
+ Int_t zBin = fHistograms->GetZBin(ePos->Vz());
+ Int_t phiBin = fHistograms->GetPhiBin(particle->Phi());
+
+ TVector3 vtxPos(ePos->Vx(),ePos->Vy(),ePos->Vz());
+
+ TString nameMCMappingPhiR="";
+ nameMCMappingPhiR.Form("MC_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
+ fHistograms->FillHistogram(nameMCMappingPhiR, ePos->Vz(), particle->Eta());
+
+ TString nameMCMappingPhi="";
+ nameMCMappingPhi.Form("MC_Conversion_Mapping_Phi%02d",phiBin);
+ // fHistograms->FillHistogram(nameMCMappingPhi, particle->Eta());
+ fHistograms->FillHistogram(nameMCMappingPhi, ePos->Vz(), particle->Eta());
+
+ TString nameMCMappingR="";
+ nameMCMappingR.Form("MC_Conversion_Mapping_R%02d",rBin);
+ // fHistograms->FillHistogram(nameMCMappingR, particle->Eta());
+ fHistograms->FillHistogram(nameMCMappingR,ePos->Vz(), particle->Eta());
+
+ TString nameMCMappingPhiInR="";
+ nameMCMappingPhiInR.Form("MC_Conversion_Mapping_Phi_in_R_%02d",rBin);
+ // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
+ fHistograms->FillHistogram(nameMCMappingPhiInR, vtxPos.Phi());
+
+ TString nameMCMappingZInR="";
+ nameMCMappingZInR.Form("MC_Conversion_Mapping_Z_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameMCMappingZInR,ePos->Vz() );
+
+
+ TString nameMCMappingPhiInZ="";
+ nameMCMappingPhiInZ.Form("MC_Conversion_Mapping_Phi_in_Z_%02d",zBin);
+ // fHistograms->FillHistogram(nameMCMappingPhiInR, tmpPhi);
+ fHistograms->FillHistogram(nameMCMappingPhiInZ, vtxPos.Phi());
+
+ TString nameMCMappingRInZ="";
+ nameMCMappingRInZ.Form("MC_Conversion_Mapping_R_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameMCMappingRInZ,ePos->R() );
+
+ if(particle->Pt() > fLowPtMapping && particle->Pt()< fHighPtMapping){
+ TString nameMCMappingMidPtPhiInR="";
+ nameMCMappingMidPtPhiInR.Form("MC_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameMCMappingMidPtPhiInR, vtxPos.Phi());
+
+ TString nameMCMappingMidPtZInR="";
+ nameMCMappingMidPtZInR.Form("MC_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameMCMappingMidPtZInR,ePos->Vz() );
+
+
+ TString nameMCMappingMidPtPhiInZ="";
+ nameMCMappingMidPtPhiInZ.Form("MC_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameMCMappingMidPtPhiInZ, vtxPos.Phi());
+
+ TString nameMCMappingMidPtRInZ="";
+ nameMCMappingMidPtRInZ.Form("MC_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameMCMappingMidPtRInZ,ePos->R() );
+
+ }
+
+ //end mapping
+
+ fHistograms->FillHistogram("MC_Conversion_R",ePos->R());
+ fHistograms->FillHistogram("MC_Conversion_ZR",ePos->Vz(),ePos->R());
+ fHistograms->FillHistogram("MC_Conversion_XY",ePos->Vx(),ePos->Vy());
+ fHistograms->FillHistogram("MC_Conversion_OpeningAngle",GetMCOpeningAngle(ePos, eNeg));
+ fHistograms->FillHistogram("MC_ConvGamma_E_AsymmetryP",particle->P(),eNeg->P()/particle->P());
+ fHistograms->FillHistogram("MC_ConvGamma_P_AsymmetryP",particle->P(),ePos->P()/particle->P());
+
+
+ if(particle->GetMother(0) < 0){ // no mother = direct gamma, still inside converted
+ fHistograms->FillHistogram("MC_ConvDirectGamma_Energy",particle->Energy());
+ fHistograms->FillHistogram("MC_ConvDirectGamma_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_ConvDirectGamma_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_ConvDirectGamma_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_ConvDirectGamma_Rapid", rapidity);
+
+ } // end direct gamma
+ else{ // mother exits
+ /* if( fStack->Particle(particle->GetMother(0))->GetPdgCode()==10441 ||//chic0
+ fStack->Particle(particle->GetMother(0))->GetPdgCode()==20443 ||//psi2S
+ fStack->Particle(particle->GetMother(0))->GetPdgCode()==445 //chic2
+ ){
+ fMCGammaChic.push_back(particle);
+ }
+ */
+ } // end if mother exits
+ } // end if particle is a photon
+
+
+
+ // process motherparticles (2 gammas as daughters)
+ // the motherparticle had already to pass the R and the eta cut, but no line cut.
+ // the line cut is just valid for the conversions!
+
+ if(particle->GetNDaughters() == 2){
+
+ TParticle* daughter0 = (TParticle*)fStack->Particle(particle->GetFirstDaughter());
+ TParticle* daughter1 = (TParticle*)fStack->Particle(particle->GetLastDaughter());
+
+ if(daughter0->GetPdgCode() != 22 || daughter1->GetPdgCode() != 22) continue; //check for gamma gamma daughters
+
+ // Check the acceptance for both gammas
+ Bool_t gammaEtaCut = kTRUE;
+ if(TMath::Abs(daughter0->Eta()) > fV0Reader->GetEtaCut() || TMath::Abs(daughter1->Eta()) > fV0Reader->GetEtaCut() ) gammaEtaCut = kFALSE;
+ Bool_t gammaRCut = kTRUE;
+ if(daughter0->R() > fV0Reader->GetMaxRCut() || daughter1->R() > fV0Reader->GetMaxRCut() ) gammaRCut = kFALSE;
+
+
+
+ // check for conversions now -> have to pass eta, R and line cut!
+ Bool_t daughter0Electron = kFALSE;
+ Bool_t daughter0Positron = kFALSE;
+ Bool_t daughter1Electron = kFALSE;
+ Bool_t daughter1Positron = kFALSE;
+
+ if(daughter0->GetNDaughters() >= 2){ // first gamma
+ for(Int_t TrackIndex=daughter0->GetFirstDaughter();TrackIndex<=daughter0->GetLastDaughter();TrackIndex++){
+ TParticle *tmpDaughter = fStack->Particle(TrackIndex);
+ if(tmpDaughter->GetUniqueID() == 5){
+ if(tmpDaughter->GetPdgCode() == 11){
+ if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
+ if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
+ if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
+ daughter0Electron = kTRUE;
+ }
+ }
+
+ }
+ }
+ else if(tmpDaughter->GetPdgCode() == -11){
+ if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
+ if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
+ if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
+ daughter0Positron = kTRUE;
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+ if(daughter1->GetNDaughters() >= 2){ // second gamma
+ for(Int_t TrackIndex=daughter1->GetFirstDaughter();TrackIndex<=daughter1->GetLastDaughter();TrackIndex++){
+ TParticle *tmpDaughter = fStack->Particle(TrackIndex);
+ if(tmpDaughter->GetUniqueID() == 5){
+ if(tmpDaughter->GetPdgCode() == 11){
+ if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
+ if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
+ if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
+ daughter1Electron = kTRUE;
+ }
+ }
+
+ }
+ }
+ else if(tmpDaughter->GetPdgCode() == -11){
+ if( TMath::Abs(tmpDaughter->Eta()) <= fV0Reader->GetEtaCut() ){
+ if( ( TMath::Abs(tmpDaughter->Vz()) * fV0Reader->GetLineCutZRSlope()) - fV0Reader->GetLineCutZValue() < tmpDaughter->R() ){
+ if(tmpDaughter->R()< fV0Reader->GetMaxRCut()){
+ daughter1Positron = kTRUE;
+ }
+ }
+
+ }
+
+ }
+ }
+ }
+ }
+
+
+
+
+ if(particle->GetPdgCode()==111){ //Pi0
+ if( iTracks >= fStack->GetNprimary()){
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Rapid", rapidity);
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_R", particle->R());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_ZR", particle->Vz(),particle->R());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_XY", particle->Vx(),particle->Vy());//only fill from one daughter to avoid multiple filling
+
+ if(gammaEtaCut && gammaRCut){
+ //if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
+ if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Secondaries_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
+ }
+ }
+ }
+ else{
+ fHistograms->FillHistogram("MC_Pi0_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Rapid", rapidity);
+ fHistograms->FillHistogram("MC_Pi0_Phi", tmpPhi);
+ fHistograms->FillHistogram("MC_Pi0_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_Pi0_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_Pi0_R", particle->R());
+ fHistograms->FillHistogram("MC_Pi0_ZR", particle->Vz(),particle->R());
+ fHistograms->FillHistogram("MC_Pi0_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
+ fHistograms->FillHistogram("MC_Pi0_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
+
+ if(gammaEtaCut && gammaRCut){
+ // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
+ fHistograms->FillHistogram("MC_Pi0_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
+ if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
+ fHistograms->FillHistogram("MC_Pi0_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Pi0_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
+ fHistograms->FillHistogram("MC_Pi0_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
+ }
+ }
+ }
+ }
+
+ if(particle->GetPdgCode()==221){ //Eta
+ fHistograms->FillHistogram("MC_Eta_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_Eta_Rapid", rapidity);
+ fHistograms->FillHistogram("MC_Eta_Phi",tmpPhi);
+ fHistograms->FillHistogram("MC_Eta_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_Eta_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_Eta_R", particle->R());
+ fHistograms->FillHistogram("MC_Eta_ZR", particle->Vz(),particle->R());
+ fHistograms->FillHistogram("MC_Eta_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
+ fHistograms->FillHistogram("MC_Eta_XY", particle->Vx(), particle->Vy());//only fill from one daughter to avoid multiple filling
+
+ if(gammaEtaCut && gammaRCut){
+ // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
+ fHistograms->FillHistogram("MC_Eta_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Eta_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
+ if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
+ fHistograms->FillHistogram("MC_Eta_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Eta_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
+ fHistograms->FillHistogram("MC_Eta_ZR_ConvGamma_withinAcceptance", particle->Vz(),particle->R());
+ }
+
+ }
+
+ }
+
+ // all motherparticles with 2 gammas as daughters
+ fHistograms->FillHistogram("MC_Mother_R", particle->R());
+ fHistograms->FillHistogram("MC_Mother_ZR", particle->Vz(),particle->R());
+ fHistograms->FillHistogram("MC_Mother_XY", particle->Vx(),particle->Vy());
+ fHistograms->FillHistogram("MC_Mother_Mass", particle->GetCalcMass());
+ fHistograms->FillHistogram("MC_Mother_GammaDaughter_OpeningAngle", GetMCOpeningAngle(daughter0,daughter1));
+ fHistograms->FillHistogram("MC_Mother_Energy", particle->Energy());
+ fHistograms->FillHistogram("MC_Mother_Pt", particle->Pt());
+ fHistograms->FillHistogram("MC_Mother_Eta", particle->Eta());
+ fHistograms->FillHistogram("MC_Mother_Rapid", rapidity);
+ fHistograms->FillHistogram("MC_Mother_Phi",tmpPhi);
+ fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt",particle->GetMass(),particle->Pt());
+
+ if(gammaEtaCut && gammaRCut){
+ // if(TMath::Abs(daughter0->Eta()) <= fV0Reader->GetEtaCut() && TMath::Abs(daughter1->Eta()) <= fV0Reader->GetEtaCut() ){
+ fHistograms->FillHistogram("MC_Mother_Pt_Eta_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Mother_Pt_Rapid_withinAcceptance", particle->Pt(),rapidity);
+ fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_withinAcceptance",particle->GetMass(),particle->Pt());
+ if(daughter0Electron && daughter0Positron && daughter1Electron && daughter1Positron){
+ fHistograms->FillHistogram("MC_Mother_Pt_Eta_ConvGamma_withinAcceptance", particle->Pt(),particle->Eta());
+ fHistograms->FillHistogram("MC_Mother_Pt_Rapid_ConvGamma_withinAcceptance", particle->Pt(),rapidity);
+ fHistograms->FillHistogram("MC_Mother_InvMass_vs_Pt_ConvGamma_withinAcceptance",particle->GetMass(),particle->Pt());
+
+ }
+
+
+ } // end passed R and eta cut
+
+ } // end if(particle->GetNDaughters() == 2)
+
+ }// end for (Int_t iTracks = 0; iTracks < fStack->GetNtrack(); iTracks++)
+
+} // end ProcessMCData
+
+
+
+void AliAnalysisTaskGammaConversion::FillNtuple(){
+ //Fills the ntuple with the different values
+
+ if(fGammaNtuple == NULL){
+ return;
+ }
+ Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
+ for(Int_t i=0;i<numberOfV0s;i++){
+ 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};
+ AliESDv0 * cV0 = fV0Reader->GetV0(i);
+ Double_t negPID=0;
+ Double_t posPID=0;
+ fV0Reader->GetPIDProbability(negPID,posPID);
+ values[0]=cV0->GetOnFlyStatus();
+ values[1]=fV0Reader->CheckForPrimaryVertex();
+ values[2]=negPID;
+ values[3]=posPID;
+ values[4]=fV0Reader->GetX();
+ values[5]=fV0Reader->GetY();
+ values[6]=fV0Reader->GetZ();
+ values[7]=fV0Reader->GetXYRadius();
+ values[8]=fV0Reader->GetMotherCandidateNDF();
+ values[9]=fV0Reader->GetMotherCandidateChi2();
+ values[10]=fV0Reader->GetMotherCandidateEnergy();
+ values[11]=fV0Reader->GetMotherCandidateEta();
+ values[12]=fV0Reader->GetMotherCandidatePt();
+ values[13]=fV0Reader->GetMotherCandidateMass();
+ values[14]=fV0Reader->GetMotherCandidateWidth();
+ // values[15]=fV0Reader->GetMotherMCParticle()->Pt(); MOVED TO THE END, HAS TO BE CALLED AFTER HasSameMother NB: still has the same entry in the array
+ values[16]=fV0Reader->GetOpeningAngle();
+ values[17]=fV0Reader->GetNegativeTrackEnergy();
+ values[18]=fV0Reader->GetNegativeTrackPt();
+ values[19]=fV0Reader->GetNegativeTrackEta();
+ values[20]=fV0Reader->GetNegativeTrackPhi();
+ values[21]=fV0Reader->GetPositiveTrackEnergy();
+ values[22]=fV0Reader->GetPositiveTrackPt();
+ values[23]=fV0Reader->GetPositiveTrackEta();
+ values[24]=fV0Reader->GetPositiveTrackPhi();
+ values[25]=fV0Reader->HasSameMCMother();
+ if(values[25] != 0){
+ values[26]=fV0Reader->GetMotherMCParticlePDGCode();
+ values[15]=fV0Reader->GetMotherMCParticle()->Pt();
+ }
+ fTotalNumberOfAddedNtupleEntries++;
+ fGammaNtuple->Fill(values);
+ }
+ fV0Reader->ResetV0IndexNumber();
+
+}
+
+void AliAnalysisTaskGammaConversion::ProcessV0sNoCut(){
+ // Process all the V0's without applying any cuts to it
+
+ Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
+ for(Int_t i=0;i<numberOfV0s;i++){
+ /*AliESDv0 * cV0 = */fV0Reader->GetV0(i);
+
+ if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
+ continue;
+ }
+
+ // if( !fV0Reader->GetV0(i)->GetOnFlyStatus()){
+ if( !fV0Reader->CheckV0FinderStatus(i)){
+ continue;
+ }
+
+
+ if( !((fV0Reader->GetNegativeESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ||
+ !((fV0Reader->GetPositiveESDTrack())->GetStatus() & AliESDtrack::kTPCrefit) ){
+ continue;
+ }
+
+ if(fDoMCTruth){
+
+ if(fV0Reader->HasSameMCMother() == kFALSE){
+ continue;
+ }
+
+ TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
+ TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
+
+ if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
+ continue;
+ }
+ if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
+ continue;
+ }
+
+ if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){
+ continue;
+ }
+
+ if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
+
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
+
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
+
+ fHistograms->FillHistogram("ESD_NoCutConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
+ fHistograms->FillHistogram("ESD_NoCutConversion_R", fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_NoCutConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_NoCutConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
+
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
+
+ //store MCTruth properties
+ fHistograms->FillHistogram("ESD_NoCutConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
+ fHistograms->FillHistogram("ESD_NoCutConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
+ fHistograms->FillHistogram("ESD_NoCutConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
+ }
+ }
+ }
+ fV0Reader->ResetV0IndexNumber();
+}
+
+void AliAnalysisTaskGammaConversion::ProcessV0s(){
+ // see header file for documentation
+
+ if(fWriteNtuple == kTRUE){
+ FillNtuple();
+ }
+
+ Int_t nSurvivingV0s=0;
+ while(fV0Reader->NextV0()){
+ nSurvivingV0s++;
+
+
+ //-------------------------- filling v0 information -------------------------------------
+ fHistograms->FillHistogram("ESD_Conversion_R", fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_Conversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_Conversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
+ fHistograms->FillHistogram("ESD_Conversion_OpeningAngle", fV0Reader->GetOpeningAngle());
+
+ fHistograms->FillHistogram("ESD_E_Energy", fV0Reader->GetNegativeTrackEnergy());
+ fHistograms->FillHistogram("ESD_E_Pt", fV0Reader->GetNegativeTrackPt());
+ fHistograms->FillHistogram("ESD_E_Eta", fV0Reader->GetNegativeTrackEta());
+ fHistograms->FillHistogram("ESD_E_Phi", fV0Reader->GetNegativeTrackPhi());
+
+ fHistograms->FillHistogram("ESD_P_Energy", fV0Reader->GetPositiveTrackEnergy());
+ fHistograms->FillHistogram("ESD_P_Pt", fV0Reader->GetPositiveTrackPt());
+ fHistograms->FillHistogram("ESD_P_Eta", fV0Reader->GetPositiveTrackEta());
+ fHistograms->FillHistogram("ESD_P_Phi", fV0Reader->GetPositiveTrackPhi());
+
+ fHistograms->FillHistogram("ESD_ConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
+ fHistograms->FillHistogram("ESD_ConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
+ fHistograms->FillHistogram("ESD_ConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
+ fHistograms->FillHistogram("ESD_ConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
+ fHistograms->FillHistogram("ESD_ConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
+ fHistograms->FillHistogram("ESD_ConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
+ fHistograms->FillHistogram("ESD_ConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_ConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
+ fHistograms->FillHistogram("ESD_ConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
+ fHistograms->FillHistogram("ESD_ConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
+
+ fHistograms->FillHistogram("ESD_ConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_ConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
+
+ fHistograms->FillHistogram("ESD_ConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
+ fHistograms->FillHistogram("ESD_ConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
+ fHistograms->FillHistogram("ESD_ConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
+ fHistograms->FillHistogram("ESD_ConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
+
+ fHistograms->FillHistogram("ESD_ConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_ConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_ConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
+ fHistograms->FillHistogram("ESD_ConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
+
+
+
+ // begin mapping
+ Int_t rBin = fHistograms->GetRBin(fV0Reader->GetXYRadius());
+ Int_t zBin = fHistograms->GetZBin(fV0Reader->GetZ());
+ Int_t phiBin = fHistograms->GetPhiBin(fV0Reader->GetNegativeTrackPhi());
+ TVector3 vtxConv(fV0Reader->GetX(),fV0Reader->GetY(), fV0Reader->GetZ());
+
+ Double_t motherCandidateEta= fV0Reader->GetMotherCandidateEta();
+
+ TString nameESDMappingPhiR="";
+ nameESDMappingPhiR.Form("ESD_Conversion_Mapping_Phi%02d_R%02d",phiBin,rBin);
+ fHistograms->FillHistogram(nameESDMappingPhiR, fV0Reader->GetZ(), motherCandidateEta);
+
+ TString nameESDMappingPhi="";
+ nameESDMappingPhi.Form("ESD_Conversion_Mapping_Phi%02d",phiBin);
+ fHistograms->FillHistogram(nameESDMappingPhi, fV0Reader->GetZ(), motherCandidateEta);
+
+ TString nameESDMappingR="";
+ nameESDMappingR.Form("ESD_Conversion_Mapping_R%02d",rBin);
+ fHistograms->FillHistogram(nameESDMappingR, fV0Reader->GetZ(), motherCandidateEta);
+
+ TString nameESDMappingPhiInR="";
+ nameESDMappingPhiInR.Form("ESD_Conversion_Mapping_Phi_in_R_%02d",rBin);
+ // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
+ fHistograms->FillHistogram(nameESDMappingPhiInR, vtxConv.Phi());
+
+ TString nameESDMappingZInR="";
+ nameESDMappingZInR.Form("ESD_Conversion_Mapping_Z_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameESDMappingZInR, fV0Reader->GetZ());
+
+ TString nameESDMappingPhiInZ="";
+ nameESDMappingPhiInZ.Form("ESD_Conversion_Mapping_Phi_in_Z_%02d",zBin);
+ // fHistograms->FillHistogram(nameESDMappingPhiInR, fV0Reader->GetMotherCandidatePhi());
+ fHistograms->FillHistogram(nameESDMappingPhiInZ, vtxConv.Phi());
+
+ TString nameESDMappingRInZ="";
+ nameESDMappingRInZ.Form("ESD_Conversion_Mapping_R_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameESDMappingRInZ, fV0Reader->GetXYRadius());
+
+ if(fV0Reader->GetMotherCandidatePt() > fLowPtMapping && fV0Reader->GetMotherCandidatePt()< fHighPtMapping){
+ TString nameESDMappingMidPtPhiInR="";
+ nameESDMappingMidPtPhiInR.Form("ESD_Conversion_Mapping_MidPt_Phi_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameESDMappingMidPtPhiInR, vtxConv.Phi());
+
+ TString nameESDMappingMidPtZInR="";
+ nameESDMappingMidPtZInR.Form("ESD_Conversion_Mapping_MidPt_Z_in_R_%02d",rBin);
+ fHistograms->FillHistogram(nameESDMappingMidPtZInR, fV0Reader->GetZ());
+
+ TString nameESDMappingMidPtPhiInZ="";
+ nameESDMappingMidPtPhiInZ.Form("ESD_Conversion_Mapping_MidPt_Phi_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameESDMappingMidPtPhiInZ, vtxConv.Phi());
+
+ TString nameESDMappingMidPtRInZ="";
+ nameESDMappingMidPtRInZ.Form("ESD_Conversion_Mapping_MidPt_R_in_Z_%02d",zBin);
+ fHistograms->FillHistogram(nameESDMappingMidPtRInZ, fV0Reader->GetXYRadius());
+
+ }
+
+
+ // end mapping
+
+ new((*fKFReconstructedGammasTClone)[fKFReconstructedGammasTClone->GetEntriesFast()]) AliKFParticle(*fV0Reader->GetMotherCandidateKFCombination());
+
+ // fKFReconstructedGammas.push_back(*fV0Reader->GetMotherCandidateKFCombination());
+ fElectronv1.push_back(fV0Reader->GetCurrentV0()->GetPindex());
+ fElectronv2.push_back(fV0Reader->GetCurrentV0()->GetNindex());
+
+ //----------------------------------- checking for "real" conversions (MC match) --------------------------------------
+ if(fDoMCTruth){
+
+ if(fV0Reader->HasSameMCMother() == kFALSE){
+ continue;
+ }
+ TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
+ TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
+
+ if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
+ continue;
+ }
+
+ if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
+ continue;
+ }
+
+ if(negativeMC->GetUniqueID() != 5 || positiveMC->GetUniqueID() !=5){// check if the daughters come from a conversion
+ continue;
+ }
+
+ if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
+
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Pt", fV0Reader->GetMotherCandidatePt());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Energy", fV0Reader->GetMotherCandidateEnergy());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Eta", fV0Reader->GetMotherCandidateEta());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Phi", fV0Reader->GetMotherCandidatePhi());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Mass", fV0Reader->GetMotherCandidateMass());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Width", fV0Reader->GetMotherCandidateWidth());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Chi2", fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_NDF", fV0Reader->GetMotherCandidateNDF());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Eta", fV0Reader->GetMotherCandidatePt(),fV0Reader->GetMotherCandidateEta());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Rapid", fV0Reader->GetMotherCandidateRapidity());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLength", /*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetNegativeTrackLength()*/fV0Reader->GetNegativeNTPCClusters(),fV0Reader->GetMotherCandidateMass());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_TrackLengthVSInvMass",/*fV0Reader->GetPositiveTrackLength()*/fV0Reader->GetPositiveNTPCClusters(),fV0Reader->GetMotherCandidateMass());
+
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Pt_Chi2", fV0Reader->GetMotherCandidatePt(), fV0Reader->GetMotherCandidateChi2());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_Eta_Chi2", fV0Reader->GetMotherCandidateEta(), fV0Reader->GetMotherCandidateChi2());
+
+
+ fHistograms->FillHistogram("ESD_TrueConversion_XY", fV0Reader->GetX(),fV0Reader->GetY());
+ fHistograms->FillHistogram("ESD_TrueConversion_R", fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_TrueConversion_ZR", fV0Reader->GetZ(),fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("ESD_TrueConversion_OpeningAngle", fV0Reader->GetOpeningAngle());
+
+ fHistograms->FillHistogram("ESD_TrueConvGamma_CosPointingAngle", fV0Reader->GetCosPointingAngle());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_DcaDaughters", fV0Reader->GetDcaDaughters());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_NormDcaDistDaughters", fV0Reader->GetNormDcaDistDaughters());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_LikelihoodAP", fV0Reader->GetLikelihoodAP());
+
+ fHistograms->FillHistogram("ESD_TrueConvGamma_E_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetNegativeTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_P_AsymmetryP",fV0Reader->GetMotherCandidateP(),fV0Reader->GetPositiveTrackP()/fV0Reader->GetMotherCandidateP());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_E_dEdxP",fV0Reader->GetNegativeTrackP(),fV0Reader->GetNegativeTrackTPCdEdx());
+ fHistograms->FillHistogram("ESD_TrueConvGamma_P_dEdxP",fV0Reader->GetPositiveTrackP(),fV0Reader->GetPositiveTrackTPCdEdx());
+
+
+
+ //store MCTruth properties
+ fHistograms->FillHistogram("ESD_TrueConvGamma_MC_Pt_Eta", fV0Reader->GetMotherMCParticle()->Pt(),fV0Reader->GetMotherMCParticle()->Eta());
+ fHistograms->FillHistogram("ESD_TrueConversion_MC_ZR", negativeMC->Vz(),negativeMC->R());
+ fHistograms->FillHistogram("ESD_TrueConversion_MC_XY", negativeMC->Vx(),negativeMC->Vy());
+
+ //resolution
+ Double_t mcpt = fV0Reader->GetMotherMCParticle()->Pt();
+ Double_t esdpt = fV0Reader->GetMotherCandidatePt();
+ Double_t resdPt = 0;
+ if(mcpt > 0){
+ resdPt = ((esdpt - mcpt)/mcpt)*100;
+ }
+ else if(mcpt < 0){
+ cout<<"Pt of MC particle is negative, this will cause wrong calculation of resPt"<<endl;
+ }
+
+ fHistograms->FillHistogram("Resolution_dPt", mcpt, resdPt);
+ fHistograms->FillHistogram("Resolution_MC_Pt", mcpt);
+ fHistograms->FillHistogram("Resolution_ESD_Pt", esdpt);
+
+ Double_t resdZ = 0;
+ if(fV0Reader->GetNegativeMCParticle()->Vz() != 0){
+ resdZ = ((fV0Reader->GetZ() -fV0Reader->GetNegativeMCParticle()->Vz())/fV0Reader->GetNegativeMCParticle()->Vz())*100;
+ }
+
+ fHistograms->FillHistogram("Resolution_dZ", fV0Reader->GetNegativeMCParticle()->Vz(), resdZ);
+ fHistograms->FillHistogram("Resolution_MC_Z", fV0Reader->GetNegativeMCParticle()->Vz());
+ fHistograms->FillHistogram("Resolution_ESD_Z", fV0Reader->GetZ());
+
+ Double_t resdR = 0;
+ if(fV0Reader->GetNegativeMCParticle()->R() != 0){
+ resdR = ((fV0Reader->GetXYRadius() - fV0Reader->GetNegativeMCParticle()->R())/fV0Reader->GetNegativeMCParticle()->R())*100;
+ }
+
+ fHistograms->FillHistogram("Resolution_dR", fV0Reader->GetNegativeMCParticle()->R(), resdR);
+ fHistograms->FillHistogram("Resolution_MC_R", fV0Reader->GetNegativeMCParticle()->R());
+ fHistograms->FillHistogram("Resolution_ESD_R", fV0Reader->GetXYRadius());
+ fHistograms->FillHistogram("Resolution_dR_dPt", resdR, resdPt);
+ }//if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22)
+ }//if(fDoMCTruth)
+ }//while(fV0Reader->NextV0)
+ fHistograms->FillHistogram("ESD_NumberOfSurvivingV0s", nSurvivingV0s);
+ fHistograms->FillHistogram("ESD_NumberOfV0s", fV0Reader->GetNumberOfV0s());
+ fHistograms->FillHistogram("ESD_NumberOfContributorsVtx", fV0Reader->GetNumberOfContributorsVtx());
+
+ fV0Reader->ResetV0IndexNumber();
+
+}
+
+void AliAnalysisTaskGammaConversion::FillAODWithConversionGammas(){
+ // Fill AOD with reconstructed Gamma
+
+ for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
+ // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
+ //Create AOD particle object from AliKFParticle
+
+ /* AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
+ //You could add directly AliKFParticle objects to the AOD, avoiding dependences with PartCorr
+ //but this means that I have to work a little bit more in my side.
+ //AODPWG4Particle objects are simpler and lighter, I think
+ AliAODPWG4Particle gamma = AliAODPWG4Particle(gammakf->Px(),gammakf->Py(),gammakf->Pz(), gammakf->E());
+ gamma.SetLabel(-1);//How to get the MC label of the reconstructed gamma?
+ gamma.SetCaloLabel(-1,-1); //How to get the MC label of the 2 electrons that form the gamma?
+ gamma.SetDetector("CTS"); //tag the gamma as reconstructed in the central barrel
+ gamma.SetPdg(AliCaloPID::kPhotonConv); //photon id
+ gamma.SetTag(-1); //Here I usually put a flag saying that montecarlo says it is prompt, decay fragmentation photon, or hadrons or whatever
+
+ //Add it to the aod list
+ Int_t i = fAODBranch->GetEntriesFast();
+ new((*fAODBranch)[i]) AliAODPWG4Particle(gamma);
+ */
+ // AliKFParticle * gammakf = &fKFReconstructedGammas[gammaIndex];
+ AliKFParticle * gammakf = (AliKFParticle *)fKFReconstructedGammasTClone->At(gammaIndex);
+ AliGammaConversionAODObject aodObject;
+ aodObject.SetPx(gammakf->GetPx());
+ aodObject.SetPy(gammakf->GetPy());
+ aodObject.SetPz(gammakf->GetPz());
+ aodObject.SetLabel1(fElectronv1[gammaIndex]);
+ aodObject.SetLabel2(fElectronv2[gammaIndex]);
+ Int_t i = fAODBranch->GetEntriesFast();
+ new((*fAODBranch)[i]) AliGammaConversionAODObject(aodObject);
+ }
+
+}
+
+
+void AliAnalysisTaskGammaConversion::ProcessGammasForNeutralMesonAnalysis(){
+ // see header file for documentation
+
+ // for(UInt_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammas.size();firstGammaIndex++){
+ // for(UInt_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammas.size();secondGammaIndex++){
+ for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
+ for(Int_t secondGammaIndex=firstGammaIndex+1;secondGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();secondGammaIndex++){
+
+ // AliKFParticle * twoGammaDecayCandidateDaughter0 = &fKFReconstructedGammas[firstGammaIndex];
+ // AliKFParticle * twoGammaDecayCandidateDaughter1 = &fKFReconstructedGammas[secondGammaIndex];
+
+ AliKFParticle * twoGammaDecayCandidateDaughter0 = (AliKFParticle *)fKFReconstructedGammasTClone->At(firstGammaIndex);
+ AliKFParticle * twoGammaDecayCandidateDaughter1 = (AliKFParticle *)fKFReconstructedGammasTClone->At(secondGammaIndex);
+
+ if(fElectronv1[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv1[firstGammaIndex]==fElectronv2[secondGammaIndex]){
+ continue;
+ }
+ if(fElectronv2[firstGammaIndex]==fElectronv1[secondGammaIndex] || fElectronv2[firstGammaIndex]==fElectronv2[secondGammaIndex]){
+ continue;
+ }
+
+ AliKFParticle *twoGammaCandidate = new AliKFParticle(*twoGammaDecayCandidateDaughter0,*twoGammaDecayCandidateDaughter1);
+
+ Double_t massTwoGammaCandidate = 0.;
+ Double_t widthTwoGammaCandidate = 0.;
+ Double_t chi2TwoGammaCandidate =10000.;
+ twoGammaCandidate->GetMass(massTwoGammaCandidate,widthTwoGammaCandidate);
+ if(twoGammaCandidate->GetNDF()>0){
+ chi2TwoGammaCandidate = twoGammaCandidate->GetChi2()/twoGammaCandidate->GetNDF();
+
+ if(chi2TwoGammaCandidate>0 && chi2TwoGammaCandidate<fV0Reader->GetChi2CutMeson()){
+
+ TVector3 momentumVectorTwoGammaCandidate(twoGammaCandidate->GetPx(),twoGammaCandidate->GetPy(),twoGammaCandidate->GetPz());
+ TVector3 spaceVectorTwoGammaCandidate(twoGammaCandidate->GetX(),twoGammaCandidate->GetY(),twoGammaCandidate->GetZ());
+
+ Double_t openingAngleTwoGammaCandidate = twoGammaDecayCandidateDaughter0->GetAngle(*twoGammaDecayCandidateDaughter1);
+ Double_t rapidity;
+ if(twoGammaCandidate->GetE() - twoGammaCandidate->GetPz() == 0 || twoGammaCandidate->GetE() + twoGammaCandidate->GetPz() == 0){
+ rapidity=0;
+ }
+ else{
+ rapidity = 0.5*(TMath::Log((twoGammaCandidate->GetE() +twoGammaCandidate->GetPz()) / (twoGammaCandidate->GetE()-twoGammaCandidate->GetPz())));
+ }
+
+ if(openingAngleTwoGammaCandidate < fMinOpeningAngleGhostCut){
+ delete twoGammaCandidate;
+ continue; // minimum opening angle to avoid using ghosttracks
+ }
+
+ fHistograms->FillHistogram("ESD_Mother_GammaDaughter_OpeningAngle", openingAngleTwoGammaCandidate);
+ fHistograms->FillHistogram("ESD_Mother_Energy", twoGammaCandidate->GetE());
+ fHistograms->FillHistogram("ESD_Mother_Pt", momentumVectorTwoGammaCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Mother_Eta", momentumVectorTwoGammaCandidate.Eta());
+ fHistograms->FillHistogram("ESD_Mother_Rapid", rapidity);
+ fHistograms->FillHistogram("ESD_Mother_Phi", spaceVectorTwoGammaCandidate.Phi());
+ fHistograms->FillHistogram("ESD_Mother_Mass", massTwoGammaCandidate);
+ fHistograms->FillHistogram("ESD_Mother_R", spaceVectorTwoGammaCandidate.Pt()); // Pt in Space == R!!!
+ fHistograms->FillHistogram("ESD_Mother_ZR", twoGammaCandidate->GetZ(), spaceVectorTwoGammaCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Mother_XY", twoGammaCandidate->GetX(), twoGammaCandidate->GetY());
+ fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Mother_InvMass",massTwoGammaCandidate);
+
+ if ( TMath::Abs(twoGammaDecayCandidateDaughter0->GetEta())<0.9 && TMath::Abs(twoGammaDecayCandidateDaughter1->GetEta())<0.9 ){
+ fHistograms->FillHistogram("ESD_Mother_InvMass_vs_Pt_Fiducial",massTwoGammaCandidate ,momentumVectorTwoGammaCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Mother_InvMass_Fiducial",massTwoGammaCandidate);
+ }
+ }
+ }
+ delete twoGammaCandidate;
+ }
+ }
+}
+
+void AliAnalysisTaskGammaConversion::CalculateBackground(){
+ // see header file for documentation
+
+
+ TClonesArray * currentEventV0s = fV0Reader->GetCurrentEventGoodV0s();
+
+ for(Int_t nEventsInBG=0;nEventsInBG <fV0Reader->GetNBGEvents();nEventsInBG++){
+ AliGammaConversionKFVector * previousEventV0s = fV0Reader->GetBGGoodV0s(nEventsInBG);
+ for(Int_t iCurrent=0;iCurrent<currentEventV0s->GetEntriesFast();iCurrent++){
+ AliKFParticle currentEventGoodV0 = *(AliKFParticle *)(currentEventV0s->At(iCurrent));
+ for(UInt_t iPrevious=0;iPrevious<previousEventV0s->size();iPrevious++){
+ AliKFParticle previousGoodV0 = (AliKFParticle)(*(previousEventV0s->at(iPrevious)));
+
+ AliKFParticle *backgroundCandidate = new AliKFParticle(currentEventGoodV0,previousGoodV0);
+
+ Double_t massBG =0.;
+ Double_t widthBG = 0.;
+ Double_t chi2BG =10000.;
+ backgroundCandidate->GetMass(massBG,widthBG);
+ if(backgroundCandidate->GetNDF()>0){
+ chi2BG = backgroundCandidate->GetChi2()/backgroundCandidate->GetNDF();
+ if(chi2BG>0 && chi2BG<fV0Reader->GetChi2CutMeson()){
+
+ TVector3 momentumVectorbackgroundCandidate(backgroundCandidate->GetPx(),backgroundCandidate->GetPy(),backgroundCandidate->GetPz());
+ TVector3 spaceVectorbackgroundCandidate(backgroundCandidate->GetX(),backgroundCandidate->GetY(),backgroundCandidate->GetZ());
+
+ Double_t openingAngleBG = currentEventGoodV0.GetAngle(previousGoodV0);
+
+ Double_t rapidity;
+ if(backgroundCandidate->GetE() - backgroundCandidate->GetPz() == 0 || backgroundCandidate->GetE() + backgroundCandidate->GetPz() == 0) rapidity=0;
+ else rapidity = 0.5*(TMath::Log((backgroundCandidate->GetE() +backgroundCandidate->GetPz()) / (backgroundCandidate->GetE()-backgroundCandidate->GetPz())));
+
+
+
+
+ if(openingAngleBG < fMinOpeningAngleGhostCut ){
+ delete backgroundCandidate;
+ continue; // minimum opening angle to avoid using ghosttracks
+ }
+
+ fHistograms->FillHistogram("ESD_Background_GammaDaughter_OpeningAngle", openingAngleBG);
+ fHistograms->FillHistogram("ESD_Background_Energy", backgroundCandidate->GetE());
+ fHistograms->FillHistogram("ESD_Background_Pt", momentumVectorbackgroundCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Background_Eta", momentumVectorbackgroundCandidate.Eta());
+ fHistograms->FillHistogram("ESD_Background_Rapidity", rapidity);
+ fHistograms->FillHistogram("ESD_Background_Phi", spaceVectorbackgroundCandidate.Phi());
+ fHistograms->FillHistogram("ESD_Background_Mass", massBG);
+ fHistograms->FillHistogram("ESD_Background_R", spaceVectorbackgroundCandidate.Pt()); // Pt in Space == R!!!!
+ fHistograms->FillHistogram("ESD_Background_ZR", backgroundCandidate->GetZ(), spaceVectorbackgroundCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Background_XY", backgroundCandidate->GetX(), backgroundCandidate->GetY());
+ fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt",massBG,momentumVectorbackgroundCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Background_InvMass",massBG);
+
+ if ( TMath::Abs(currentEventGoodV0.GetEta())<0.9 && TMath::Abs(previousGoodV0.GetEta())<0.9 ){
+ fHistograms->FillHistogram("ESD_Background_InvMass_vs_Pt_Fiducial",massBG,momentumVectorbackgroundCandidate.Pt());
+ fHistograms->FillHistogram("ESD_Background_InvMass_Fiducial",massBG);
+ }
+
+ }
+ }
+ delete backgroundCandidate;
+ }
+ }
+ }
+}
+
+
+void AliAnalysisTaskGammaConversion::ProcessGammasForGammaJetAnalysis(){
+ //ProcessGammasForGammaJetAnalysis
+
+ Double_t distIsoMin;
+
+ CreateListOfChargedParticles();
+
+
+ // for(UInt_t gammaIndex=0;gammaIndex<fKFReconstructedGammas.size();gammaIndex++){
+ for(Int_t gammaIndex=0;gammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();gammaIndex++){
+ AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
+ TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
+
+ if( momentumVectorCurrentGamma.Pt()> fMinPtForGammaJet){
+ distIsoMin=GetMinimumDistanceToCharge(gammaIndex);
+
+ if (distIsoMin > fMinIsoConeSize && fLeadingChargedIndex>=0){
+ CalculateJetCone(gammaIndex);
+ }
+ }
+ }
+}
+
+void AliAnalysisTaskGammaConversion::CreateListOfChargedParticles(){
+ // CreateListOfChargedParticles
+
+ fESDEvent = fV0Reader->GetESDEvent();
+ for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
+ AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
+
+ if(!curTrack){
+ continue;
+ }
+
+ if(fEsdTrackCuts->AcceptTrack(curTrack) ){
+ new((*fChargedParticles)[fChargedParticles->GetEntriesFast()]) AliESDtrack(*curTrack);
+ // fChargedParticles.push_back(curTrack);
+ fChargedParticlesId.push_back(iTracks);
+ }
+ }
+}
+void AliAnalysisTaskGammaConversion::CalculateJetCone(Int_t gammaIndex){
+ // CaculateJetCone
+
+ Double_t cone;
+ Double_t coneSize=0.3;
+ Double_t ptJet=0;
+
+ // AliKFParticle * currentGamma = &fKFReconstructedGammas[gammaIndex];
+ AliKFParticle * currentGamma = (AliKFParticle*)fKFReconstructedGammasTClone->At(gammaIndex);
+
+ TVector3 momentumVectorCurrentGamma(currentGamma->GetPx(),currentGamma->GetPy(),currentGamma->GetPz());
+
+ AliESDtrack* leadingCharged = (AliESDtrack*)(fChargedParticles->At(fLeadingChargedIndex));
+
+ Double_t momLeadingCharged[3];
+ leadingCharged->GetConstrainedPxPyPz(momLeadingCharged);
+
+ TVector3 momentumVectorLeadingCharged(momLeadingCharged[0],momLeadingCharged[1],momLeadingCharged[2]);
+
+ Double_t phi1=momentumVectorLeadingCharged.Phi();
+ Double_t eta1=momentumVectorLeadingCharged.Eta();
+ Double_t phi3=momentumVectorCurrentGamma.Phi();
+
+ for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
+ AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
+ Int_t chId = fChargedParticlesId[iCh];
+ if(fLeadingChargedIndex==chId || fLeadingChargedIndex==chId) continue;
+ Double_t mom[3];
+ curTrack->GetConstrainedPxPyPz(mom);
+ TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
+ Double_t phi2=momentumVectorChargedParticle.Phi();
+ Double_t eta2=momentumVectorChargedParticle.Eta();
+
+
+ cone=100.;
+ if( TMath::Abs(phi2 - phi1) <= ( TMath::TwoPi()-coneSize) ){
+ cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-phi1),2) );
+ }else{
+ if( (phi2 - phi1)> TMath::TwoPi()-coneSize ){
+ cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2-TMath::TwoPi()-phi1),2) );
+ }
+ if( (phi2 - phi1)< -(TMath::TwoPi()-coneSize) ){
+ cone = TMath::Sqrt( TMath::Power((eta2-eta1),2)+ TMath::Power((phi2+TMath::TwoPi()-phi1),2) );
+ }
+ }
+
+ if(cone <coneSize&& momentumVectorChargedParticle.Pt()>fMinPtJetCone ){
+ ptJet+= momentumVectorChargedParticle.Pt();
+ Double_t ffzHdrGam = momentumVectorChargedParticle.Pt()/momentumVectorCurrentGamma.Pt();
+ Double_t imbalanceHdrGam=-momentumVectorChargedParticle.Dot(momentumVectorCurrentGamma)/momentumVectorCurrentGamma.Mag2();
+ fHistograms->FillHistogram("ESD_FFzHdrGam",ffzHdrGam);
+ fHistograms->FillHistogram("ESD_ImbalanceHdrGam",imbalanceHdrGam);
+
+ }
+
+ Double_t dphiHdrGam=phi3-phi2;
+ if ( dphiHdrGam < (-TMath::PiOver2())){
+ dphiHdrGam+=(TMath::TwoPi());
+ }
+
+ if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
+ dphiHdrGam-=(TMath::TwoPi());
+ }
+
+ if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
+ fHistograms->FillHistogram("ESD_dphiHdrGamIsolated",dphiHdrGam);
+ }
+ }//track loop
+
+
+}
+
+Double_t AliAnalysisTaskGammaConversion::GetMinimumDistanceToCharge(Int_t indexHighestPtGamma){
+ // GetMinimumDistanceToCharge
+
+ Double_t fIsoMin=100.;
+ Double_t ptLeadingCharged=-1.;
+
+ fLeadingChargedIndex=-1;
+
+ AliKFParticle * gammaHighestPt = (AliKFParticle*)fKFReconstructedGammasTClone->At(indexHighestPtGamma);
+ TVector3 momentumVectorgammaHighestPt(gammaHighestPt->GetPx(),gammaHighestPt->GetPy(),gammaHighestPt->GetPz());
+
+ Double_t phi1=momentumVectorgammaHighestPt.Phi();
+ Double_t eta1=momentumVectorgammaHighestPt.Eta();
+
+ for(Int_t iCh=0;iCh<fChargedParticles->GetEntriesFast();iCh++){
+ AliESDtrack* curTrack = (AliESDtrack*)(fChargedParticles->At(iCh));
+ Int_t chId = fChargedParticlesId[iCh];
+ if(fElectronv1[indexHighestPtGamma]==chId || fElectronv2[indexHighestPtGamma]==chId) continue;
+ Double_t mom[3];
+ curTrack->GetConstrainedPxPyPz(mom);
+ TVector3 momentumVectorChargedParticle(mom[0],mom[1],mom[2]);
+ Double_t phi2=momentumVectorChargedParticle.Phi();
+ Double_t eta2=momentumVectorChargedParticle.Eta();
+ Double_t iso=pow( (pow( (eta1-eta2),2)+ pow((phi1-phi2),2)),0.5 );
+
+ if(momentumVectorChargedParticle.Pt()>fMinPtIsoCone ){
+ if (iso<fIsoMin){
+ fIsoMin=iso;
+ }
+ }
+
+ Double_t dphiHdrGam=phi1-phi2;
+ if ( dphiHdrGam < (-TMath::PiOver2())){
+ dphiHdrGam+=(TMath::TwoPi());
+ }
+
+ if ( dphiHdrGam > (3.*TMath::PiOver2()) ){
+ dphiHdrGam-=(TMath::TwoPi());
+ }
+ if (momentumVectorChargedParticle.Pt()>fMinPtGamChargedCorr){
+ fHistograms->FillHistogram("ESD_dphiHdrGam",dphiHdrGam);
+ }
+
+ if (dphiHdrGam>0.9*TMath::Pi() && dphiHdrGam<1.1*TMath::Pi()){
+ if (momentumVectorChargedParticle.Pt()> ptLeadingCharged && momentumVectorChargedParticle.Pt()>0.1*momentumVectorgammaHighestPt.Pt()){
+ ptLeadingCharged=momentumVectorChargedParticle.Pt();
+ fLeadingChargedIndex=iCh;
+ }
+ }
+
+ }//track loop
+ fHistograms->FillHistogram("ESD_MinimumIsoDistance",fIsoMin);
+ return fIsoMin;
+
+}
+
+Int_t AliAnalysisTaskGammaConversion::GetIndexHighestPtGamma(){
+ //GetIndexHighestPtGamma
+
+ Int_t indexHighestPtGamma=-1;
+ //Double_t
+ fGammaPtHighest = -100.;
+
+ for(Int_t firstGammaIndex=0;firstGammaIndex<fKFReconstructedGammasTClone->GetEntriesFast();firstGammaIndex++){
+ AliKFParticle * gammaHighestPtCandidate = (AliKFParticle*)fKFReconstructedGammasTClone->At(firstGammaIndex);
+ TVector3 momentumVectorgammaHighestPtCandidate(gammaHighestPtCandidate->GetPx(),gammaHighestPtCandidate->GetPy(),gammaHighestPtCandidate->GetPz());
+ if (momentumVectorgammaHighestPtCandidate.Pt() > fGammaPtHighest){
+ fGammaPtHighest=momentumVectorgammaHighestPtCandidate.Pt();
+ //gammaHighestPt = gammaHighestPtCandidate;
+ indexHighestPtGamma=firstGammaIndex;
+ }
+ }
+
+ return indexHighestPtGamma;
+
+}
+
+
+void AliAnalysisTaskGammaConversion::Terminate(Option_t */*option*/)
+{
+ // Terminate analysis
+ //
+ AliDebug(1,"Do nothing in Terminate");
+}
+
+void AliAnalysisTaskGammaConversion::UserCreateOutputObjects()
+{
+ //AOD
+ if(fAODBranch==NULL){
+ fAODBranch = new TClonesArray("AliGammaConversionAODObject", 0);
+ }
+ fAODBranch->Delete();
+ fAODBranch->SetName(fAODBranchName);
+ AddAODBranch("TClonesArray", &fAODBranch);
+
+ // Create the output container
+ if(fOutputContainer != NULL){
+ delete fOutputContainer;
+ fOutputContainer = NULL;
+ }
+ if(fOutputContainer == NULL){
+ fOutputContainer = new TList();
+ fOutputContainer->SetOwner(kTRUE);
+ }
+
+ //Adding the histograms to the output container
+ fHistograms->GetOutputContainer(fOutputContainer);
+
+
+ if(fWriteNtuple){
+ if(fGammaNtuple == NULL){
+ 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);
+ }
+ if(fNeutralMesonNtuple == NULL){
+ fNeutralMesonNtuple = new TNtuple("NeutralMesonNtuple","NeutralMesonNtuple","test");
+ }
+ TList * ntupleTList = new TList();
+ ntupleTList->SetOwner(kTRUE);
+ ntupleTList->SetName("Ntuple");
+ ntupleTList->Add((TNtuple*)fGammaNtuple);
+ fOutputContainer->Add(ntupleTList);
+ }
+
+ fOutputContainer->SetName(GetName());
+}
+
+Double_t AliAnalysisTaskGammaConversion::GetMCOpeningAngle(TParticle* const daughter0, TParticle* const daughter1) const{
+ //helper function
+ TVector3 v3D0(daughter0->Px(),daughter0->Py(),daughter0->Pz());
+ TVector3 v3D1(daughter1->Px(),daughter1->Py(),daughter1->Pz());
+ return v3D0.Angle(v3D1);
+}
+
+void AliAnalysisTaskGammaConversion::CheckV0Efficiency(){
+ // see header file for documentation
+
+ vector<Int_t> indexOfGammaParticle;
+
+ fStack = fV0Reader->GetMCStack();
+
+ if(fV0Reader->CheckForPrimaryVertex() == kFALSE){
+ return; // aborts if the primary vertex does not have contributors.
+ }
+
+ for (Int_t iTracks = 0; iTracks < fStack->GetNprimary(); iTracks++) {
+ TParticle* particle = (TParticle *)fStack->Particle(iTracks);
+ if(particle->GetPdgCode()==22){ //Gamma
+ if(particle->GetNDaughters() >= 2){
+ TParticle* electron=NULL;
+ TParticle* positron=NULL;
+ for(Int_t daughterIndex=particle->GetFirstDaughter();daughterIndex<=particle->GetLastDaughter();daughterIndex++){
+ TParticle *tmpDaughter = fStack->Particle(daughterIndex);
+ if(tmpDaughter->GetUniqueID() == 5){
+ if(tmpDaughter->GetPdgCode() == 11){
+ electron = tmpDaughter;
+ }
+ else if(tmpDaughter->GetPdgCode() == -11){
+ positron = tmpDaughter;
+ }
+ }
+ }
+ if(electron!=NULL && positron!=0){
+ if(electron->R()<160){
+ indexOfGammaParticle.push_back(iTracks);
+ }
+ }
+ }
+ }
+ }
+
+ Int_t nFoundGammas=0;
+ Int_t nNotFoundGammas=0;
+
+ Int_t numberOfV0s = fV0Reader->GetNumberOfV0s();
+ for(Int_t i=0;i<numberOfV0s;i++){
+ fV0Reader->GetV0(i);
+
+ if(fV0Reader->HasSameMCMother() == kFALSE){
+ continue;
+ }
+
+ TParticle * negativeMC = (TParticle*)fV0Reader->GetNegativeMCParticle();
+ TParticle * positiveMC = (TParticle*)fV0Reader->GetPositiveMCParticle();
+
+ if(TMath::Abs(negativeMC->GetPdgCode())!=11 || TMath::Abs(positiveMC->GetPdgCode())!=11){
+ continue;
+ }
+ if(negativeMC->GetPdgCode()==positiveMC->GetPdgCode()){
+ continue;
+ }
+
+ if(fV0Reader->GetMotherMCParticle()->GetPdgCode() == 22){
+ //TParticle * v0Gamma = fV0Reader->GetMotherMCParticle();
+ for(UInt_t mcIndex=0;mcIndex<indexOfGammaParticle.size();mcIndex++){
+ if(negativeMC->GetFirstMother()==indexOfGammaParticle[mcIndex]){
+ nFoundGammas++;
+ }
+ else{
+ nNotFoundGammas++;
+ }
+ }
+ }
+ }
+}
+
+
+void AliAnalysisTaskGammaConversion::ProcessGammaElectronsForChicAnalysis(){
+ // see header file for documantation
+
+ fESDEvent = fV0Reader->GetESDEvent();
+
+
+ TClonesArray * vESDeNegTemp = new TClonesArray("AliESDtrack",0);
+ TClonesArray * vESDePosTemp = new TClonesArray("AliESDtrack",0);
+ TClonesArray * vESDxNegTemp = new TClonesArray("AliESDtrack",0);
+ TClonesArray * vESDxPosTemp = new TClonesArray("AliESDtrack",0);
+ TClonesArray * vESDeNegNoJPsi = new TClonesArray("AliESDtrack",0);
+ TClonesArray * vESDePosNoJPsi = new TClonesArray("AliESDtrack",0);
+
+ /*
+ vector <AliESDtrack*> vESDeNegTemp(0);
+ vector <AliESDtrack*> vESDePosTemp(0);
+ vector <AliESDtrack*> vESDxNegTemp(0);
+ vector <AliESDtrack*> vESDxPosTemp(0);
+ vector <AliESDtrack*> vESDeNegNoJPsi(0);
+ vector <AliESDtrack*> vESDePosNoJPsi(0);
+ */
+
+
+ fHistograms->FillTable("Table_Electrons",0);//Count number of Events
+
+ for(Int_t iTracks = 0; iTracks < fESDEvent->GetNumberOfTracks(); iTracks++){
+ AliESDtrack* curTrack = fESDEvent->GetTrack(iTracks);
+
+ if(!curTrack){
+ //print warning here
+ continue;
+ }
+
+ double p[3];if(!curTrack->GetConstrainedPxPyPz(p))continue;
+ double r[3];curTrack->GetConstrainedXYZ(r);
+
+ TVector3 rXYZ(r);
+
+ fHistograms->FillTable("Table_Electrons",4);//Count number of ESD tracks
+
+ Bool_t flagKink = kTRUE;
+ Bool_t flagTPCrefit = kTRUE;
+ Bool_t flagTRDrefit = kTRUE;
+ Bool_t flagITSrefit = kTRUE;
+ Bool_t flagTRDout = kTRUE;
+ Bool_t flagVertex = kTRUE;
+
+
+ //Cuts ---------------------------------------------------------------
+
+ if(curTrack->GetKinkIndex(0) > 0){
+ fHistograms->FillHistogram("Table_Electrons",5);//Count kink
+ flagKink = kFALSE;
+ }
+
+ ULong_t trkStatus = curTrack->GetStatus();
+
+ ULong_t tpcRefit = (trkStatus & AliESDtrack::kTPCrefit);
+
+ if(!tpcRefit){
+ fHistograms->FillHistogram("Table_Electrons",9);//Count not TPCrefit
+ flagTPCrefit = kFALSE;
+ }
+
+ ULong_t itsRefit = (trkStatus & AliESDtrack::kITSrefit);
+ if(!itsRefit){
+ fHistograms->FillHistogram("Table_Electrons",10);//Count not ITSrefit
+ flagITSrefit = kFALSE;
+ }
+
+ ULong_t trdRefit = (trkStatus & AliESDtrack::kTRDrefit);
+
+ if(!trdRefit){
+ fHistograms->FillHistogram("Table_Electrons",8); //Count not TRDrefit
+ flagTRDrefit = kFALSE;
+ }
+
+ ULong_t trdOut = (trkStatus & AliESDtrack::kTRDout);
+
+ if(!trdOut) {
+ fHistograms->FillHistogram("Table_Electrons",7); //Count not TRDout
+ flagTRDout = kFALSE;
+ }
+
+ double nsigmaToVxt = GetSigmaToVertex(curTrack);
+
+ if(nsigmaToVxt > 3){
+ fHistograms->FillHistogram("Table_Electrons",6); //Count Tracks with number of sigmas > 3
+ flagVertex = kFALSE;
+ }
+
+ if(! (flagKink && flagTPCrefit && flagITSrefit && flagTRDrefit && flagTRDout && flagVertex ) ) continue;
+ fHistograms->FillHistogram("Table_Electrons",11);//Count Tracks passed Cuts
+
+
+ Stat_t pid, weight;
+ GetPID(curTrack, pid, weight);
+
+ if(pid!=0){
+ fHistograms->FillHistogram("Table_Electrons",12); //Count Tracks with pid != 0
+ }
+
+ if(pid == 0){
+ fHistograms->FillHistogram("Table_Electrons",13); //Count Tracks with pid != 0
+ }
+
+
+
+
+
+
+ TLorentzVector curElec;
+ curElec.SetXYZM(p[0],p[1],p[2],fElectronMass);
+
+
+ if(fDoMCTruth){
+ Int_t labelMC = TMath::Abs(curTrack->GetLabel());
+ TParticle* curParticle = fStack->Particle(labelMC);
+ if(curTrack->GetSign() > 0){
+ if( pid == 0){
+ fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
+ fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
+ }
+ else{
+ fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
+ fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
+ }
+ }
+ }
+
+
+ if(curTrack->GetSign() > 0){
+
+ // vESDxPosTemp.push_back(curTrack);
+ new((*vESDxPosTemp)[vESDxPosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
+
+ if( pid == 0){
+
+ fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
+ fHistograms->FillHistogram("ESD_ElectronPosPt",curElec.Pt());
+ // fHistograms->FillHistogram("MC_ElectronPosNegPt",curParticle->Pt());
+ fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
+ // fHistograms->FillHistogram("MC_ElectronPosNegEta",curParticle->Eta());
+ // vESDePosTemp.push_back(curTrack);
+ new((*vESDePosTemp)[vESDePosTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
+
+ }
+
+ }
+ else {
+
+ new((*vESDxNegTemp)[vESDxNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
+
+ if( pid == 0){
+
+ fHistograms->FillHistogram("ESD_ElectronPosNegPt",curElec.Pt());
+ fHistograms->FillHistogram("ESD_ElectronNegPt",curElec.Pt());
+ fHistograms->FillHistogram("ESD_ElectronPosNegEta",curElec.Eta());
+ new((*vESDeNegTemp)[vESDeNegTemp->GetEntriesFast()]) AliESDtrack(*curTrack);
+
+ }
+
+ }
+
+ }
+
+
+ Bool_t ePosJPsi = kFALSE;
+ Bool_t eNegJPsi = kFALSE;
+ Bool_t ePosPi0 = kFALSE;
+ Bool_t eNegPi0 = kFALSE;
+
+ UInt_t iePosJPsi=0,ieNegJPsi=0,iePosPi0=0,ieNegPi0=0;
+
+ for(Int_t iNeg=0; iNeg < vESDeNegTemp->GetEntriesFast(); iNeg++){
+ if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetPdgCode() == 11)
+ if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0) > -1){
+ Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(iNeg)))->GetLabel()))->GetMother(0);
+ TParticle* partMother = fStack ->Particle(labelMother);
+ if (partMother->GetPdgCode() == 111){
+ ieNegPi0 = iNeg;
+ eNegPi0 = kTRUE;
+ }
+ if(partMother->GetPdgCode() == 443){ //Mother JPsi
+ fHistograms->FillTable("Table_Electrons",14);
+ ieNegJPsi = iNeg;
+ eNegJPsi = kTRUE;
+ }
+ else{
+ // vESDeNegNoJPsi.push_back(vESDeNegTemp[iNeg]);
+ new((*vESDeNegNoJPsi)[vESDeNegNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDeNegTemp->At(iNeg)));
+ // cout<<"ESD No Positivo JPsi "<<endl;
+ }
+
+ }
+ }
+
+ for(Int_t iPos=0; iPos < vESDePosTemp->GetEntriesFast(); iPos++){
+ if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetPdgCode() == -11)
+ if(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0) > -1){
+ Int_t labelMother = fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iPos)))->GetLabel()))->GetMother(0);
+ TParticle* partMother = fStack ->Particle(labelMother);
+ if (partMother->GetPdgCode() == 111){
+ iePosPi0 = iPos;
+ ePosPi0 = kTRUE;
+ }
+ if(partMother->GetPdgCode() == 443){ //Mother JPsi
+ fHistograms->FillTable("Table_Electrons",15);
+ iePosJPsi = iPos;
+ ePosJPsi = kTRUE;
+ }
+ else{
+ // vESDePosNoJPsi.push_back(vESDePosTemp[iPos]);
+ new((*vESDePosNoJPsi)[vESDePosNoJPsi->GetEntriesFast()]) AliESDtrack(*(AliESDtrack*)(vESDePosTemp->At(iPos)));
+ // cout<<"ESD No Negativo JPsi "<<endl;
+ }
+
+ }
+ }
+
+ if( eNegJPsi && ePosJPsi ){
+ TVector3 tempeNegV,tempePosV;
+ tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->Pz());
+ tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->Pz());
+ fHistograms->FillTable("Table_Electrons",16);
+ fHistograms->FillHistogram("ESD_ElectronPosNegJPsiAngle",tempeNegV.Angle(tempePosV));
+ fHistograms->FillHistogram("MC_ElectronPosNegJPsiAngle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegJPsi)))->GetLabel())),
+ fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosJPsi)))->GetLabel()))));
+ }
+
+ if( eNegPi0 && ePosPi0 ){
+ TVector3 tempeNegV,tempePosV;
+ tempeNegV.SetXYZ(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Px(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Py(),((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->Pz());
+ tempePosV.SetXYZ(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Px(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Py(),((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->Pz());
+ fHistograms->FillHistogram("ESD_ElectronPosNegPi0Angle",tempeNegV.Angle(tempePosV));
+ fHistograms->FillHistogram("MC_ElectronPosNegPi0Angle",GetMCOpeningAngle(fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDeNegTemp->At(ieNegPi0)))->GetLabel())),
+ fStack->Particle(TMath::Abs(((AliESDtrack*)(vESDePosTemp->At(iePosPi0)))->GetLabel()))));
+ }
+
+
+ FillAngle("ESD_eNegePosAngleBeforeCut",GetTLorentzVector(vESDeNegTemp),GetTLorentzVector(vESDePosTemp));
+
+ CleanWithAngleCuts(*vESDeNegTemp,*vESDePosTemp,*fKFReconstructedGammasTClone);
+
+ // vector <TLorentzVector> vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectron);
+ // vector <TLorentzVector> vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectron);
+
+ TClonesArray vCurrentTLVeNeg = GetTLorentzVector(fCurrentEventNegElectronTClone);
+ TClonesArray vCurrentTLVePos = GetTLorentzVector(fCurrentEventPosElectronTClone);
+
+
+ FillAngle("ESD_eNegePosAngleAfterCut",vCurrentTLVeNeg,vCurrentTLVePos);
+
+
+
+
+ //FillAngle("ESD_eNegePosAngleAfterCut",CurrentTLVeNeg,CurrentTLVePos);
+
+
+ FillElectronInvMass("ESD_InvMass_ePluseMinus",vCurrentTLVeNeg,vCurrentTLVePos);
+ FillElectronInvMass("ESD_InvMass_xPlusxMinus",GetTLorentzVector(vESDxNegTemp),GetTLorentzVector(vESDxPosTemp));
+
+
+
+ FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusChiC","ESD_InvMass_GammaePluseMinusChiCDiff",*fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
+
+ FillGammaElectronInvMass("ESD_InvMass_GammaePluseMinusPi0","ESD_InvMass_GammaePluseMinusPi0Diff",
+ *fKFReconstructedGammasCutTClone,vCurrentTLVeNeg,vCurrentTLVePos);
+
+ //BackGround
+
+ //Like Sign e+e-
+ ElectronBackground("ESD_ENegBackground",vCurrentTLVeNeg);
+ ElectronBackground("ESD_EPosBackground",vCurrentTLVePos);
+ ElectronBackground("ESD_EPosENegBackground",vCurrentTLVeNeg);
+ ElectronBackground("ESD_EPosENegBackground",vCurrentTLVePos);
+
+ // Like Sign e+e- no JPsi
+ ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDeNegNoJPsi));
+ ElectronBackground("ESD_EPosENegNoJPsiBG",GetTLorentzVector(vESDePosNoJPsi));
+
+ //Mixed Event
+
+ if( fCurrentEventPosElectronTClone->GetEntriesFast() > 0 && fCurrentEventNegElectronTClone->GetEntriesFast() > 0 && fKFReconstructedGammasCutTClone->GetEntriesFast() > 0 ){
+ FillGammaElectronInvMass("ESD_EPosENegGammaBackgroundMX","ESD_EPosENegGammaBackgroundMXDiff",
+ *fKFReconstructedGammasCutTClone,*fPreviousEventTLVNegElectronTClone,*fPreviousEventTLVPosElectronTClone);
+ *fPreviousEventTLVNegElectronTClone = vCurrentTLVeNeg;
+ *fPreviousEventTLVPosElectronTClone = vCurrentTLVePos;
+
+ }
+
+ /*
+ //Photons P
+ Double_t vtx[3];
+ vtx[0]=0;vtx[1]=0;vtx[2]=0;
+ for(UInt_t i=0;i<fKFReconstructedGammasChic.size();i++){
+
+ // if(fMCGammaChicTempCut[i]->GetMother(0) < 0) continue;
+
+
+
+ Int_t tempLabel = fStack->Particle(fMCGammaChicTempCut[i]->GetMother(0))->GetPdgCode();
+ // cout<<" Label Pedro Gonzalez " <<tempLabel <<endl;
+
+ // cout<<" Label Distance"<<fKFReconstructedGammasChic[i].GetDistanceFromVertex(vtx)<<endl;
+
+ if( tempLabel == 10441 || tempLabel == 20443 || tempLabel == 445 )
+
+ fHistograms->FillHistogram("ESD_PhotonsMomentum",fKFReconstructedGammasChic[i].GetMomentum());
+
+
+ }
+
+
+ */
+
+
+ vESDeNegTemp->Delete();
+ vESDePosTemp->Delete();
+ vESDxNegTemp->Delete();
+ vESDxPosTemp->Delete();
+ vESDeNegNoJPsi->Delete();
+ vESDePosNoJPsi->Delete();
+
+ delete vESDeNegTemp;
+ delete vESDePosTemp;
+ delete vESDxNegTemp;
+ delete vESDxPosTemp;
+ delete vESDeNegNoJPsi;
+ delete vESDePosNoJPsi;
+}
+
+/*
+ void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,vector <TLorentzVector> tlVeNeg, vector <TLorentzVector> tlVePos){
+ //see header file for documentation
+ for( UInt_t iNeg=0; iNeg < tlVeNeg.size(); iNeg++){
+ for (UInt_t iPos=0; iPos < tlVePos.size(); iPos++){
+ fHistograms->FillHistogram(histoName.Data(),tlVeNeg[iNeg].Vect().Angle(tlVePos[iPos].Vect()));
+ }
+ }
+ }
+*/
+void AliAnalysisTaskGammaConversion::FillAngle(TString histoName,TClonesArray const tlVeNeg, TClonesArray const tlVePos){
+ //see header file for documentation
+ for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++){
+ for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
+ fHistograms->FillHistogram(histoName.Data(),((TLorentzVector*)(tlVeNeg.At(iNeg)))->Vect().Angle(((TLorentzVector*)(tlVePos.At(iPos)))->Vect()));
+ }
+ }
+}
+void AliAnalysisTaskGammaConversion::FillElectronInvMass(TString histoName, TClonesArray const eNeg, TClonesArray const ePos){
+ //see header file for documentation
+ for( Int_t n=0; n < eNeg.GetEntriesFast(); n++){
+ TLorentzVector en = (*(TLorentzVector*)(eNeg.At(n)));
+ for (Int_t p=0; p < ePos.GetEntriesFast(); p++){
+ TLorentzVector ep = (*(TLorentzVector*)(ePos.At(p)));
+ TLorentzVector np = ep + en;
+ fHistograms->FillHistogram(histoName.Data(),np.M());
+ }
+ }
+}
+
+void AliAnalysisTaskGammaConversion::FillGammaElectronInvMass(TString histoMass,TString histoDiff,TClonesArray const fKFGammas,
+ TClonesArray const tlVeNeg,TClonesArray const tlVePos)
+{
+ //see header file for documentation
+
+ for( Int_t iNeg=0; iNeg < tlVeNeg.GetEntriesFast(); iNeg++ ){
+
+ for (Int_t iPos=0; iPos < tlVePos.GetEntriesFast(); iPos++){
+
+ TLorentzVector xy = *((TLorentzVector *)(tlVePos.At(iPos))) + *((TLorentzVector *)(tlVeNeg.At(iNeg)));
+
+ for (Int_t iGam=0; iGam < fKFGammas.GetEntriesFast(); iGam++){
+
+ // AliKFParticle * gammaCandidate = &fKFGammas[iGam];
+ AliKFParticle * gammaCandidate = (AliKFParticle *)(fKFGammas.At(iGam));
+ TLorentzVector g;
+
+ g.SetXYZM(gammaCandidate->GetPx(),gammaCandidate->GetPy(),gammaCandidate->GetPz(),fGammaMass);
+ TLorentzVector xyg = xy + g;
+ fHistograms->FillHistogram(histoMass.Data(),xyg.M());
+ fHistograms->FillHistogram(histoDiff.Data(),(xyg.M()-xy.M()));
+ }
+ }
+ }
+
+}
+void AliAnalysisTaskGammaConversion::ElectronBackground(TString hBg, TClonesArray e)
+{
+ // see header file for documentation
+ for(Int_t i=0; i < e.GetEntriesFast(); i++)
+ {
+ for (Int_t j=i+1; j < e.GetEntriesFast(); j++)
+ {
+ TLorentzVector ee = (*(TLorentzVector*)(e.At(i))) + (*(TLorentzVector*)(e.At(j)));
+
+ fHistograms->FillHistogram(hBg.Data(),ee.M());
+ }
+ }
+}
+
+
+void AliAnalysisTaskGammaConversion::CleanWithAngleCuts(TClonesArray const negativeElectrons,
+ TClonesArray const positiveElectrons,
+ TClonesArray const gammas){
+ // see header file for documentation
+
+ UInt_t sizeN = negativeElectrons.GetEntriesFast();
+ UInt_t sizeP = positiveElectrons.GetEntriesFast();
+ UInt_t sizeG = gammas.GetEntriesFast();
+
+
+
+ vector <Bool_t> xNegBand(sizeN);
+ vector <Bool_t> xPosBand(sizeP);
+ vector <Bool_t> gammaBand(sizeG);
+
+
+ for(UInt_t iNeg=0; iNeg < sizeN; iNeg++) xNegBand[iNeg]=kTRUE;
+ for(UInt_t iPos=0; iPos < sizeP; iPos++) xPosBand[iPos]=kTRUE;
+ for(UInt_t iGam=0; iGam < sizeG; iGam++) gammaBand[iGam]=kTRUE;
+
+
+ for(UInt_t iPos=0; iPos < sizeP; iPos++){
+
+ Double_t aP[3];
+ ((AliESDtrack*)(positiveElectrons.At(iPos)))->GetConstrainedPxPyPz(aP);
+
+ TVector3 ePosV(aP[0],aP[1],aP[2]);
+
+ for(UInt_t iNeg=0; iNeg < sizeN; iNeg++){
+
+ Double_t aN[3];
+ ((AliESDtrack*)(negativeElectrons.At(iNeg)))->GetConstrainedPxPyPz(aN);
+ TVector3 eNegV(aN[0],aN[1],aN[2]);
+
+ if(ePosV.Angle(eNegV) < 0.05){ //e+e- from gamma
+ xPosBand[iPos]=kFALSE;
+ xNegBand[iNeg]=kFALSE;
+ }
+
+ for(UInt_t iGam=0; iGam < sizeG; iGam++){
+ AliKFParticle* gammaCandidate = (AliKFParticle*)gammas.At(iGam);
+ TVector3 gammaCandidateVector(gammaCandidate->Px(),gammaCandidate->Py(),gammaCandidate->Pz());
+ if(ePosV.Angle(gammaCandidateVector) < 0.05 || eNegV.Angle(gammaCandidateVector) < 0.05)
+ gammaBand[iGam]=kFALSE;
+ }
+ }
+ }
+
+
+
+
+ for(UInt_t iPos=0; iPos < sizeP; iPos++){
+ if(xPosBand[iPos]){
+ new((*fCurrentEventPosElectronTClone)[fCurrentEventPosElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(positiveElectrons.At(iPos))));
+ // fCurrentEventPosElectron.push_back(positiveElectrons[iPos]);
+ }
+ }
+ for(UInt_t iNeg=0;iNeg < sizeN; iNeg++){
+ if(xNegBand[iNeg]){
+ new((*fCurrentEventNegElectronTClone)[fCurrentEventNegElectronTClone->GetEntriesFast()]) AliESDtrack((*(AliESDtrack*)(negativeElectrons.At(iNeg))));
+ // fCurrentEventNegElectron.push_back(negativeElectrons[iNeg]);
+ }
+ }
+ for(UInt_t iGam=0; iGam < sizeG; iGam++){
+ if(gammaBand[iGam]){
+ new((*fKFReconstructedGammasCutTClone)[fKFReconstructedGammasCutTClone->GetEntriesFast()]) AliKFParticle((*(AliKFParticle*)(gammas.At(iGam))));
+ //fKFReconstructedGammasCut.push_back(*(AliKFParticle*)gammas->At(iGam));
+ }
+ }
+}
+
+
+void AliAnalysisTaskGammaConversion::GetPID(AliESDtrack *track, Stat_t &pid, Stat_t &weight)
+{
+ // see header file for documentation
+ pid = -1;
+ weight = -1;
+
+ double wpart[5];
+ double wpartbayes[5];
+
+ //get probability of the diffenrent particle types
+ track->GetESDpid(wpart);
+
+ // Tentative particle type "concentrations"
+ double c[5]={0.01, 0.01, 0.85, 0.10, 0.05};
+
+ //Bayes' formula
+ double rcc = 0.;
+ for (int i = 0; i < 5; i++)
+ {
+ rcc+=(c[i] * wpart[i]);
+ }
+
+
+
+ for (int i=0; i<5; i++) {
+ 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)
+ wpartbayes[i] = c[i] * wpart[i] / rcc;
+ }
+ }
+
+
+
+ Float_t max=0.;
+ int ipid=-1;
+ //find most probable particle in ESD pid
+ //0:Electron - 1:Muon - 2:Pion - 3:Kaon - 4:Proton
+ for (int i = 0; i < 5; i++)
+ {
+ if (wpartbayes[i] > max)
+ {
+ ipid = i;
+ max = wpartbayes[i];
+ }
+ }
+
+ pid = ipid;
+ weight = max;
+}
+double AliAnalysisTaskGammaConversion::GetSigmaToVertex(AliESDtrack* t)
+{
+ // Calculates the number of sigma to the vertex.
+
+ Float_t b[2];
+ Float_t bRes[2];
+ Float_t bCov[3];
+ t->GetImpactParameters(b,bCov);
+ if (bCov[0]<=0 || bCov[2]<=0) {
+ AliDebug(1, "Estimated b resolution lower or equal zero!");
+ bCov[0]=0; bCov[2]=0;
+ }
+ bRes[0] = TMath::Sqrt(bCov[0]);
+ bRes[1] = TMath::Sqrt(bCov[2]);
+
+ // -----------------------------------
+ // How to get to a n-sigma cut?
+ //
+ // The accumulated statistics from 0 to d is
+ //
+ // -> Erf(d/Sqrt(2)) for a 1-dim gauss (d = n_sigma)
+ // -> 1 - Exp(-d**2) for a 2-dim gauss (d*d = dx*dx + dy*dy != n_sigma)
+ //
+ // It means that for a 2-dim gauss: n_sigma(d) = Sqrt(2)*ErfInv(1 - Exp((-x**2)/2)
+ // Can this be expressed in a different way?
+
+ if (bRes[0] == 0 || bRes[1] ==0)
+ return -1;
+
+ double d = TMath::Sqrt(TMath::Power(b[0]/bRes[0],2) + TMath::Power(b[1]/bRes[1],2));
+
+ // stupid rounding problem screws up everything:
+ // if d is too big, TMath::Exp(...) gets 0, and TMath::ErfInverse(1) that should be infinite, gets 0 :(
+ if (TMath::Exp(-d * d / 2) < 1e-10)
+ return 1000;
+
+
+ d = TMath::ErfInverse(1 - TMath::Exp(-d * d / 2)) * TMath::Sqrt(2);
+ return d;
+}
+
+//vector <TLorentzVector> AliAnalysisTaskGammaConversion::GetTLorentzVector(vector <AliESDtrack*> esdTrack){
+TClonesArray AliAnalysisTaskGammaConversion::GetTLorentzVector(TClonesArray *const esdTrack){
+ //Return TLoresntz vector of track?
+ // vector <TLorentzVector> tlVtrack(0);
+ TClonesArray array("TLorentzVector",0);
+
+ for(Int_t itrack=0; itrack < esdTrack->GetEntriesFast(); itrack++){
+ double p[3];
+ //esdTrack[itrack]->GetConstrainedPxPyPz(p);
+ ((AliESDtrack*)(esdTrack->At(itrack)))->GetConstrainedPxPyPz(p);
+ TLorentzVector currentTrack;
+ currentTrack.SetXYZM(p[0],p[1],p[2],fElectronMass);
+ new((array)[array.GetEntriesFast()]) TLorentzVector(currentTrack);
+ // tlVtrack.push_back(currentTrack);
+ }
+
+ return array;
+
+ // return tlVtrack;
+}
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff