From: pchrist <pchrist@f7af4fe6-9843-0410-8265-dc069ae4e863>
Date: Mon, 8 Feb 2010 14:32:51 +0000 (+0000)
Subject: Implementing better the kSigma2 approach - dEdx(measured)/dEdx(theory)
X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=f81f1b195228c75423c4f838f2be6de0f27b9004

Implementing better the kSigma2 approach - dEdx(measured)/dEdx(theory)
---

diff --git a/PWG2/SPECTRA/AliProtonAnalysisBase.cxx b/PWG2/SPECTRA/AliProtonAnalysisBase.cxx
index 3d883289ec0..9c9982977f3 100644
--- a/PWG2/SPECTRA/AliProtonAnalysisBase.cxx
+++ b/PWG2/SPECTRA/AliProtonAnalysisBase.cxx
@@ -30,6 +30,7 @@
 #include <AliPID.h>
 #include <AliVertexerTracks.h>
 #include <AliESDpid.h>
+#include <AliTPCPIDResponse.h>
 class AliLog;
 class AliESDVertex;
 
@@ -917,13 +918,28 @@ Bool_t AliProtonAnalysisBase::IsProton(AliESDtrack *track) {
       gP = tpcTrack->P();
       gEta = tpcTrack->Eta();
     }
-    //We start the P slices at >0.3GeV/c with a bining of 50MeV/c ==> Int_t(0.3001/0.05) = 6
-    Int_t nbinP = Int_t(gP/0.05) - 6;
-    Double_t tpcSignal = track->GetTPCsignal();
-    Double_t dEdxTheory = Bethe(gP/9.38270000000000048e-01);
-    Double_t dEdxSigma = fdEdxSigma[nbinP];
-    Double_t nsigma = TMath::Abs(tpcSignal - dEdxTheory)/(tpcSignal*(dEdxSigma/TMath::Sqrt(track->GetTPCsignalN())));
-    if(nsigma <= fNSigma) 
+    Double_t fAlephParameters[5];
+    if(fAnalysisMC) {
+      fAlephParameters[0] = 2.15898e+00/50.;
+      fAlephParameters[1] = 1.75295e+01;
+      fAlephParameters[2] = 3.40030e-09;
+      fAlephParameters[3] = 1.96178e+00;
+      fAlephParameters[4] = 3.91720e+00;
+    }
+    else {
+      fAlephParameters[0] = 0.0283086;
+      fAlephParameters[1] = 2.63394e+01;
+      fAlephParameters[2] = 5.04114e-11;
+      fAlephParameters[3] = 2.12543e+00;
+      fAlephParameters[4] = 4.88663e+00;
+    }
+
+    AliESDpid *fESDpid = new AliESDpid(); 
+    AliTPCPIDResponse tpcResponse = fESDpid->GetTPCResponse(); 
+    tpcResponse.SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
+    Double_t normalizeddEdx = TMath::Log(track->GetTPCsignal()/tpcResponse.GetExpectedSignal(gP,AliPID::kProton));
+
+    if(normalizeddEdx >= -0.15)
       return kTRUE;
   }
 
diff --git a/PWG2/SPECTRA/macros/drawProtonResults.C b/PWG2/SPECTRA/macros/drawProtonResults.C
index 652ff794cd6..ffffd295b88 100644
--- a/PWG2/SPECTRA/macros/drawProtonResults.C
+++ b/PWG2/SPECTRA/macros/drawProtonResults.C
@@ -139,7 +139,7 @@ void drawQAPlots(const char* analysisOutput,
   //2D de/dx vs P
   TH2F *gHistdEdxP = dynamic_cast<TH2F *>(fQA2DList->At(0));
   gHistdEdxP->SetStats(kFALSE);
-  drawdEdx(gHistdEdxP);
+  drawdEdx(gHistdEdxP,0);
   TH2F *gHistProtonsdEdxP = dynamic_cast<TH2F *>(fQA2DList->At(1));
   gHistProtonsdEdxP->SetStats(kFALSE);
 
@@ -201,7 +201,7 @@ void drawQAPlots(const char* analysisOutput,
   //2D de/dx vs P
   TH2F *gHistZP = dynamic_cast<TH2F *>(fQA2DList->At(2));
   gHistZP->SetStats(kFALSE);
-  //drawdEdx(gHistdEdxP);
+  drawdEdx(gHistZP,1);
   TH2F *gHistProtonsZP = dynamic_cast<TH2F *>(fQA2DList->At(3));
   gHistProtonsZP->SetStats(kFALSE);
 
@@ -314,6 +314,11 @@ void drawQAPlots(const char* analysisOutput,
   
   cdEdx->cd(2)->SetLogx(); gHistProtonsdEdxP->Draw("col");
 
+  TCanvas *cZdEdx = new TCanvas("cZdEdx","Normalized dE/dx (TPC)",500,0,700,400);
+  cZdEdx->SetFillColor(10); cZdEdx->SetHighLightColor(10); cZdEdx->Divide(2,1);
+  cZdEdx->cd(1); gHistZP->Draw("col");
+  cZdEdx->cd(2); gHistProtonsZP->Draw("col");
+
   TCanvas *cEtaPhi = new TCanvas("cEtaPhi",
 				 "eta-phi",
 				 0,0,700,400);
@@ -703,8 +708,10 @@ void PrintYields(TH1 *h) {
 }
 
 //___________________________________________________//
-void drawdEdx(TH2F *gHistdEdxP) {
+void drawdEdx(TH2F *gHistdEdxP, Int_t iMode) {
   //Draws the dE/dx distributions for the different momentum bins
+  //iMode == 0: dEdx vs P
+  //iMode == 1: normalized(dEdx) vs P
   TString title;
   TH1D *gHist[100];
   Int_t iCounter = 0;
@@ -717,16 +724,28 @@ void drawdEdx(TH2F *gHistdEdxP) {
 
   for(Int_t iBin = 1; iBin <= gHistdEdxP->GetNbinsX(); iBin++) {
     if((binMax > 0.41)&&(binMin < 1.01)) {
-      title = "P: "; title += binMin; title += " - "; 
-      title += binMax; title += "GeV/c";
-      canvasTitle = "dedxMomentumSlice."; canvasTitle += iCounter;
-      canvasTitle += ".gif";
-      c[iCounter] = new TCanvas(title.Data(),title.Data(),0,0,500,500);
+      if(iMode == 0) {
+	title = "(dEdx)P: "; title += binMin; title += " - "; 
+	title += binMax; title += "GeV/c";
+	canvasTitle = "dedxMomentumSlice."; canvasTitle += iCounter;
+	canvasTitle += ".gif";
+	c[iCounter] = new TCanvas(title.Data(),title.Data(),0,0,500,500);
+      }
+      if(iMode == 1) {
+	title = "(normdEdx)P: "; title += binMin; title += " - "; 
+	title += binMax; title += "GeV/c";
+	canvasTitle = "normdedxMomentumSlice."; canvasTitle += iCounter;
+	canvasTitle += ".gif";
+	c[iCounter] = new TCanvas(title.Data(),title.Data(),500,0,500,500);
+      }
       c[iCounter]->SetFillColor(10); c[iCounter]->SetHighLightColor(10); 
       gHist[iCounter] = gHistdEdxP->ProjectionY(title.Data(),iBin,iBin+1);
       gHist[iCounter]->SetTitle(title.Data());
       gHist[iCounter]->SetStats(kFALSE);
-      gHist[iCounter]->GetXaxis()->SetRangeUser(0.0,300.);
+      if(iMode == 0)
+	gHist[iCounter]->GetXaxis()->SetRangeUser(0.0,300.);
+      if(iMode == 1)
+	gHist[iCounter]->GetXaxis()->SetRangeUser(-2.0,2.0);
       if(gHist[iCounter]->GetEntries() != 0)
 	c[iCounter]->SetLogy();
       gHist[iCounter]->Draw();