From de6f809032712511748076a5d9a2d2d7f3eb3ff4 Mon Sep 17 00:00:00 2001 From: kleinb Date: Wed, 15 Apr 2009 14:23:29 +0000 Subject: [PATCH] Added PID task, some fixes for coding conventions --- PWG4/CMake_libPWG4JetTasks.txt | 1 + PWG4/JetTasks/AliAnaESDSpectraQA.cxx | 6 +- PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx | 54 +- PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx | 926 ++++++++++++++++++ PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h | 97 ++ PWG4/JetTasks/AliAnalysisTaskUE.cxx | 30 +- PWG4/JetTasks/AliAnalysisTaskUE.h | 2 +- PWG4/PWG4JetTasksLinkDef.h | 2 +- PWG4/libPWG4JetTasks.pkg | 2 +- 9 files changed, 1081 insertions(+), 39 deletions(-) create mode 100644 PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx create mode 100644 PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h diff --git a/PWG4/CMake_libPWG4JetTasks.txt b/PWG4/CMake_libPWG4JetTasks.txt index 2ad98b59a56..59b1bc08897 100644 --- a/PWG4/CMake_libPWG4JetTasks.txt +++ b/PWG4/CMake_libPWG4JetTasks.txt @@ -5,6 +5,7 @@ set(SRCS JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx + JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx ) diff --git a/PWG4/JetTasks/AliAnaESDSpectraQA.cxx b/PWG4/JetTasks/AliAnaESDSpectraQA.cxx index 303738cb089..ffafff5c6e5 100644 --- a/PWG4/JetTasks/AliAnaESDSpectraQA.cxx +++ b/PWG4/JetTasks/AliAnaESDSpectraQA.cxx @@ -176,14 +176,14 @@ void AliAnaESDSpectraQA::Exec(Option_t *) { Int_t nTracks = fESD->GetNumberOfTracks(); AliDebug(2,Form("nTracks %d", nTracks)); printf("nTracks %d\n", nTracks); - static Int_t Mult = 0; - Mult = 0; // Need extra init bc of static + static Int_t fMult = 0; + fMult = 0; // Need extra init bc of static for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) { AliESDtrack *track = fESD->GetTrack(iTrack); hists *curTypeHists = 0; if (fTrackCuts->AcceptTrack(track)) { - Mult++; + fMult++; Float_t dca2D, dcaZ; track->GetImpactParameters(dca2D,dcaZ); diff --git a/PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx b/PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx index a9da2688dc7..5a2947d4b7a 100644 --- a/PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx @@ -260,7 +260,7 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects() binLimitsPt[iPt] = 0.0; } else {// 1.0 - binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2; + binLimitsPt[iPt] = binLimitsPt[iPt-1] + 2.5; } } @@ -316,22 +316,22 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects() fh2PtFGen[ij] = new TH2F(Form("fh2PtFGen_j%d",ij),"Pt Found vs. gen;p_{T,rec} (GeV/c);p_{T,gen} (GeV/c)", nBinPt,binLimitsPt,nBinPt,binLimitsPt); - fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};phi_{gen}", + fh2PhiFGen[ij] = new TH2F(Form("fh2PhiFGen_j%d",ij),"#phi Found vs. gen;#phi_{rec};#phi_{gen}", nBinPhi,binLimitsPhi,nBinPhi,binLimitsPhi); - fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};eta_{gen}", + fh2EtaFGen[ij] = new TH2F(Form("fh2EtaFGen_j%d",ij),"#eta Found vs. gen;#eta_{rec};#eta_{gen}", nBinEta,binLimitsEta,nBinEta,binLimitsEta); - fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-phi_{rec}", + fh2PtGenDeltaPhi[ij] = new TH2F(Form("fh2PtGenDeltaPhi_j%d",ij),"delta phi vs. P_{T,gen};p_{T,gen} (GeV/c);#phi_{gen}-#phi_{rec}", nBinPt,binLimitsPt,100,-1.0,1.0); - fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-eta_{rec}", + fh2PtGenDeltaEta[ij] = new TH2F(Form("fh2PtGenDeltaEta_j%d",ij),"delta eta vs. p_{T,gen};p_{T,gen} (GeV/c);#eta_{gen}-#eta_{rec}", nBinPt,binLimitsPt,100,-1.0,1.0); - fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#Delta{R} to lower energy jets j > i;p_{T,rec,j};#Delta R", + fh2PtRecDeltaR[ij] = new TH2F(Form("fh2PtRecDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,rec,j};#Delta R", nBinPt,binLimitsPt,60,0,6.0); - fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#Delta{R} to lower energy jets j > i;p_{T,gen,j};#Delta R", + fh2PtGenDeltaR[ij] = new TH2F(Form("fh2PtGenDeltaR_j%d",ij),"#DeltaR to lower energy jets j > i;p_{T,gen,j};#Delta R", nBinPt,binLimitsPt,60,0,6.0); @@ -363,7 +363,7 @@ void AliAnalysisTaskJetSpectrum::UserCreateOutputObjects() - fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,rec} (GeV/c)", + fh3GenEtaPhiPt[ij] = new TH3F(Form("fh3GenEtaPhiPt_j%d",ij),"Gen eta, phi, pt; #eta; #phi; p_{T,} (GeV/c)", nBinEta,binLimitsEta,nBinPhi,binLimitsPhi,nBinPt,binLimitsPt); } @@ -532,7 +532,23 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/) nTrials = pythiaGenHeader->Trials(); ptHard = pythiaGenHeader->GetPtHard(); - + int iProcessType = pythiaGenHeader->ProcessType(); + // 11 f+f -> f+f + // 12 f+barf -> f+barf + // 13 f+barf -> g+g + // 28 f+g -> f+g + // 53 g+g -> f+barf + // 68 g+g -> g+g + /* + if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType); + // if(iProcessType != 13 && iProcessType != 68){ // allow only glue + if(iProcessType != 11 && iProcessType != 12 && iProcessType != 53){ // allow only quark + // if(iProcessType != 28){ // allow only -> f+g + PostData(1, fHistList); + return; + } + */ + if (fDebug > 10)Printf("%d iProcessType %d",__LINE__, iProcessType); if(fDebug>20)AliAnalysisHelperJetTasks::PrintStack(mcEvent); @@ -675,10 +691,14 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/) fh3PtRecGenHard[ir]->Fill(ptRec,ptGen,ptHard,eventW); fh3PtRecGenHardNoW[ir]->Fill(ptRec,ptGen,ptHard,1); ///////////////////////////////////////////////////// + + // Double_t eRec = recJets[ir].E(); + // Double_t eGen = genJets[ig].E(); + // CKB use p_T not Energy + // TODO recname variabeles and histos Double_t eRec = recJets[ir].E(); Double_t eGen = genJets[ig].E(); - fh2Efficiency->Fill(eGen, eRec/eGen); if (eGen>=0. && eGen<=250.){ @@ -688,19 +708,17 @@ void AliAnalysisTaskJetSpectrum::UserExec(Option_t */*option*/) // loop over tracks for (Int_t it = 0; it< nTracks; it++){ // if (fAOD->GetTrack(it)->E() > eGen) continue; // CKB. Not allowed! cannot cut on gen properties in real events! - Double_t phiTrack = fAOD->GetTrack(it)->Phi(); - if (phiTrack<0) phiTrack+=TMath::Pi()*2.; - Double_t etaTrack = fAOD->GetTrack(it)->Eta(); - Float_t deta = etaRec - etaTrack; - Float_t dphi = TMath::Abs(phiRec - phiTrack); - Float_t r = TMath::Sqrt(deta*deta + dphi*dphi); // find leading particle - if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){ + // if (r<0.4 && fAOD->GetTrack(it)->E()>eLeading){ + // TODO implement esd filter flag to be the same as in the jet finder + // allow also for MC particles... + Float_t r = recJets[ir].DeltaR(fAOD->GetTrack(it)); + if (r<0.4 && fAOD->GetTrack(it)->Pt()>ptleading){ eLeading = fAOD->GetTrack(it)->E(); ptleading = fAOD->GetTrack(it)->Pt(); } // if (fAOD->GetTrack(it)->Pt()>0.03*eGen && fAOD->GetTrack(it)->E()<=eGen && r<0.7) // CKB cannot cut on gen properties - if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->E()<=eRec && r<0.7) + if (fAOD->GetTrack(it)->Pt()>0.03*eRec && fAOD->GetTrack(it)->Pt()<=eRec && r<0.7) nPart++; } if (fDebug > 10)Printf("%s:%d",(char*)__FILE__,__LINE__); diff --git a/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx b/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx new file mode 100644 index 00000000000..f43fe918732 --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx @@ -0,0 +1,926 @@ +/************************************************************************** + * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * + * * + * Author: The ALICE Off-line Project. * + * Contributors are mentioned in the code where appropriate. * + * * + * 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. * + **************************************************************************/ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "AliAnalysisManager.h" +#include "AliAnalysisFilter.h" +#include "AliESDInputHandler.h" +#include "AliESDEvent.h" +#include "AliESDVertex.h" +#include "AliMCEventHandler.h" +#include "AliMCEvent.h" +#include "AliStack.h" +#include "AliAODInputHandler.h" +#include "AliAODEvent.h" +#include "AliAODVertex.h" +#include "AliAODMCParticle.h" +#include "AliLog.h" + +#include "AliAnalysisTaskPWG4PidDetEx.h" + +// STL includes +#include +using namespace std; + +// +// Analysis class example for PID using detector signals. +// Works on ESD and AOD; has a flag for MC analysis +// +// Alexandru.Dobrin@hep.lu.se +// Peter.Christiansen@hep.lu.se +// + +ClassImp(AliAnalysisTaskPWG4PidDetEx) + +const Double_t AliAnalysisTaskPWG4PidDetEx::fgkCtau = 370; // distance for kaon decay +const Double_t AliAnalysisTaskPWG4PidDetEx::fgkPionMass = TDatabasePDG::Instance()->GetParticle(211)->Mass(); +const Double_t AliAnalysisTaskPWG4PidDetEx::fgkKaonMass = TDatabasePDG::Instance()->GetParticle(321)->Mass(); +const Double_t AliAnalysisTaskPWG4PidDetEx::fgkProtonMass = TDatabasePDG::Instance()->GetParticle(2212)->Mass(); + +//_____________________________________________________________________________ +AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx(): + AliAnalysisTaskSE(), + fESD(0), + fAOD(0), + fListOfHists(0), + fEtaCut(0.9), + fPtCut(0.5), + fXbins(20), + fXmin(0.), + fTOFCutP(2.), + fTrackFilter(0x0), + fAnalysisMC(kTRUE), + fAnalysisType("ESD"), + fTriggerMode(kMB2), + fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0), + fdNdPt(0), fMassAll(0), + fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0), + fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0), + fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0), + fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0) +{ + // Default constructor (should not be used) +} + +//______________________________________________________________________________ +AliAnalysisTaskPWG4PidDetEx::AliAnalysisTaskPWG4PidDetEx(const char *name): + AliAnalysisTaskSE(name), + fESD(0), + fAOD(0), + fListOfHists(0), + fEtaCut(0.9), + fPtCut(0.5), + fXbins(20), + fXmin(0.), + fTOFCutP(2.), + fTrackFilter(0x0), + fAnalysisMC(kTRUE), + fAnalysisType("ESD"), + fTriggerMode(kMB2), + fEvents(0), fEffTot(0), fEffPID(0), fAccP(0), fAccPt(0), fKaonDecayCorr(0), + fdNdPt(0), fMassAll(0), + fdNdPtPion(0), fMassPion(0), fdEdxTPCPion(0), fbgTPCPion(0), + fdNdPtKaon(0), fMassKaon(0), fdEdxTPCKaon(0), fbgTPCKaon(0), + fdNdPtProton(0), fMassProton(0), fdEdxTPCProton(0), fbgTPCProton(0), + fdNdPtMC(0), fdNdPtMCPion(0), fdNdPtMCKaon(0), fdNdPtMCProton(0) +{ + // Normal constructor + // Input slot #0 works with a TChain + DefineInput(0, TChain::Class()); + // Output slot #0 writes into a TList + DefineOutput(0, TList::Class()); + +} + +//_____________________________________________________________________________ +AliAnalysisTaskPWG4PidDetEx::~AliAnalysisTaskPWG4PidDetEx() +{ + // Destructor + // histograms are in the output list and deleted when the output + // list is deleted by the TSelector dtor + if (fListOfHists) { + delete fListOfHists; + fListOfHists = 0; + } +} + +//______________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::ConnectInputData(Option_t *) +{ + // Connect AOD here + // Called once + if (fDebug > 1) AliInfo("ConnectInputData() \n"); + + TTree* tree = dynamic_cast (GetInputData(0)); + if (!tree) { + Printf("ERROR: Could not read chain from input slot 0"); + } + else { + if(fAnalysisType == "ESD") { + AliESDInputHandler *esdH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + if (!esdH) { + Printf("ERROR: Could not get ESDInputHandler"); + } else + fESD = esdH->GetEvent(); + } + else if(fAnalysisType == "AOD") { + AliAODInputHandler *aodH = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()); + if (!aodH) { + Printf("ERROR: Could not get AODInputHandler"); + } else + fAOD = aodH->GetEvent(); + } + } + +} + +//______________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::CreateOutputObjects() +{ + // Create histograms + // Called once + if (fDebug > 1) AliInfo("CreateOutPutData() \n"); + + OpenFile(0); + fListOfHists = new TList(); + + fEvents = new TH1I("fEvents","Number of analyzed events; Events; Counts", 1, 0, 1); + fListOfHists->Add(fEvents); + + fEffTot = new TH1F ("fEffTot","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP); + fEffTot->Sumw2(); + fListOfHists->Add(fEffTot); + + fEffPID = new TH1F ("fEffPID","Efficiency; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP); + fEffPID->Sumw2(); + fListOfHists->Add(fEffPID); + + fAccP = new TH1F ("fAccP","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP); + fAccP->Sumw2(); + fListOfHists->Add(fAccP); + + fAccPt = new TH1F ("fAccPt","Acceptance; p_{T} [GeV/c]; Counts", fXbins, fXmin, fTOFCutP); + fAccPt->Sumw2(); + fListOfHists->Add(fAccPt); + + fKaonDecayCorr = new TProfile("fKaonDecayCorr","Kaon decay correction vs p_{T}; p_{T} [GeV/c]; Kaon decay correction", fXbins, fXmin, fTOFCutP); + fListOfHists->Add(fKaonDecayCorr); + + fdNdPt = new TH1F("fdNdPt","p_{T} distribution; p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100 , 0, 10); + fdNdPt->Sumw2(); + fListOfHists->Add(fdNdPt); + + fMassAll = new TH2F("fMassAll","Mass^{2} vs momentum (ToF); p [GeV/c]; M^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2); + fListOfHists->Add(fMassAll); + + //Pion + fdNdPtPion = new TH1F("fdNdPtPion","p_{T} distribution (Pions); p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtPion->Sumw2(); + fListOfHists->Add(fdNdPtPion); + + fMassPion = new TH2F("fMassPion","Mass squared for pions after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2); + fListOfHists->Add(fMassPion); + + fdEdxTPCPion = new TH2F("fdEdxTPCPion","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200); + fListOfHists->Add(fdEdxTPCPion); + + fbgTPCPion = new TH2F("fbgTPCPion", "Energy loss vs #beta#gamma (Pions);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200); + fListOfHists->Add(fbgTPCPion); + + //Kaon + fdNdPtKaon = new TH1F("fdNdPtKaon","p_{T} distribution (Kaons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtKaon->Sumw2(); + fListOfHists->Add(fdNdPtKaon); + + fMassKaon = new TH2F("fMassKaon","Mass squared for kaons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2); + fListOfHists->Add(fMassKaon); + + fdEdxTPCKaon = new TH2F("fdEdxTPCKaon","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200); + fListOfHists->Add(fdEdxTPCKaon); + + fbgTPCKaon = new TH2F("fbgTPCKaon", "Energy loss vs #beta#gamma (Kaons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200); + fListOfHists->Add(fbgTPCKaon); + + //Proton + fdNdPtProton = new TH1F("fdNdPtProton","p_{T} distribution (Protons);p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtProton->Sumw2(); + fListOfHists->Add(fdNdPtProton); + + fMassProton = new TH2F("fMassProton","Mass squared for protons after cuts vs pmean; p [GeV/c]; Mass^{2} [GeV^{2}/c^{4}]", 100, -5, 5, 100, -0.2, 1.2); + fListOfHists->Add(fMassProton); + + fdEdxTPCProton = new TH2F("fdEdxTPCProton","Energy loss vs momentum; p [GeV/c]; dE/dx [a. u.]", 40, -2, 2, 100, 0, 200); + fListOfHists->Add(fdEdxTPCProton); + + fbgTPCProton = new TH2F("fbgTPCProton", "Energy loss vs #beta#gamma (Protons);#beta#gamma; dE/dx [a. u.]", 100, 0, 15, 100, 0, 200); + fListOfHists->Add(fbgTPCProton); + + + //MC + if(fAnalysisMC){ + fdNdPtMC = new TH1F("fdNdPtMC","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", 100, 0, 10); + fdNdPtMC->SetLineColor(2); + fdNdPtMC->Sumw2(); + fListOfHists->Add(fdNdPtMC); + + fdNdPtMCPion = new TH1F("fdNdPtMCPion","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtMCPion->SetLineColor(2); + fdNdPtMCPion->Sumw2(); + fListOfHists->Add(fdNdPtMCPion); + + fdNdPtMCKaon = new TH1F("fdNdPtMCKaon","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtMCKaon->SetLineColor(2); + fdNdPtMCKaon->Sumw2(); + fListOfHists->Add(fdNdPtMCKaon); + + fdNdPtMCProton = new TH1F("fdNdPtMCProton","p_{T} distribution;p_{T} [GeV/c]; #frac{1}{2#pip_{T}} #frac{1}{N_{ev}} #frac{dN}{dp_{T}}", fXbins, fXmin, fTOFCutP); + fdNdPtMCProton->SetLineColor(2); + fdNdPtMCProton->Sumw2(); + fListOfHists->Add(fdNdPtMCProton); + } + +} + +//______________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::Exec(Option_t *) +{ + // Main loop + // Called for each event + if (fDebug > 1) AliInfo("Exec() \n" ); + + if(fAnalysisType == "ESD") { + if (!fESD) { + Printf("ERROR: fESD not available"); + return; + } + if(IsEventTriggered(fESD, fTriggerMode)) { + Printf("PWG4 PID ESD analysis"); + AnalyzeESD(fESD); + }//triggered event + }//ESD analysis + + else if(fAnalysisType == "AOD") { + if (!fAOD) { + Printf("ERROR: fAOD not available"); + return; + } + if(IsEventTriggered(fAOD, fTriggerMode)) { + Printf("PWG4 PID AOD analysis"); + AnalyzeAOD(fAOD); + }//triggered event + }//AOD analysis + + // Post output data. + PostData(0, fListOfHists); +} + +//________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::AnalyzeESD(AliESDEvent* esd) +{ + // Get vertex + const AliESDVertex* primaryVertex = esd->GetPrimaryVertex(); + const Double_t vertexResZ = primaryVertex->GetZRes(); + + // Select only events with a reconstructed vertex + if(vertexResZ<5.0) { + const Int_t nESDTracks = esd->GetNumberOfTracks(); + for(Int_t i = 0; i < nESDTracks; i++) { + AliESDtrack* trackESD = esd->GetTrack(i); + + //Cuts + UInt_t selectInfo = 0; + if (fTrackFilter) { + selectInfo = fTrackFilter->IsSelected(trackESD); + if (!selectInfo) continue; + } + + if ((trackESD->Pt() < fPtCut) || (TMath::Abs(trackESD->Eta()) > fEtaCut )) + continue; + + //Corrections + fEffTot->Fill(trackESD->Pt()); + if (trackESD->GetTOFsignal() !=0) + fEffPID->Fill(trackESD->Pt()); + + if (trackESD->P() < fTOFCutP) fAccP->Fill(trackESD->Pt()); + if (trackESD->Pt() < fTOFCutP) fAccPt->Fill(trackESD->Pt()); + + //Analysis + fdNdPt->Fill(trackESD->Pt(), 1.0/trackESD->Pt()); + + //TOF + if ((trackESD->GetIntegratedLength() == 0) || (trackESD->GetTOFsignal() == 0)) + continue; + + fMassAll->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD)); + + if ((MassSquared(trackESD) < 0.15) && (MassSquared(trackESD) > -0.15) && (trackESD->P() < fTOFCutP)){ + fdNdPtPion->Fill(trackESD->Pt(), 1.0/trackESD->Pt()); + fMassPion->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD)); + fdEdxTPCPion->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal()); + fbgTPCPion->Fill(trackESD->P()/fgkPionMass,trackESD->GetTPCsignal()); + } + + if ((MassSquared(trackESD) > 0.15) && (MassSquared(trackESD) < 0.45) && (trackESD->P() < fTOFCutP)){ + fdNdPtKaon->Fill(trackESD->Pt(), 1.0/trackESD->Pt()); + fMassKaon->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD)); + fdEdxTPCKaon->Fill(trackESD->Charge()*trackESD->P(), trackESD->GetTPCsignal()); + fbgTPCKaon->Fill(trackESD->P()/fgkKaonMass, trackESD->GetTPCsignal()); + //Kaon decay correction + fKaonDecayCorr->Fill(trackESD->Pt(), TMath::Exp(KaonDecay(trackESD))); + } + + if ((MassSquared(trackESD) > 0.75) && (MassSquared(trackESD) < 1.05) && (trackESD->P() < fTOFCutP)){ + fdNdPtProton->Fill(trackESD->Pt(), 1.0/trackESD->Pt()); + fMassProton->Fill(trackESD->Charge()*trackESD->P(), MassSquared(trackESD)); + fdEdxTPCProton->Fill(trackESD->Charge()*trackESD->P(),trackESD->GetTPCsignal()); + fbgTPCProton->Fill(trackESD->P()/fgkProtonMass,trackESD->GetTPCsignal()); + } + }//ESD track loop + + //MC + if(fAnalysisMC){ + AliMCEventHandler* mcHandler = dynamic_cast (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); + if (!mcHandler) { + Printf("ERROR: Could not retrieve MC event handler"); + return; + } + + AliMCEvent* mcEvent = mcHandler->MCEvent(); + if (!mcEvent) { + Printf("ERROR: Could not retrieve MC event"); + return; + } + + AliStack* mcStack = mcEvent->Stack(); + if (!mcStack) { + Printf("ERROR: Could not retrieve MC stack"); + return; + } + + const Int_t nTracksMC = mcStack->GetNtrack(); + for (Int_t iTracks = 0; iTracks < nTracksMC; iTracks++) { + //Cuts + if(!(mcStack->IsPhysicalPrimary(iTracks))) + continue; + + TParticle* mcTrack = mcStack->Particle(iTracks); + + Double_t charge = mcTrack->GetPDG()->Charge(); + if (charge == 0) + continue; + + if ((mcTrack->Pt() < fPtCut) || (TMath::Abs(mcTrack->Eta()) > fEtaCut )) + continue; + + //Analysis + fdNdPtMC->Fill(mcTrack->Pt(), 1.0/mcTrack->Pt()); + + if ((mcTrack->GetPdgCode() == 211) || (mcTrack->GetPdgCode() == -211)) + fdNdPtMCPion->Fill(mcTrack->Pt(),1.0/mcTrack->Pt()); + + if ((mcTrack->GetPdgCode() == 321) || (mcTrack->GetPdgCode() == -321)) + fdNdPtMCKaon->Fill(mcTrack->Pt(),1.0/mcTrack->Pt()); + + if ((mcTrack->GetPdgCode() == 2212) || (mcTrack->GetPdgCode() == -2212)) + fdNdPtMCProton->Fill(mcTrack->Pt(),1.0/mcTrack->Pt()); + }//MC track loop + }//if MC + fEvents->Fill(0); + }//if vertex + +} + +//________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::AnalyzeAOD(AliAODEvent* aod) +{ + // Get vertex + AliAODVertex* primaryVertex = aod->GetPrimaryVertex(); + const Double_t vertexResZ = TMath::Sqrt(primaryVertex->RotatedCovMatrixZZ()); + + // Select only events with a reconstructed vertex + if (vertexResZ < 5) { + //AOD + Int_t nGoodTracks = aod->GetNumberOfTracks(); + + for(Int_t iTracks = 0; iTracks < nGoodTracks; iTracks++) { + AliAODTrack* trackAOD = aod->GetTrack(iTracks); + //Cuts + if (!(trackAOD->IsPrimaryCandidate())) + continue; + + if ((trackAOD->Pt() < fPtCut) || (TMath::Abs(trackAOD->Eta()) > fEtaCut )) + continue; + + //Corrections + fEffTot->Fill(trackAOD->Pt()); + if (trackAOD->GetDetPid()->GetTOFsignal() !=0 ) + fEffPID->Fill(trackAOD->Pt()); + + if (trackAOD->P() < fTOFCutP) fAccP->Fill(trackAOD->Pt()); + if (trackAOD->Pt() < fTOFCutP) fAccPt->Fill(trackAOD->Pt()); + + //Analysis + fdNdPt->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt()); + + //TOF + if ((IntegratedLength(trackAOD) == 0) || (trackAOD->GetDetPid()->GetTOFsignal() == 0)) + continue; + + fMassAll->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD)); + + if ((MassSquared(trackAOD) < 0.15) && (MassSquared(trackAOD) > -0.15) && (trackAOD->P() < fTOFCutP)){ + fdNdPtPion->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt()); + fMassPion->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD)); + fdEdxTPCPion->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal()); + fbgTPCPion->Fill(trackAOD->P()/fgkPionMass,trackAOD->GetDetPid()->GetTPCsignal()); + } + + if ((MassSquared(trackAOD) > 0.15) && (MassSquared(trackAOD) < 0.45) && (trackAOD->P() < fTOFCutP)){ + fdNdPtKaon->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt()); + fMassKaon->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD)); + fdEdxTPCKaon->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal()); + fbgTPCKaon->Fill(trackAOD->P()/fgkKaonMass,trackAOD->GetDetPid()->GetTPCsignal()); + //Kaon decay correction + fKaonDecayCorr->Fill(trackAOD->Pt(), TMath::Exp(KaonDecay(trackAOD))); + } + + if ((MassSquared(trackAOD) > 0.75) && (MassSquared(trackAOD) < 1.05) && (trackAOD->P() < fTOFCutP)){ + fdNdPtProton->Fill(trackAOD->Pt(), 1.0/trackAOD->Pt()); + fMassProton->Fill(trackAOD->Charge()*trackAOD->P(), MassSquared(trackAOD)); + fdEdxTPCProton->Fill(trackAOD->Charge()*trackAOD->P(),trackAOD->GetDetPid()->GetTPCsignal()); + fbgTPCProton->Fill(trackAOD->P()/fgkProtonMass,trackAOD->GetDetPid()->GetTPCsignal()); + } + }//AOD track loop + + //MC + if(fAnalysisMC){ + TClonesArray* farray = (TClonesArray*)aod->FindListObject("mcparticles"); + Int_t ntrks = farray->GetEntries(); + for(Int_t i =0 ; i < ntrks; i++){ + AliAODMCParticle* trk = (AliAODMCParticle*)farray->At(i); + //Cuts + if (!(trk->IsPhysicalPrimary())) + continue; + + if (trk->Charge() == 0) + continue; + + if ((trk->Pt() < fPtCut) || (TMath::Abs(trk->Eta()) > fEtaCut )) + continue; + + //Analysis + fdNdPtMC->Fill(trk->Pt(), 1.0/trk->Pt()); + + if ((trk->GetPdgCode() == 211) || (trk->GetPdgCode() == -211)) + fdNdPtMCPion->Fill(trk->Pt(),1.0/trk->Pt()); + + if ((trk->GetPdgCode() == 321) || (trk->GetPdgCode() == -321)) + fdNdPtMCKaon->Fill(trk->Pt(),1.0/trk->Pt()); + + if ((trk->GetPdgCode() == 2212) || (trk->GetPdgCode() == -2212)) + fdNdPtMCProton->Fill(trk->Pt(),1.0/trk->Pt()); + }//MC track loop + }//if MC + fEvents->Fill(0); + }//if vertex resolution + +} + +//________________________________________________________________________ +Bool_t AliAnalysisTaskPWG4PidDetEx::IsEventTriggered(AliVEvent* ev, TriggerMode trigger) +{ + //adapted from PWG2 (AliAnalysisTaskProton) + + ULong64_t triggerMask = ev->GetTriggerMask(); + + // definitions from p-p.cfg + ULong64_t spdFO = (1 << 14); + ULong64_t v0left = (1 << 11); + ULong64_t v0right = (1 << 12); + + switch (trigger) { + case kMB1: + if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right))) + return kTRUE; + break; + + case kMB2: + if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right))) + return kTRUE; + break; + + case kSPDFASTOR: + if (triggerMask & spdFO) + return kTRUE; + break; + + }//switch + + return kFALSE; +} + +//_____________________________________________________________________________ +Double_t AliAnalysisTaskPWG4PidDetEx::IntegratedLength(AliVTrack* track) const +{ + Double_t intTime [5]; + for (Int_t i = 0; i < 5; i++) intTime[i] = -100.; + Double_t timeElectron = 0, intLength = 0; + + AliAODTrack* trackAOD = (AliAODTrack*)track; + trackAOD->GetDetPid()->GetIntegratedTimes(intTime); + timeElectron = intTime[0]; + intLength = TMath::C()*timeElectron; + + return intLength; +} + +//_____________________________________________________________________________ +Double_t AliAnalysisTaskPWG4PidDetEx::MassSquared(AliVTrack* track) const +{ + Double_t beta = -10, mass = -10; + + if(fAnalysisType == "ESD"){ + AliESDtrack* trackESD = (AliESDtrack*)track; + beta = trackESD->GetIntegratedLength()/trackESD->GetTOFsignal()/TMath::C(); + mass = trackESD->P()*trackESD->P()*(1./(beta*beta) - 1.0); + } + + if(fAnalysisType == "AOD"){ + AliAODTrack* trackAOD = (AliAODTrack*)track; + beta =IntegratedLength(trackAOD)/trackAOD->GetDetPid()->GetTOFsignal()/TMath::C(); + mass = trackAOD->P()*trackAOD->P()*(1./(beta*beta) - 1.0); + } + + return mass; +} + +//_____________________________________________________________________________ +Double_t AliAnalysisTaskPWG4PidDetEx::KaonDecay(AliVTrack* track) const +{ + Double_t decay = -10; + + if(fAnalysisType == "ESD"){ + AliESDtrack* trackESD = (AliESDtrack*)track; + decay = trackESD->GetIntegratedLength()*fgkKaonMass/fgkCtau/trackESD->P(); + } + + if (fAnalysisType == "AOD"){ + AliAODTrack* trackAOD = (AliAODTrack*)track; + decay = IntegratedLength(trackAOD)*fgkKaonMass/fgkCtau/trackAOD->P(); + } + + return decay; +} + +//_____________________________________________________________________________ +void AliAnalysisTaskPWG4PidDetEx::Terminate(Option_t *) +{ + // Terminate loop + if (fDebug > 1) Printf("Terminate()"); + + // + // The followig code has now been moved to drawPid.C + // + +// fListOfHists = dynamic_cast (GetOutputData(0)); +// if (!fListOfHists) { +// Printf("ERROR: fListOfHists not available"); +// return; +// } + +// TH1I* hevents = dynamic_cast (fListOfHists->FindObject("fEvents")); +// Int_t nEvents = (Int_t) hevents->GetBinContent(1); + +// const Float_t normalization = 2.0*fEtaCut*2.0*TMath::Pi(); + +// //----------------- +// TCanvas* c1 = new TCanvas("c1", "c1"); +// c1->cd(); + +// TH2F* hMassAll = dynamic_cast (fListOfHists->FindObject("fMassAll")); +// hMassAll->SetLineColor(1); +// hMassAll->SetMarkerColor(1); +// hMassAll->DrawCopy(); + +// TH2F* hMassPion = dynamic_cast (fListOfHists->FindObject("fMassPion")); +// hMassPion->SetLineColor(2); +// hMassPion->SetMarkerColor(2); +// hMassPion->SetMarkerStyle(7); +// hMassPion->DrawCopy("same"); + +// TH2F* hMassKaon = dynamic_cast (fListOfHists->FindObject("fMassKaon")); +// hMassKaon->SetLineColor(3); +// hMassKaon->SetMarkerColor(3); +// hMassKaon->SetMarkerStyle(7); +// hMassKaon->DrawCopy("same"); + +// TH2F* hMassProton = dynamic_cast (fListOfHists->FindObject("fMassProton")); +// hMassProton->SetLineColor(4); +// hMassProton->SetMarkerColor(4); +// hMassProton->SetMarkerStyle(7); +// hMassProton->DrawCopy("same"); + +// TLegend* legend1 = new TLegend(0.8, 0.8, 1, 1); +// legend1->SetBorderSize(0); +// legend1->SetFillColor(0); +// legend1->AddEntry(hMassAll, "All","LP"); +// legend1->AddEntry(hMassPion, "Pions","LP"); +// legend1->AddEntry(hMassKaon, "Kaons","LP"); +// legend1->AddEntry(hMassProton, "Protons","LP"); +// legend1->Draw(); + +// c1->Update(); + +// //----------------- +// TCanvas* c2 = new TCanvas("c2", "c2"); +// c2->cd(); + +// TH2F* hdEdxTPCPion = dynamic_cast (fListOfHists->FindObject("fdEdxTPCPion")); +// hdEdxTPCPion->SetTitle("dE/dx vs p (TPC)"); +// hdEdxTPCPion->SetLineColor(2); +// hdEdxTPCPion->SetMarkerColor(2); +// hdEdxTPCPion->SetMarkerStyle(7); +// hdEdxTPCPion->DrawCopy(); + +// TH2F* hdEdxTPCKaon = dynamic_cast (fListOfHists->FindObject("fdEdxTPCKaon")); +// hdEdxTPCKaon->SetLineColor(3); +// hdEdxTPCKaon->SetMarkerColor(3); +// hdEdxTPCKaon->SetMarkerStyle(7); +// hdEdxTPCKaon->DrawCopy("same"); + +// TH2F* hdEdxTPCProton = dynamic_cast (fListOfHists->FindObject("fdEdxTPCProton")); +// hdEdxTPCProton->SetLineColor(4); +// hdEdxTPCProton->SetMarkerColor(4); +// hdEdxTPCProton->SetMarkerStyle(7); +// hdEdxTPCProton->DrawCopy("same"); + +// TLegend* legend2 = new TLegend(0.66, 0.66, 0.88, 0.88); +// legend2->SetBorderSize(0); +// legend2->SetFillColor(0); +// legend2->AddEntry(hdEdxTPCPion, "Pions","LP"); +// legend2->AddEntry(hdEdxTPCKaon, "Kaons","LP"); +// legend2->AddEntry(hdEdxTPCProton, "Protons","LP"); +// legend2->Draw(); + +// c2->Update(); + +// //----------------- +// TCanvas* c3 = new TCanvas("c3", "c3"); +// c3->cd(); + +// TH2F* hdEdxTPCbgPion = dynamic_cast (fListOfHists->FindObject("fbgTPCPion")); +// hdEdxTPCbgPion->SetTitle("dE/dx vs #beta#gamma (TPC)"); +// hdEdxTPCbgPion->SetLineColor(2); +// hdEdxTPCbgPion->SetMarkerColor(2); +// hdEdxTPCbgPion->SetMarkerStyle(7); +// hdEdxTPCbgPion->DrawCopy(); + +// TH2F* hdEdxTPCbgKaon = dynamic_cast (fListOfHists->FindObject("fbgTPCKaon")); +// hdEdxTPCbgKaon->SetLineColor(3); +// hdEdxTPCbgKaon->SetMarkerColor(3); +// hdEdxTPCbgKaon->SetMarkerStyle(7); +// hdEdxTPCbgKaon->DrawCopy("same"); + +// TH2F* hdEdxTPCbgProton = dynamic_cast (fListOfHists->FindObject("fbgTPCProton")); +// hdEdxTPCbgProton->SetLineColor(4); +// hdEdxTPCbgProton->SetMarkerColor(4); +// hdEdxTPCbgProton->SetMarkerStyle(7); +// hdEdxTPCbgProton->DrawCopy("same"); + +// TLegend* legend3 = new TLegend(0.66, 0.66, 0.88, 0.88); +// legend3->SetBorderSize(0); +// legend3->SetFillColor(0); +// legend3->AddEntry(hdEdxTPCbgPion, "Pions","LP"); +// legend3->AddEntry(hdEdxTPCbgKaon, "Kaons","LP"); +// legend3->AddEntry(hdEdxTPCbgProton, "Protons","LP"); +// legend3->Draw(); + +// c3->Update(); + +// //----------------- +// TCanvas* c4 = new TCanvas("c4", "c4", 100, 100, 500, 900); +// c4->Divide(1,3); + +// c4->cd(1); +// TH1F* hAccepatncePt = dynamic_cast (fListOfHists->FindObject("fAccPt")); +// TH1F* hAcceptanceP = dynamic_cast (fListOfHists->FindObject("fAccP")); +// hAccepatncePt->Divide(hAcceptanceP); +// hAccepatncePt->SetTitle("Acceptance correction"); +// hAccepatncePt->GetYaxis()->SetTitle("Acceptance correction"); +// hAccepatncePt->DrawCopy(); + +// c4->cd(2); +// TH1F* hEfficiencyPID = dynamic_cast (fListOfHists->FindObject("fEffPID")); +// TH1F* hEfficiencyTot = dynamic_cast (fListOfHists->FindObject("fEffTot")); +// hEfficiencyPID->Divide(hEfficiencyTot); +// hEfficiencyPID->SetTitle("Efficiency correction"); +// hEfficiencyPID->GetYaxis()->SetTitle("Efficiency correction"); +// hEfficiencyPID->DrawCopy(); + +// c4->cd(3); +// TProfile* hKDecayCorr = dynamic_cast (fListOfHists->FindObject("fKaonDecayCorr")); +// hKDecayCorr->SetTitle("Kaon decay correction"); +// hKDecayCorr->GetYaxis()->SetTitle("Kaon decay correction"); +// hKDecayCorr->GetXaxis()->SetTitle(" p_{T} [GeV/c]"); +// hKDecayCorr->DrawCopy(); + +// c4->Update(); + +// //--------------------- +// TCanvas* c5 = new TCanvas("c5", "c5", 100, 100, 600, 900); +// c5->Divide(1,2); + +// c5->cd(1); +// TH1F* hPtPionNoCorr = dynamic_cast (fListOfHists->FindObject("fdNdPtPion")); +// hPtPionNoCorr->Scale(1.0/nEvents/normalization/hPtPionNoCorr->GetBinWidth(1)); +// hPtPionNoCorr->SetTitle("p_{T} distribution (no corrections)"); +// hPtPionNoCorr->SetLineColor(2); +// hPtPionNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1); +// hPtPionNoCorr->DrawCopy(); + +// TH1F* hPtKaonNoCorr = dynamic_cast (fListOfHists->FindObject("fdNdPtKaon")); +// hPtKaonNoCorr->Scale(1.0/nEvents/normalization/hPtKaonNoCorr->GetBinWidth(1)); +// hPtKaonNoCorr->SetLineColor(3); +// hPtKaonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1); +// hPtKaonNoCorr->DrawCopy("same"); + +// TH1F* hPtProtonNoCorr = dynamic_cast (fListOfHists->FindObject("fdNdPtProton")); +// hPtProtonNoCorr->Scale(1.0/nEvents/normalization/hPtProtonNoCorr->GetBinWidth(1)); +// hPtProtonNoCorr->SetLineColor(4); +// hPtProtonNoCorr->GetYaxis()->SetRangeUser(1E-3,1E1); +// hPtProtonNoCorr->DrawCopy("same"); + +// TH1F* hPt = dynamic_cast (fListOfHists->FindObject("fdNdPt")); +// hPt->Scale(1.0/nEvents/normalization/hPt->GetBinWidth(1)); +// hPt->GetYaxis()->SetRangeUser(1E-3,1E1); +// hPt->DrawCopy("same"); + +// TLegend* legend4 = new TLegend(0.63, 0.63, 0.88, 0.88); +// legend4->SetBorderSize(0); +// legend4->SetFillColor(0); +// legend4->AddEntry(hPt, "p_{T} dist (all)","L"); +// legend4->AddEntry(hPtPionNoCorr, "p_{T} dist (Pions)","L"); +// legend4->AddEntry(hPtKaonNoCorr, "p_{T} dist (Kaons)","L"); +// legend4->AddEntry(hPtProtonNoCorr, "p_{T} dist (Protons)","L"); +// legend4->Draw(); + +// gPad->SetLogy(); + + +// c5->cd(2); +// TH1F* hPtPionCorr = static_cast(hPtPionNoCorr->Clone()); +// hPtPionCorr->SetTitle("p_{T} distribution (with corrections)"); +// hPtPionCorr->Multiply(hAccepatncePt); +// hPtPionCorr->Divide(hEfficiencyPID); +// hPtPionCorr->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPtPionCorr->DrawCopy(); + +// TH1F* hPtKaonCorr = static_cast(hPtKaonNoCorr->Clone()); +// hPtKaonCorr->Multiply(hAccepatncePt); +// hPtKaonCorr->Divide(hEfficiencyPID); +// hPtKaonCorr->Multiply(hKDecayCorr); +// hPtKaonCorr->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPtKaonCorr->DrawCopy("same"); + +// TH1F* hPtProtonCorr = static_cast(hPtProtonNoCorr->Clone()); +// hPtProtonCorr->Multiply(hAccepatncePt); +// hPtProtonCorr->Divide(hEfficiencyPID); +// hPtProtonCorr->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPtProtonCorr->DrawCopy("same"); + +// hPt->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPt->DrawCopy("same"); + +// TLegend* legend5 = new TLegend(0.63, 0.63, 0.88, 0.88); +// legend5->SetBorderSize(0); +// legend5->SetFillColor(0); +// legend5->AddEntry(hPt, "p_{T} dist (all)","L"); +// legend5->AddEntry(hPtPionCorr, "p_{T} dist (Pions)","L"); +// legend5->AddEntry(hPtKaonCorr, "p_{T} dist (Kaons)","L"); +// legend5->AddEntry(hPtProtonCorr, "p_{T} dist (Protons)","L"); +// legend5->Draw(); + +// gPad->SetLogy(); + +// c5->Update(); + +// //----------------- +// if (fAnalysisMC){ +// TCanvas* c6 = new TCanvas("c6", "c6", 100, 100, 1200, 800); +// c6->Divide(2,2); + +// c6->cd(1); +// TH1F* hPt_clone = static_cast(hPt->Clone()); +// hPt_clone->GetYaxis()->SetRangeUser(1E-5,1E1); +// hPt_clone->SetTitle("p_{T} distribution (all)"); +// hPt_clone->SetLineColor(1); +// hPt_clone->DrawCopy(); + +// TH1F* hPtMC = dynamic_cast (fListOfHists->FindObject("fdNdPtMC")); +// hPtMC->Scale(1.0/nEvents/normalization/hPtMC->GetBinWidth(1)); +// hPtMC->GetYaxis()->SetRangeUser(1E-5,1E1); +// hPtMC->DrawCopy("same"); + +// TLegend* legend6 = new TLegend(0.57, 0.57, 0.87, 0.87); +// legend6->SetBorderSize(0); +// legend6->SetFillColor(0); +// legend6->AddEntry(hPt_clone, "p_{T} dist (Rec)", "L"); +// legend6->AddEntry(hPtMC, "p_{T} dist (MC)", "L"); +// legend6->Draw(); + +// gPad->SetLogy(); + +// c6->cd(2); +// TH1F* hPtPion_clone = static_cast(hPtPionCorr->Clone()); +// hPtPion_clone->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPtPion_clone->SetTitle("p_{T} distribution (Pions)"); +// hPtPion_clone->SetLineColor(1); +// hPtPion_clone->DrawCopy(); + +// TH1F* hPtMCPion = dynamic_cast (fListOfHists->FindObject("fdNdPtMCPion")); +// hPtMCPion->Scale(1.0/nEvents/normalization/hPtMCPion->GetBinWidth(1)); +// hPtMCPion->GetYaxis()->SetRangeUser(1E-2,1E1); +// hPtMCPion->DrawCopy("same"); + +// TLegend* legend7 = new TLegend(0.57, 0.57, 0.87, 0.87); +// legend7->SetBorderSize(0); +// legend7->SetFillColor(0); +// legend7->AddEntry(hPtPion_clone, "p_{T} dist (Rec)", "L"); +// legend7->AddEntry(hPtMCPion, "p_{T} dist (MC)", "L"); +// legend7->Draw(); + +// gPad->SetLogy(); + +// c6->cd(3); +// TH1F* hPtKaon_clone = static_cast(hPtKaonCorr->Clone()); +// hPtKaon_clone->GetYaxis()->SetRangeUser(1E-2,1E0); +// hPtKaon_clone->SetLineColor(1); +// hPtKaon_clone->DrawCopy(); + +// TH1F* hPtMCKaon = dynamic_cast (fListOfHists->FindObject("fdNdPtMCKaon")); +// hPtMCKaon->Scale(1.0/nEvents/normalization/hPtMCKaon->GetBinWidth(1)); +// hPtMCKaon->GetYaxis()->SetRangeUser(1E-2,1E0); +// hPtMCKaon->DrawCopy("same"); + +// TLegend* legend8 = new TLegend(0.57, 0.57, 0.87, 0.87); +// legend8->SetBorderSize(0); +// legend8->SetFillColor(0); +// legend8->AddEntry(hPtKaon_clone, "p_{T} dist (Rec)", "L"); +// legend8->AddEntry(hPtMCKaon, "p_{T} dist (MC)", "L"); +// legend8->Draw(); + +// gPad->SetLogy(); + +// c6->cd(4); +// TH1F* hPtProton_clone = static_cast(hPtProtonCorr->Clone()); +// hPtProton_clone->GetYaxis()->SetRangeUser(1E-2,1E-1); +// hPtProton_clone->SetLineColor(1); +// hPtProton_clone->DrawCopy(); + +// TH1F* hPtMCProton = dynamic_cast (fListOfHists->FindObject("fdNdPtMCProton")); +// hPtMCProton->Scale(1.0/nEvents/normalization/hPtMCProton->GetBinWidth(1)); +// hPtMCProton->GetYaxis()->SetRangeUser(1E-2,1E-1); +// hPtMCProton->DrawCopy("same"); + +// TLegend* legend9 = new TLegend(0.2, 0.25, 0.5, 0.55); +// legend9->SetBorderSize(0); +// legend9->SetFillColor(0); +// legend9->AddEntry(hPtProton_clone, "p_{T} dist (Rec)", "L"); +// legend9->AddEntry(hPtMCProton, "p_{T} dist (MC)", "L"); +// legend9->Draw(); + +// gPad->SetLogy(); + +// c6->Update(); +// } + +} diff --git a/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h b/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h new file mode 100644 index 00000000000..e452308a919 --- /dev/null +++ b/PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h @@ -0,0 +1,97 @@ +#ifndef ALIANALYSISTASKPWG4PidDetEx_H +#define ALIANALYSISTASKPWG4PidDetEx_H + +#include +#include +#include +#include + +#include "AliAnalysisTaskSE.h" +#include "AliAODEvent.h" +#include "AliESDEvent.h" +#include "AliAnalysisFilter.h" + +class AliAnalysisTaskPWG4PidDetEx : public AliAnalysisTaskSE { + public: + enum TriggerMode {kMB1, kMB2, kSPDFASTOR}; + + AliAnalysisTaskPWG4PidDetEx(); + AliAnalysisTaskPWG4PidDetEx(const char *name); + virtual ~AliAnalysisTaskPWG4PidDetEx(); + + virtual void ConnectInputData(Option_t *); + virtual void CreateOutputObjects(); + virtual void Exec(Option_t *option); + virtual void Terminate(Option_t *); + + virtual void SetTrackFilter(AliAnalysisFilter* trackF) {fTrackFilter = trackF;} + void SetAnalysisType(const char* analysisType) {fAnalysisType = analysisType;} + void SetTriggerMode(TriggerMode triggermode) {fTriggerMode = triggermode;} + void SetMC(Bool_t analysisMC){fAnalysisMC = analysisMC;} + void SetPtCut(Double_t ptCut){fPtCut = ptCut;} + void SetEtaCut(Double_t etaCut){fEtaCut = etaCut;} + + Bool_t IsEventTriggered(AliVEvent* ev, TriggerMode trigger); + + void AnalyzeESD(AliESDEvent* esd); + void AnalyzeAOD(AliAODEvent* aod); + + private: + AliAnalysisTaskPWG4PidDetEx(const AliAnalysisTaskPWG4PidDetEx&); + AliAnalysisTaskPWG4PidDetEx& operator=(const AliAnalysisTaskPWG4PidDetEx&); + + Double_t IntegratedLength(AliVTrack* track) const; + Double_t MassSquared (AliVTrack* track) const; + Double_t KaonDecay (AliVTrack* track) const; + + static const Double_t fgkCtau; // distance for kaon decay + static const Double_t fgkPionMass; // pion mass + static const Double_t fgkKaonMass; // kaon mass + static const Double_t fgkProtonMass; // proton mass + + + AliESDEvent* fESD; //! ESD object + AliAODEvent* fAOD; //! AOD object + TList* fListOfHists; //! Output list of histograms + Double_t fEtaCut; // Eta cut used to select particles + Double_t fPtCut; // pT cut used to select particles + Int_t fXbins; // #bins for Pt histos range + Double_t fXmin; // min X value for histo range + Double_t fTOFCutP; // max X value for histo range; also the p cut used in TOF for PID + + AliAnalysisFilter* fTrackFilter; // Track Filter + Bool_t fAnalysisMC; // Flag for MC analysis + TString fAnalysisType; // "ESD" or "AOD" + TriggerMode fTriggerMode; // Trigger mode + + + //Histograms + TH1I* fEvents; //! #analyzed events + TH1F* fEffTot; //! pT for all charged particles + TH1F* fEffPID; //! pT for charged particles with TOF signal + TH1F* fAccP; //! pT for charged particles with p < fTOFCutP + TH1F* fAccPt; //! pT for charged particles with pT < fTOFCutP + TProfile* fKaonDecayCorr; //! decay correction for Kaons + TH1F* fdNdPt; //! pT dist (Rec) + TH2F* fMassAll; //! mass calculated from TOF vs p + TH1F* fdNdPtPion; //! pT for pions identified with TOF + TH2F* fMassPion; //! mass for pions identified with TOF + TH2F* fdEdxTPCPion; //! dE/dx vs p (TPC) for pions identified with TOF + TH2F* fbgTPCPion; //! dE/dx vs betagamma (TPC) for pions identified with TOF + TH1F* fdNdPtKaon; //! pT for kaons identified with TOF + TH2F* fMassKaon; //! mass for kaons identified with TOF + TH2F* fdEdxTPCKaon; //! dE/dx vs p (TPC) for kaons identified with TOF + TH2F* fbgTPCKaon; //! dE/dx vs betagamma (TPC) for kaons identified with TOF + TH1F* fdNdPtProton; //! pT for protons identified with TOF + TH2F* fMassProton; //! mass for protons identified with TOF + TH2F* fdEdxTPCProton; //! dE/dx vs p (TPC) for protons identified with TOF + TH2F* fbgTPCProton; //! dE/dx vs betagamma (TPC) for protons identified with TOF + TH1F* fdNdPtMC; //! pT dist (MC) + TH1F* fdNdPtMCPion; //! pT dist for pions (MC) + TH1F* fdNdPtMCKaon; //! pT dist for kaons (MC) + TH1F* fdNdPtMCProton; //! pT dist for protons (MC) + + ClassDef(AliAnalysisTaskPWG4PidDetEx, 1); //Analysis task for PWG4 PID using detector signals +}; + +#endif diff --git a/PWG4/JetTasks/AliAnalysisTaskUE.cxx b/PWG4/JetTasks/AliAnalysisTaskUE.cxx index e00296166aa..982c91e6f1c 100644 --- a/PWG4/JetTasks/AliAnalysisTaskUE.cxx +++ b/PWG4/JetTasks/AliAnalysisTaskUE.cxx @@ -102,7 +102,7 @@ fhMinRegAvePt(0x0), fhMinRegSumPt(0x0), fhMinRegMaxPtPart(0x0), fhMinRegSumPtvsMult(0x0), -fhdNdEta_PhiDist(0x0), +fhdNdEtaPhiDist(0x0), fhFullRegPartPtDistVsEt(0x0), fhTransRegPartPtDistVsEt(0x0), fhRegionSumPtMaxVsEt(0x0), @@ -361,7 +361,7 @@ void AliAnalysisTaskUE::AnalyseUE() Double_t deltaPhi = jetVect[0].DeltaPhi(partVect)+k270rad; if( deltaPhi > 2.*TMath::Pi() ) deltaPhi-= 2.*TMath::Pi(); - fhdNdEta_PhiDist->Fill( deltaPhi ); + fhdNdEtaPhiDist->Fill( deltaPhi ); fhFullRegPartPtDistVsEt->Fill( part->Pt(), maxPtJet1 ); Int_t region = IsTrackInsideRegion( jetVect, &partVect ); @@ -570,11 +570,11 @@ TObjArray* AliAnalysisTaskUE::FindChargedParticleJets() Double_t r = TMath::Sqrt( (jet->Eta()-track1->Eta())*(jet->Eta()-track1->Eta()) + dphi*dphi ); if( r < fConeRadius ) { - Double_t Pt = jet->E()+track1->Pt(); // Scalar sum of Pt + Double_t fPt = jet->E()+track1->Pt(); // Scalar sum of Pt // recalculating the centroid - Double_t eta = jet->Eta()*jet->E()/Pt + track1->Eta()/Pt; - Double_t phi = jet->Phi()*jet->E()/Pt + track1->Phi()/Pt; - jet->SetPtEtaPhiE( 1., eta, phi, Pt ); + Double_t eta = jet->Eta()*jet->E()/fPt + track1->Eta()/fPt; + Double_t phi = jet->Phi()*jet->E()/fPt + track1->Phi()/fPt; + jet->SetPtEtaPhiE( 1., eta, phi, fPt ); tracks.Remove( track1 ); } } @@ -717,11 +717,11 @@ void AliAnalysisTaskUE::CreateHistos() fhMinRegSumPtvsMult->Sumw2(); fListOfHistos->Add( fhMinRegSumPtvsMult ); // At(7); - fhdNdEta_PhiDist = new TH1F("hdNdEta_PhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi()); - fhdNdEta_PhiDist->SetXTitle("#Delta#phi"); - fhdNdEta_PhiDist->SetYTitle("dN_{ch}/d#etad#phi"); - fhdNdEta_PhiDist->Sumw2(); - fListOfHistos->Add( fhdNdEta_PhiDist ); // At(8) + fhdNdEtaPhiDist = new TH1F("hdNdEtaPhiDist", "Charge particle density |#eta|< 1 vs #Delta#phi", 120, 0., 2.*TMath::Pi()); + fhdNdEtaPhiDist->SetXTitle("#Delta#phi"); + fhdNdEtaPhiDist->SetYTitle("dN_{ch}/d#etad#phi"); + fhdNdEtaPhiDist->Sumw2(); + fListOfHistos->Add( fhdNdEtaPhiDist ); // At(8) // Can be use to get part pt distribution for differente Jet Pt bins fhFullRegPartPtDistVsEt = new TH2F("hFullRegPartPtDistVsEt", "dN/dP_{T} |#eta|<1 vs Leading Jet P_{T}", @@ -826,7 +826,7 @@ void AliAnalysisTaskUE::Terminate(Option_t */*option*/) return; } fhEleadingPt = (TH1F*)fListOfHistos->At(1); - fhdNdEta_PhiDist = (TH1F*)fListOfHistos->At(8); + fhdNdEtaPhiDist = (TH1F*)fListOfHistos->At(8); fhRegionSumPtMaxVsEt = (TH1F*)fListOfHistos->At(11); fhRegionSumPtMinVsEt = (TH1F*)fListOfHistos->At(12); fhRegionMultMaxVsEt = (TH1F*)fListOfHistos->At(14); @@ -895,9 +895,9 @@ void AliAnalysisTaskUE::Terminate(Option_t */*option*/) gPad->SetLogy(); c2->cd(2); - fhdNdEta_PhiDist->SetMarkerStyle(20); - fhdNdEta_PhiDist->SetMarkerColor(2); - fhdNdEta_PhiDist->DrawCopy("p"); + fhdNdEtaPhiDist->SetMarkerStyle(20); + fhdNdEtaPhiDist->SetMarkerColor(2); + fhdNdEtaPhiDist->DrawCopy("p"); gPad->SetLogy(); c2->cd(3); diff --git a/PWG4/JetTasks/AliAnalysisTaskUE.h b/PWG4/JetTasks/AliAnalysisTaskUE.h index 5a8a904611e..39a0db0c7a3 100644 --- a/PWG4/JetTasks/AliAnalysisTaskUE.h +++ b/PWG4/JetTasks/AliAnalysisTaskUE.h @@ -141,7 +141,7 @@ class AliAnalysisTaskUE : public AliAnalysisTask TH1F* fhMinRegMaxPtPart; //! TH1F* fhMinRegSumPtvsMult; //! - TH1F* fhdNdEta_PhiDist; //! + TH1F* fhdNdEtaPhiDist; //! TH2F* fhFullRegPartPtDistVsEt; //! TH2F* fhTransRegPartPtDistVsEt; //! diff --git a/PWG4/PWG4JetTasksLinkDef.h b/PWG4/PWG4JetTasksLinkDef.h index 4ef44c8e73d..be8fc196f58 100644 --- a/PWG4/PWG4JetTasksLinkDef.h +++ b/PWG4/PWG4JetTasksLinkDef.h @@ -8,7 +8,7 @@ #pragma link C++ class AliAnalysisTaskJetSpectrum+; #pragma link C++ class AliAnalysisHelperJetTasks+; #pragma link C++ class AliAnaESDSpectraQA+; +#pragma link C++ class AliAnalysisTaskPWG4PidDetEx+; #pragma link C++ class AliJetSpectrumUnfolding+; - #endif diff --git a/PWG4/libPWG4JetTasks.pkg b/PWG4/libPWG4JetTasks.pkg index ebf8ecb218c..a68e5154308 100644 --- a/PWG4/libPWG4JetTasks.pkg +++ b/PWG4/libPWG4JetTasks.pkg @@ -1,6 +1,6 @@ #-*- Mode: Makefile -*- -SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx JetTasks/AliJetSpectrumUnfolding.cxx +SRCS = JetTasks/AliAnalysisTaskUE.cxx JetTasks/AliAnalysisTaskJetSpectrum.cxx JetTasks/AliAnalysisHelperJetTasks.cxx JetTasks/AliAnaESDSpectraQA.cxx JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx JetTasks/AliJetSpectrumUnfolding.cxx HDRS:= $(SRCS:.cxx=.h) -- 2.39.3