ITS/ParametrizeITSBetheBloch.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include <TCanvas.h>
   3 #include <TGraph.h>
   4 #include <TGraphErrors.h>
   5 #include <TGraphAsymmErrors.h>
   6 #include <TH1.h>
   7 #include <TF1.h>
   8 #include <TH2.h>
   9 #include <TFile.h>
  10 #include <TList.h>
  11 #include <TString.h>
  12 #include <TDirectoryFile.h>
  13 #include <TLatex.h>
  14 #include <TStyle.h>
  15 #include <TMath.h>
  16 #include <TDatabasePDG.h>
  17 #include "AliITSPIDResponse.h"
  18 #endif
  19
  20 /*  $Id$    */
  21
  22
  23 // Macro to extract the parameterizations of dE/dx in ITS vs. p
  24 // using Phobos BB parameterization
  25
  26
  27 Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par);
  28 Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par);
  29 Double_t BetheBlochAleph(const Double_t *x, const Double_t *par);
  30
  31
  32 void ParametrizeITSBetheBloch(Bool_t optFromITSsa=kFALSE){
  33
  34   TFile *fil=new TFile("AnalysisResults.root");
  35   TH2F *hsigvsp=0x0;
  36   TString dfName="PWG2SpectraITSsa";
  37   TString listName="clistITSsaMult-1to-1";
  38   TString hisName="fHistDEDX";
  39
  40   if(optFromITSsa){
  41     dfName="ITSsaTracks";
  42     listName="clistITSsaTracks";
  43     hisName="hdedxvsPt4clsITSpureSA";
  44   }
  45   TDirectoryFile* df=(TDirectoryFile*)fil->Get(dfName.Data());
  46   if(!df){
  47     printf("ERROR: directory file %s not found\n",dfName.Data());
  48     return;
  49   }
  50   TList* l=(TList*)df->Get(listName.Data());
  51   if(!l){
  52     printf("ERROR: TList %s not found\n",listName.Data());
  53     return;
  54   }
  55   hsigvsp=(TH2F*)l->FindObject(hisName.Data());
  56   if(!hsigvsp){
  57     printf("ERROR: historgam %s not found\n",hisName.Data());
  58     return;
  59   }
  60   hsigvsp->GetXaxis()->SetTitle("momentum (GeV/c)");
  61   hsigvsp->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
  62
  63   Double_t minpPionFit=0.1;
  64   Double_t maxpPionFit=0.8;
  65   Double_t minbgGlobFit=0.3;
  66   Double_t maxbgGlobFit=2.5;
  67   Double_t minbgPions=2.;
  68   Double_t maxbgPions=10.;
  69   Double_t minpKaons=0.2;
  70   Double_t maxpKaons=0.4;
  71   Double_t maxbgKaons=0.7;
  72   Double_t minpProtons=0.35;
  73   Double_t maxpProtons=0.57;
  74   Double_t maxbgProtons=0.7;
  75
  76   Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
  77   Double_t massKaon=TDatabasePDG::Instance()->GetParticle(321)->Mass();
  78   Double_t massProton=TDatabasePDG::Instance()->GetParticle(2212)->Mass();
  79   Double_t massDeuteron=1.8756;
  80
  81   Int_t iPad=1;
  82   Int_t iCanvas=0;
  83   Int_t nbinsX=hsigvsp->GetNbinsX();
  84   printf("Total Number of bins =%d\n",nbinsX);
  85   Int_t group=2;
  86
  87   Int_t totPads=(Int_t)((Float_t)nbinsX/(Float_t)group+0.5);
  88   Int_t totCanvas=totPads/16;
  89   TCanvas** c0=new TCanvas*[totCanvas];
  90   c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
  91   c0[iCanvas]->Divide(4,4);
  92   TCanvas** c0k=new TCanvas*[totCanvas];
  93   c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
  94   c0k[iCanvas]->Divide(4,4);
  95   TCanvas** c0p=new TCanvas*[totCanvas];
  96   c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
  97   c0p[iCanvas]->Divide(4,4);
  98
  99   TGraphErrors* gdedxpi=new TGraphErrors(0);
 100   TGraphErrors* gdedxk=new TGraphErrors(0);
 101   TGraphErrors* gdedxp=new TGraphErrors(0);
 102   TGraphErrors* gdedxbg=new TGraphErrors(0);
 103   gdedxbg->SetTitle("");
 104
 105   Int_t nPpi=0;
 106   Int_t nPk=0;
 107   Int_t nPp=0;
 108   Int_t nPbg=0;
 109   for(Int_t iBin=1; iBin<=nbinsX; iBin+=group){
 110     if(iPad>16){
 111       iCanvas++;
 112       c0[iCanvas]=new TCanvas(Form("c0_%d",iCanvas),Form("c0_%d",iCanvas));
 113       c0[iCanvas]->Divide(4,4);
 114       c0k[iCanvas]=new TCanvas(Form("c0k_%d",iCanvas),Form("c0k_%d",iCanvas));
 115       c0k[iCanvas]->Divide(4,4);
 116       c0p[iCanvas]=new TCanvas(Form("c0p_%d",iCanvas),Form("c0p_%d",iCanvas));
 117       c0p[iCanvas]->Divide(4,4);
 118       iPad=1;
 119     }
 120
 121     Float_t sum=0;
 122     Float_t sumwei=0;
 123     Double_t pmed=0.;
 124     for(Int_t ibx=iBin; ibx<iBin+group;ibx++){
 125       for(Int_t iby=1; iby<=hsigvsp->GetNbinsY();iby++){
 126         sum+=hsigvsp->GetXaxis()->GetBinCenter(ibx)*hsigvsp->GetBinContent(ibx,iby);
 127         sumwei+=hsigvsp->GetBinContent(ibx,iby);
 128       }
 129     }
 130     if(sumwei>0.) pmed=sum/sumwei;
 131     if(pmed<0.1) continue;
 132     Int_t iPmed=(Int_t)(pmed*1000);
 133     TH1D* hProj=hsigvsp->TH2F::ProjectionY(Form("proj%d",iPmed),iBin,iBin+group-1);
 134     TH1D* hKaon=(TH1D*)hProj->Clone(Form("projk%d",iPmed));
 135     c0[iCanvas]->cd(iPad);
 136     hProj->SetLineColor(2);
 137     hProj->Draw();
 138
 139     Float_t peakpos=hProj->GetBinCenter(hProj->GetMaximumBin());
 140     Float_t peaksig=0.1*peakpos;
 141     hProj->Fit("gaus","LL","",peakpos-3*peaksig,peakpos+3.*peaksig);
 142     TF1* fgaus=(TF1*)hProj->GetListOfFunctions()->FindObject("gaus");
 143     gdedxpi->SetPoint(nPpi,sum/sumwei,fgaus->GetParameter(1));
 144     gdedxpi->SetPointError(nPpi,0.,fgaus->GetParError(1));
 145     nPpi++;
 146     if(pmed/massPion>minbgPions && pmed/massPion<maxbgPions){
 147       gdedxbg->SetPoint(nPbg,pmed/massPion,fgaus->GetParameter(1));
 148       gdedxbg->SetPointError(nPbg,0.,fgaus->GetParError(1));
 149       nPbg++;
 150     }
 151     for(Int_t iBinK=1; iBinK<=hKaon->GetNbinsX(); iBinK++){
 152       Double_t cont=hKaon->GetBinContent(iBinK)-fgaus->Eval(hKaon->GetBinCenter(iBinK));
 153       if(cont<0) cont=0.;
 154       if(hKaon->GetBinCenter(iBinK)<fgaus->GetParameter(1)) cont=0;
 155       hKaon->SetBinContent(iBinK,cont);
 156     }
 157     c0k[iCanvas]->cd(iPad);
 158     hKaon->SetLineColor(4);
 159     hKaon->Draw();
 160     TH1D* hProton=(TH1D*)hKaon->Clone(Form("projp%d",iPmed));
 161     TF1* fgausk=0x0;
 162     Double_t mink=250-600*(pmed-0.2);
 163     Double_t maxk=mink+mink;
 164     hKaon->Fit("gaus","LL","",mink,maxk);
 165     fgausk=(TF1*)hKaon->GetListOfFunctions()->FindObject("gaus");
 166     if(pmed>minpKaons && pmed<maxpKaons){
 167       gdedxk->SetPoint(nPk,pmed,fgausk->GetParameter(1));
 168       gdedxk->SetPointError(nPk,0.,fgausk->GetParError(1));
 169       nPk++;
 170       if(pmed/massKaon<maxbgKaons){
 171         gdedxbg->SetPoint(nPbg,pmed/massKaon,fgausk->GetParameter(1));
 172         gdedxbg->SetPointError(nPbg,0.,fgausk->GetParError(1));
 173         nPbg++;
 174       }
 175     }
 176     if(fgausk){
 177       for(Int_t iBinP=1; iBinP<=hProton->GetNbinsX(); iBinP++){
 178         Double_t cont=hKaon->GetBinContent(iBinP)-fgausk->Eval(hKaon->GetBinCenter(iBinP));
 179         if(cont<0) cont=0.;
 180         if(hProton->GetBinCenter(iBinP)<fgausk->GetParameter(1)) cont=0;
 181         hProton->SetBinContent(iBinP,cont);
 182       }
 183     }
 184     c0p[iCanvas]->cd(iPad);
 185     hProton->SetLineColor(kGreen+1);
 186     hProton->Draw();
 187     if(pmed>minpProtons && pmed<maxpProtons){
 188       Double_t minP=300-600*(pmed-0.35);
 189       Double_t maxP=minP+minP;
 190       hProton->Fit("gaus","LL","",minP,maxP);
 191       TF1* fgausp=(TF1*)hProton->GetListOfFunctions()->FindObject("gaus");
 192       gdedxp->SetPoint(nPp,pmed,fgausp->GetParameter(1));
 193       gdedxp->SetPointError(nPp,0.,fgausp->GetParError(1));
 194       nPp++;
 195       if(pmed/massProton<maxbgProtons){
 196         gdedxbg->SetPoint(nPbg,pmed/massProton,fgausp->GetParameter(1));
 197         gdedxbg->SetPointError(nPbg,0.,fgausp->GetParError(1));
 198         nPbg++;
 199       }
 200     }
 201     ++iPad;
 202   }
 203
 204   gStyle->SetPalette(1);
 205   hsigvsp->GetXaxis()->SetRangeUser(0.07,1.8);
 206   hsigvsp->SetStats(0);
 207   TCanvas* c1=new TCanvas("c1","Pion Fit");
 208   c1->SetLogz();
 209   hsigvsp->Draw("col");
 210   gdedxpi->Draw("Psame");
 211   gdedxk->Draw("Psame");
 212   gdedxp->Draw("Psame");
 213
 214   TF1 *fPhobos=new TF1("fPhobos",BetheBlochPhobos,minpPionFit,maxpPionFit,5);
 215   fPhobos->SetParLimits(0,0.,1.E9);
 216   fPhobos->SetParLimits(1,0.,100.);
 217   fPhobos->SetParLimits(2,0.,1.);
 218   fPhobos->SetParLimits(3,0.,1.);
 219   fPhobos->SetParLimits(4,0.,1.);
 220   gdedxpi->Fit("fPhobos","R");
 221
 222   TCanvas* c1g=new TCanvas("c1g","Global Fit");
 223   TF1 *fPhobosbg=new TF1("fPhobosbg",BetheBlochPhobosBG,minbgGlobFit,maxbgGlobFit,5);
 224   fPhobosbg->SetParLimits(0,0.,1.E9);
 225   fPhobosbg->SetParLimits(1,0.,100.);
 226   fPhobosbg->SetParLimits(2,0.,1.);
 227   fPhobosbg->SetParLimits(3,0.,1.);
 228   fPhobosbg->SetParLimits(4,0.,1.);
 229   gdedxbg->Draw("AP");
 230   gdedxbg->GetXaxis()->SetTitle("#beta#gamma");
 231   gdedxbg->GetYaxis()->SetTitle("dE/dx (keV/300 #mum)");
 232   gdedxbg->Fit("fPhobosbg","R");
 233   c1g->Update();
 234
 235
 236   AliITSPIDResponse* itsresGlobFit=new AliITSPIDResponse(kFALSE);
 237   AliITSPIDResponse* itsres=new AliITSPIDResponse(kFALSE);
 238   Double_t parameters[5],parametersGlob[5];
 239
 240   printf("-------- Pion Fit Parameters -----------\n");
 241   for(Int_t iPar=0; iPar<5; iPar++){
 242     parameters[iPar]=fPhobos->GetParameter(iPar);
 243     printf("parameters[%d]=%g;\n",iPar,parameters[iPar]);
 244   }
 245   printf("-------- Global Fit Parameters -----------\n");
 246   for(Int_t iPar=0; iPar<5; iPar++){
 247     parametersGlob[iPar]=fPhobosbg->GetParameter(iPar);
 248     printf("parametersGlob[%d]=%g;\n",iPar,parametersGlob[iPar]);
 249   }
 250
 251   itsresGlobFit->SetBetheBlochParamsITSsa(parametersGlob);
 252   itsres->SetBetheBlochParamsITSsa(parameters);
 253
 254   TGraph* gPion=new TGraph(0);
 255   TGraph* gKaon=new TGraph(0);
 256   TGraph* gProton=new TGraph(0);
 257   TGraph* gPionGlobFit=new TGraph(0);
 258   TGraph* gKaonGlobFit=new TGraph(0);
 259   TGraph* gProtonGlobFit=new TGraph(0);
 260
 261   TGraphAsymmErrors* gBBpions=new TGraphAsymmErrors(0);
 262   TGraphAsymmErrors* gBBkaons=new TGraphAsymmErrors(0);
 263   TGraphAsymmErrors* gBBprotons=new TGraphAsymmErrors(0);
 264   TGraphAsymmErrors* gBBpionsGlobFit=new TGraphAsymmErrors(0);
 265   TGraphAsymmErrors* gBBkaonsGlobFit=new TGraphAsymmErrors(0);
 266   TGraphAsymmErrors* gBBprotonsGlobFit=new TGraphAsymmErrors(0);
 267
 268   Double_t numberofsigmapions=2.;
 269   Float_t resodedx[4];
 270   resodedx[0]=0.0001; //default value
 271   resodedx[1]=0.0001; //default value
 272   resodedx[2]=0.116;
 273   resodedx[3]=0.103;
 274
 275
 276   for(Int_t ip=0; ip<400; ip++){
 277     Double_t mom=0.1+0.02*ip;
 278     Double_t dedxpions=itsres->Bethe(mom,massPion,kTRUE);
 279     Double_t dedxkaons=itsres->Bethe(mom,massKaon,kTRUE);
 280     Double_t dedxprotons=itsres->Bethe(mom,massProton,kTRUE);
 281     Double_t dedxdeutons=itsres->Bethe(mom,massDeuteron,kTRUE);
 282     Double_t dedxpionsGlobFit=itsresGlobFit->Bethe(mom,massPion,kTRUE);
 283     Double_t dedxkaonsGlobFit=itsresGlobFit->Bethe(mom,massKaon,kTRUE);
 284     Double_t dedxprotonsGlobFit=itsresGlobFit->Bethe(mom,massProton,kTRUE);
 285     Double_t dedxdeutonsGlobFit=itsresGlobFit->Bethe(mom,massDeuteron,kTRUE);
 286
 287     gPion->SetPoint(ip,mom,dedxpions);
 288     gKaon->SetPoint(ip,mom,dedxkaons);
 289     gProton->SetPoint(ip,mom,dedxprotons);
 290     gPionGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
 291     gKaonGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
 292     gProtonGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
 293
 294     gBBpions->SetPoint(ip,mom,dedxpions);
 295     gBBpions->SetPointError(ip,0,0,(dedxpions*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaons)/2));
 296     gBBkaons->SetPoint(ip,mom,dedxkaons);
 297     gBBkaons->SetPointError(ip,0,0,TMath::Abs((dedxpions-dedxkaons)/2),TMath::Abs((dedxprotons-dedxkaons)/2));
 298     gBBprotons->SetPoint(ip,mom,dedxprotons);
 299     gBBprotons->SetPointError(ip,0,0,TMath::Abs((dedxprotons-dedxkaons)/2),TMath::Abs((dedxprotons-dedxdeutons)/2));
 300
 301     gBBpionsGlobFit->SetPoint(ip,mom,dedxpionsGlobFit);
 302     gBBpionsGlobFit->SetPointError(ip,0,0,(dedxpionsGlobFit*resodedx[3]*numberofsigmapions),TMath::Abs((dedxpions-dedxkaonsGlobFit)/2));
 303     gBBkaonsGlobFit->SetPoint(ip,mom,dedxkaonsGlobFit);
 304     gBBkaonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxpionsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2));
 305     gBBprotonsGlobFit->SetPoint(ip,mom,dedxprotonsGlobFit);
 306     gBBprotonsGlobFit->SetPointError(ip,0,0,TMath::Abs((dedxprotonsGlobFit-dedxkaonsGlobFit)/2),TMath::Abs((dedxprotonsGlobFit-dedxdeutonsGlobFit)/2));
 307   }
 308
 309
 310   TCanvas* c2=new TCanvas("c2","Parameterizations");
 311   c2->SetLogz();
 312   hsigvsp->SetMinimum(20);
 313   hsigvsp->Draw("col");
 314   gPion->Draw("LSAME");
 315   gKaon->Draw("LSAME");
 316   gProton->Draw("LSAME");
 317   gPionGlobFit->SetLineColor(kRed+1);
 318   gKaonGlobFit->SetLineColor(kRed+1);
 319   gProtonGlobFit->SetLineColor(kRed+1);
 320   gPionGlobFit->Draw("LSAME");
 321   gKaonGlobFit->Draw("LSAME");
 322   gProtonGlobFit->Draw("LSAME");
 323   TLatex* t1=new TLatex(0.7,0.85,"Pion Fit");
 324   t1->SetNDC();
 325   t1->Draw();
 326   TLatex* tG=new TLatex(0.7,0.75,"Global Fit");
 327   tG->SetNDC();
 328   tG->SetTextColor(kRed+1);
 329   tG->Draw();
 330
 331
 332   gBBpions->SetLineColor(2);
 333   gBBpions->SetLineWidth(2);
 334   gBBpions->SetFillColor(2);
 335   gBBpions->SetFillStyle(3001);
 336   gBBkaons->SetLineColor(3);
 337   gBBkaons->SetLineWidth(2);
 338   gBBkaons->SetFillColor(3);
 339   gBBkaons->SetFillStyle(3001);
 340   gBBprotons->SetLineColor(4);
 341   gBBprotons->SetLineWidth(2);
 342   gBBprotons->SetFillColor(4);
 343   gBBprotons->SetFillStyle(3001);
 344
 345   TCanvas* c3=new TCanvas("c3","Asymm Bands, PionFit");
 346   c3->SetLogz();
 347   c3->SetLogx();
 348   hsigvsp->SetMinimum(20);
 349   hsigvsp->Draw("col");
 350   gBBpions->Draw("L4same");
 351   gBBkaons->Draw("L4same");
 352   gBBprotons->Draw("L4same");
 353
 354   gBBpionsGlobFit->SetLineColor(2);
 355   gBBpionsGlobFit->SetLineWidth(2);
 356   gBBpionsGlobFit->SetFillColor(2);
 357   gBBpionsGlobFit->SetFillStyle(3001);
 358   gBBkaonsGlobFit->SetLineColor(3);
 359   gBBkaonsGlobFit->SetLineWidth(2);
 360   gBBkaonsGlobFit->SetFillColor(3);
 361   gBBkaonsGlobFit->SetFillStyle(3001);
 362   gBBprotonsGlobFit->SetLineColor(4);
 363   gBBprotonsGlobFit->SetLineWidth(2);
 364   gBBprotonsGlobFit->SetFillColor(4);
 365   gBBprotonsGlobFit->SetFillStyle(3001);
 366
 367   TCanvas* c3g=new TCanvas("c3g","Asymm Bands, GlobFit");
 368   c3g->SetLogz();
 369   c3g->SetLogx();
 370   hsigvsp->SetMinimum(20);
 371   hsigvsp->Draw("col");
 372   gBBpionsGlobFit->Draw("L4same");
 373   gBBkaonsGlobFit->Draw("L4same");
 374   gBBprotonsGlobFit->Draw("L4same");
 375
 376 }
 377
 378
 379
 380
 381
 382
 383 Double_t BetheBlochPhobos(const Double_t *x, const Double_t *par){
 384
 385   //
 386   // returns AliExternalTrackParam::BetheBloch normalized to
 387   // fgMIP at the minimum
 388   //
 389
 390   Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 391   Double_t bg=x[0]/massPion;
 392   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 393   Double_t gamma=bg/beta;
 394   Double_t eff=1.0;
 395   if(bg<par[2])
 396     eff=(bg-par[3])*(bg-par[3])+par[4];
 397   else
 398     eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 399
 400   Double_t bb=0.;
 401   if(gamma>=0. && beta>0.){
 402     bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 403   }
 404   return bb;
 405 }
 406
 407 Double_t BetheBlochPhobosBG(const Double_t *x, const Double_t *par){
 408
 409   //
 410   // returns AliExternalTrackParam::BetheBloch normalized to
 411   // fgMIP at the minimum
 412   //
 413
 414   Double_t bg=x[0];
 415   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 416   Double_t gamma=bg/beta;
 417   Double_t eff=1.0;
 418   if(bg<par[2])
 419     eff=(bg-par[3])*(bg-par[3])+par[4];
 420   else
 421     eff=(par[2]-par[3])*(par[2]-par[3])+par[4];
 422
 423   Double_t bb=0.;
 424   if(gamma>=0. && beta>0.){
 425     bb=(par[1]+2.0*TMath::Log(gamma)-beta*beta)*(par[0]/(beta*beta))*eff;
 426   }
 427   return bb;
 428 }
 429
 430 Double_t BetheBlochAleph(const Double_t *x, const Double_t *par){
 431
 432   //
 433   // This is the empirical ALEPH parameterization of the Bethe-Bloch formula.
 434   // It is normalized to 1 at the minimum.
 435   //
 436   // bg - beta*gamma
 437   //
 438   // The default values for the kp* parameters are for ALICE TPC.
 439   // The returned value is in MIP units
 440   //
 441
 442   Double_t massPion=TDatabasePDG::Instance()->GetParticle(211)->Mass();
 443   Double_t bg=x[0]/massPion;
 444   Double_t beta = bg/TMath::Sqrt(1.+ bg*bg);
 445
 446   Double_t aa = TMath::Power(beta,par[3]);
 447   Double_t bb = TMath::Power(1./bg,par[4]);
 448
 449   bb=TMath::Log(par[2]+bb);
 450
 451   return (par[1]-aa-bb)*par[0]/aa;
 452 }