Added PID task, some fixes for coding conventions
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 15 Apr 2009 14:23:29 +0000 (14:23 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Wed, 15 Apr 2009 14:23:29 +0000 (14:23 +0000)
PWG4/CMake_libPWG4JetTasks.txt
PWG4/JetTasks/AliAnaESDSpectraQA.cxx
PWG4/JetTasks/AliAnalysisTaskJetSpectrum.cxx
PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskPWG4PidDetEx.h [new file with mode: 0644]
PWG4/JetTasks/AliAnalysisTaskUE.cxx
PWG4/JetTasks/AliAnalysisTaskUE.h
PWG4/PWG4JetTasksLinkDef.h
PWG4/libPWG4JetTasks.pkg

index 2ad98b5..59b1bc0 100644 (file)
@@ -5,6 +5,7 @@ set(SRCS
        JetTasks/AliAnalysisTaskJetSpectrum.cxx
        JetTasks/AliAnalysisHelperJetTasks.cxx
        JetTasks/AliAnaESDSpectraQA.cxx
+       JetTasks/AliAnalysisTaskPWG4PidDetEx.cxx
        JetTasks/AliJetSpectrumUnfolding.cxx
 )
 
index 303738c..ffafff5 100644 (file)
@@ -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);
index a9da268..5a2947d 100644 (file)
@@ -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 (file)
index 0000000..f43fe91
--- /dev/null
@@ -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 <TROOT.h>
+#include <TSystem.h>
+#include <TChain.h>
+#include <TTree.h>
+#include <TFile.h>
+#include <TList.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TMath.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TF1.h>
+#include <TProfile.h>
+#include <TLegend.h>
+#include <TDatabasePDG.h>
+#include <TParticle.h>
+
+#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 <iostream>
+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<TTree*> (GetInputData(0));
+  if (!tree) {
+    Printf("ERROR: Could not read chain from input slot 0");
+  } 
+  else {
+    if(fAnalysisType == "ESD") {
+      AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());    
+      if (!esdH) {
+       Printf("ERROR: Could not get ESDInputHandler");
+      } else
+       fESD = esdH->GetEvent();
+    }
+    else if(fAnalysisType == "AOD") {
+      AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (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<AliMCEventHandler*> (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<TList*> (GetOutputData(0));
+//   if (!fListOfHists) {
+//     Printf("ERROR: fListOfHists not available");
+//     return;
+//   }
+//   TH1I* hevents = dynamic_cast<TH1I*> (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<TH2F*> (fListOfHists->FindObject("fMassAll"));
+//   hMassAll->SetLineColor(1);
+//   hMassAll->SetMarkerColor(1);
+//   hMassAll->DrawCopy();
+
+//   TH2F* hMassPion =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassPion"));
+//   hMassPion->SetLineColor(2);
+//   hMassPion->SetMarkerColor(2);
+//   hMassPion->SetMarkerStyle(7);
+//   hMassPion->DrawCopy("same");
+
+//   TH2F* hMassKaon =  dynamic_cast<TH2F*> (fListOfHists->FindObject("fMassKaon"));
+//   hMassKaon->SetLineColor(3);
+//   hMassKaon->SetMarkerColor(3);
+//   hMassKaon->SetMarkerStyle(7);
+//   hMassKaon->DrawCopy("same");
+
+//   TH2F* hMassProton =  dynamic_cast<TH2F*> (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<TH2F*> (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<TH2F*> (fListOfHists->FindObject("fdEdxTPCKaon"));
+//   hdEdxTPCKaon->SetLineColor(3);
+//   hdEdxTPCKaon->SetMarkerColor(3);
+//   hdEdxTPCKaon->SetMarkerStyle(7);
+//   hdEdxTPCKaon->DrawCopy("same");
+
+//   TH2F* hdEdxTPCProton =  dynamic_cast<TH2F*> (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<TH2F*> (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<TH2F*> (fListOfHists->FindObject("fbgTPCKaon"));
+//   hdEdxTPCbgKaon->SetLineColor(3);
+//   hdEdxTPCbgKaon->SetMarkerColor(3);
+//   hdEdxTPCbgKaon->SetMarkerStyle(7);
+//   hdEdxTPCbgKaon->DrawCopy("same");
+
+//   TH2F* hdEdxTPCbgProton =  dynamic_cast<TH2F*> (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<TH1F*> (fListOfHists->FindObject("fAccPt"));
+//   TH1F* hAcceptanceP = dynamic_cast<TH1F*> (fListOfHists->FindObject("fAccP"));
+//   hAccepatncePt->Divide(hAcceptanceP);
+//   hAccepatncePt->SetTitle("Acceptance correction");
+//   hAccepatncePt->GetYaxis()->SetTitle("Acceptance correction");
+//   hAccepatncePt->DrawCopy();
+
+//   c4->cd(2);
+//   TH1F* hEfficiencyPID = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffPID"));
+//   TH1F* hEfficiencyTot = dynamic_cast<TH1F*> (fListOfHists->FindObject("fEffTot"));
+//   hEfficiencyPID->Divide(hEfficiencyTot);
+//   hEfficiencyPID->SetTitle("Efficiency correction");
+//   hEfficiencyPID->GetYaxis()->SetTitle("Efficiency correction");
+//   hEfficiencyPID->DrawCopy();
+
+//   c4->cd(3);
+//   TProfile* hKDecayCorr = dynamic_cast<TProfile*> (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<TH1F*> (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<TH1F*> (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<TH1F*> (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<TH1F*> (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<TH1F*>(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<TH1F*>(hPtKaonNoCorr->Clone());
+//   hPtKaonCorr->Multiply(hAccepatncePt);
+//   hPtKaonCorr->Divide(hEfficiencyPID);
+//   hPtKaonCorr->Multiply(hKDecayCorr);
+//   hPtKaonCorr->GetYaxis()->SetRangeUser(1E-2,1E1);
+//   hPtKaonCorr->DrawCopy("same");
+
+//   TH1F* hPtProtonCorr = static_cast<TH1F*>(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<TH1F*>(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<TH1F*> (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<TH1F*>(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<TH1F*> (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<TH1F*>(hPtKaonCorr->Clone()); 
+//     hPtKaon_clone->GetYaxis()->SetRangeUser(1E-2,1E0);
+//     hPtKaon_clone->SetLineColor(1);
+//     hPtKaon_clone->DrawCopy();
+
+//     TH1F* hPtMCKaon = dynamic_cast<TH1F*> (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<TH1F*>(hPtProtonCorr->Clone()); 
+//     hPtProton_clone->GetYaxis()->SetRangeUser(1E-2,1E-1);
+//     hPtProton_clone->SetLineColor(1);
+//     hPtProton_clone->DrawCopy();
+
+//     TH1F* hPtMCProton = dynamic_cast<TH1F*> (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 (file)
index 0000000..e452308
--- /dev/null
@@ -0,0 +1,97 @@
+#ifndef ALIANALYSISTASKPWG4PidDetEx_H
+#define ALIANALYSISTASKPWG4PidDetEx_H
+
+#include <TList.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TProfile.h>
+
+#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
index e002961..982c91e 100644 (file)
@@ -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);      
index 5a8a904..39a0db0 100644 (file)
@@ -141,7 +141,7 @@ class  AliAnalysisTaskUE : public AliAnalysisTask
     TH1F*  fhMinRegMaxPtPart;        //!
     TH1F*  fhMinRegSumPtvsMult;      //!
     
-    TH1F*  fhdNdEta_PhiDist;         //!
+    TH1F*  fhdNdEtaPhiDist;         //!
     TH2F*  fhFullRegPartPtDistVsEt;  //!
     TH2F*  fhTransRegPartPtDistVsEt; //!
     
index 4ef44c8..be8fc19 100644 (file)
@@ -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
index ebf8ecb..a68e515 100644 (file)
@@ -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)