PWGLF/SPECTRA/IdentifiedHighPt/lib/AliHighPtDeDxBase.cxx

   1 #include "AliHighPtDeDxBase.h"
   2
   3 #include <TMath.h>
   4 #include <TStyle.h>
   5 #include <TROOT.h>
   6
   7 #ifndef __IOSTREAM__
   8 #include <iostream>
   9 #endif
  10
  11 using namespace std;
  12
  13 ClassImp(AliHighPtDeDxBase);
  14
  15 //
  16 // AliHighPtDeDxBase class
  17 //
  18 // This class contains the AliHighPtDeDxBase information
  19 //
  20
  21 //_________________________________________________________
  22 AliHighPtDeDxBase::AliHighPtDeDxBase():
  23   TNamed(),
  24   fIsMc(kFALSE),
  25   fUseRunCut(kFALSE),
  26   fRun(0),
  27   fUseEtaCut(kFALSE),
  28   fUseEtaCutAbs(kFALSE),
  29   fEtaLow(-0.8),
  30   fEtaHigh(0.8),
  31   fUseFilterCut(kFALSE),
  32   fFilter(1),
  33   fUsePhiCut(kFALSE),
  34   fPhiCutLow(0x0),
  35   fPhiCutHigh(0x0),
  36   fEventVtxStatus(-999),
  37   fEventVtxStatusMc(-999),
  38   fEventRun(0),
  39   fEventMag(0),
  40   fEventTrigger(-1),
  41   fTrackCharge(0),
  42   fTrackEta(0),
  43   fTrackP(0),
  44   fTrackPt(0),
  45   fTrackFilter(-999),
  46   fTrackPhi(-999),
  47   fTrackPhiCutVariable(-999),
  48   fTrackDeDx(-999),
  49   fTrackNcl(-999),
  50   fTrackPidMc(0),
  51   fTrackPrimaryMc(-1),
  52   hVtxStatus(0x0),
  53   hNevents(0x0),
  54   hPhi(0x0),
  55   hEta(0x0),
  56   hPt(0x0),
  57   hMeanPt(0x0),
  58   hNclVsPhiVsPtBefore(0x0),
  59   hNclVsPhiVsPtAfter(0x0)
  60 {
  61   // default constructor - do not use
  62 }
  63
  64 //_________________________________________________________
  65 AliHighPtDeDxBase::AliHighPtDeDxBase(const char* name, const char* title):
  66   TNamed(name, title),
  67   fIsMc(kFALSE),
  68   fUseRunCut(kFALSE),
  69   fRun(0),
  70   fUseEtaCut(kFALSE),
  71   fUseEtaCutAbs(kFALSE),
  72   fEtaLow(-0.8),
  73   fEtaHigh(0.8),
  74   fUseFilterCut(kFALSE),
  75   fFilter(1),
  76   fUsePhiCut(kFALSE),
  77   fPhiCutLow(0x0),
  78   fPhiCutHigh(0x0),
  79   fEventVtxStatus(-999),
  80   fEventVtxStatusMc(-999),
  81   fEventRun(0),
  82   fEventMag(0),
  83   fEventTrigger(-1),
  84   fTrackCharge(0),
  85   fTrackEta(0),
  86   fTrackP(0),
  87   fTrackPt(0),
  88   fTrackFilter(-999),
  89   fTrackPhi(-999),
  90   fTrackPhiCutVariable(-999),
  91   fTrackDeDx(-999),
  92   fTrackNcl(-999),
  93   fTrackPidMc(0),
  94   fTrackPrimaryMc(-1),
  95   hVtxStatus(0x0),
  96   hNevents(0x0),
  97   hPhi(0x0),
  98   hEta(0x0),
  99   hPt(0x0),
 100   hMeanPt(0x0),
 101   hNclVsPhiVsPtBefore(0x0),
 102   hNclVsPhiVsPtAfter(0x0)
 103 {
 104   // named constructor
 105 }
 106
 107 //_________________________________________________________
 108 AliHighPtDeDxBase::~AliHighPtDeDxBase()
 109 {
 110   delete hVtxStatus;
 111   delete hNevents;
 112   delete hEta;
 113   delete hPhi;
 114   delete hPt;
 115   delete hMeanPt;
 116   delete hNclVsPhiVsPtBefore;
 117   delete hNclVsPhiVsPtAfter;
 118 }
 119
 120 //_________________________________________________________
 121 void AliHighPtDeDxBase::Init(Int_t nPtBins, Double_t* ptBins)
 122 {
 123   //
 124   // Create histograms
 125   //
 126   hVtxStatus = new TH1D("hVtxStatus", "Vtx status - No Vtx = -1, Vtx outside cut = 0, Vtx inside = 1",
 127                         3, -1.5, 1.5);
 128   hVtxStatus->Sumw2();
 129   hVtxStatus->SetDirectory(0);
 130
 131   hNevents = new TH1D("hNevents", "N events - No Vtx = 0, Vtx OK = 1",
 132                       2, 0, 2);
 133   hNevents->Sumw2();
 134   hNevents->SetDirectory(0);
 135
 136   hEta = new TH1D("hEta", "#eta distribution; #eta; Counts",
 137                   100, -1.0, 1.0);
 138   hEta->Sumw2();
 139   hEta->SetDirectory(0);
 140
 141   hPhi = new TH1D("hPhi", "#varphi distribution; #varphi; Counts",
 142                   360, 0, TMath::TwoPi());
 143   hPhi->Sumw2();
 144   hPhi->SetDirectory(0);
 145
 146   hPt = new TH1D("hPt", "p_{T} spectrum; p_{T} [GeV/c]; Counts",
 147                  nPtBins, ptBins);
 148   hPt->Sumw2();
 149   hPt->SetDirectory(0);
 150
 151   hMeanPt = new TProfile("hMeanPt", "mean p_{T}; p_{T} [GeV/c]; mean p_{T}",
 152                          nPtBins, ptBins);
 153   hMeanPt->SetDirectory(0);
 154
 155   const Int_t nPhiBins = 50;
 156   const Double_t phiBinSize = TMath::Pi()/9.0/nPhiBins;
 157   Double_t phiBins[nPhiBins+1];
 158
 159   for(Int_t i = 0; i <= nPhiBins; i++) {
 160
 161     phiBins[i] = phiBinSize*i;
 162   }
 163
 164   const Int_t nNclBins = 45;
 165   const Double_t nclBinSize = 90.0/nNclBins;
 166   Double_t nclBins[nNclBins+1];
 167
 168   for(Int_t i = 0; i <= nNclBins; i++) {
 169
 170     nclBins[i] = 69.5 + nclBinSize*i;
 171   }
 172
 173   hNclVsPhiVsPtBefore = new TH3F("hNclVsPhiVsPtBefore", "<Ncl> vs Pt and #phi (before #phi cut); p_{T} [GeV/c]; fmod(#phi, #pi/9.0)",
 174                                  nPtBins, ptBins, nPhiBins, phiBins, nNclBins, nclBins);
 175   hNclVsPhiVsPtBefore->SetDirectory(0);
 176   hNclVsPhiVsPtAfter = new TH3F("hNclVsPhiVsPtAfter", "<Ncl> vs Pt and #phi (after #phi cut); p_{T} [GeV/c]; fmod(#phi, #pi/9.0)",
 177                                 nPtBins, ptBins, nPhiBins, phiBins, nNclBins, nclBins);
 178   hNclVsPhiVsPtAfter->SetDirectory(0);
 179 }
 180
 181 //_________________________________________________________
 182 Bool_t AliHighPtDeDxBase::EventAccepted()
 183 {
 184   if(fUseRunCut && fRun!=fEventRun)
 185     return kFALSE;
 186
 187   return kTRUE;
 188 }
 189
 190 //_________________________________________________________
 191 Bool_t AliHighPtDeDxBase::TrackAccepted()
 192 {
 193   if(fUseFilterCut && !(fTrackFilter&fFilter))
 194     return kFALSE;
 195
 196   if(fUseEtaCut && (fTrackEta<fEtaLow || fTrackEta>fEtaHigh))
 197     return kFALSE;
 198
 199   if(fUseEtaCutAbs && (TMath::Abs(fTrackEta)<fEtaLow || TMath::Abs(fTrackEta)>fEtaHigh))
 200     return kFALSE;
 201
 202   if(fUsePhiCut) {
 203
 204     fTrackPhiCutVariable = fTrackPhi;
 205     if(fEventMag<0)    // for negatve polarity field
 206       fTrackPhiCutVariable = TMath::TwoPi() - fTrackPhiCutVariable;
 207     if(fTrackCharge<0) // for negatve charge
 208       fTrackPhiCutVariable = TMath::TwoPi()-fTrackPhiCutVariable;
 209     if(fTrackPhiCutVariable < 0)
 210       cout << "Warning!!!!! phi < 0: " << fTrackPhiCutVariable << endl;
 211
 212     fTrackPhiCutVariable += TMath::Pi()/18.0; // to center gap in the middle
 213     fTrackPhiCutVariable = fmod(fTrackPhiCutVariable, TMath::Pi()/9.0);
 214
 215     hNclVsPhiVsPtBefore->Fill(fTrackPt, fTrackPhiCutVariable, fTrackNcl);
 216
 217     if(fTrackPt>2.0 && fTrackPhiCutVariable<fPhiCutHigh->Eval(fTrackPt)
 218        && fTrackPhiCutVariable>fPhiCutLow->Eval(fTrackPt))
 219       return kFALSE; // reject track
 220
 221     hNclVsPhiVsPtAfter->Fill(fTrackPt, fTrackPhiCutVariable, fTrackNcl);
 222   }
 223
 224   return kTRUE;
 225 }
 226
 227 //_________________________________________________________
 228 void AliHighPtDeDxBase::FillEventInfo()
 229 {
 230   if(fEventTrigger==1) {
 231     hVtxStatus->Fill(fEventVtxStatus);
 232     if(fEventVtxStatus != 0) {
 233       hNevents->Fill(1.0+0.5*fEventVtxStatus);
 234     }
 235   }
 236 }
 237
 238 //_________________________________________________________
 239 void AliHighPtDeDxBase::FillTrackInfo(Float_t weight)
 240 {
 241   hEta->Fill(fTrackEta, weight);
 242   hPhi->Fill(fTrackPhi, weight);
 243   hPt->Fill(fTrackPt, weight);
 244   hMeanPt->Fill(fTrackPt, fTrackPt);
 245 }
 246
 247 //_________________________________________________________
 248 void AliHighPtDeDxBase::Print(Option_t* option) const
 249 {
 250   cout << ClassName() << " : " << GetName() << endl
 251        << "Event cuts: " << endl;
 252   if(fUseRunCut)
 253     cout << "   Run cut: " << fRun << endl;
 254   else
 255     cout << "   Run cut is diabled " << endl;
 256
 257   cout << "Track cuts: " << endl;
 258   if(fUseFilterCut)
 259     cout << "   Filter cut: " << fFilter << endl;
 260   else
 261     cout << "   Filter cut is diabled " << endl;
 262   if(fUseEtaCut)
 263     cout << "   Eta range: " << fEtaLow << " - " << fEtaHigh << endl;
 264   if(fUseEtaCutAbs)
 265     cout << "   |Eta| range: " << fEtaLow << " - " << fEtaHigh << endl;
 266   else
 267     cout << "   Eta cut is diabled " << endl;
 268   if(fUsePhiCut)
 269     cout << "   Phi cut is ENABLED" << endl;
 270   else
 271     cout << "   Phi cut is diabled " << endl;
 272 }
 273
 274 //_________________________________________________________
 275 void AliHighPtDeDxBase::SetUseEtaCut(Bool_t  value)
 276 {
 277   fUseEtaCut = value;
 278
 279   if(value == kTRUE)
 280     fUseEtaCutAbs = kFALSE;
 281 }
 282
 283 //_________________________________________________________
 284 void AliHighPtDeDxBase::SetUseEtaCutAbs(Bool_t  value)
 285 {
 286   fUseEtaCutAbs = value;
 287
 288   if(value == kTRUE)
 289     fUseEtaCut = kFALSE;
 290 }
 291
 292 //_________________________________________________________
 293 TF1* AliHighPtDeDxBase::GetStandardPhiCutLow()
 294 {
 295   //  TF1* cutLow  = new TF1("StandardPhiCutLow",  "-0.01/x+pi/18.0-0.015", 0, 50);
 296   TF1* cutLow  = new TF1("StandardPhiCutLow",  "0.1/x/x+pi/18.0-0.025", 0, 50);
 297   return cutLow;
 298 }
 299
 300 //_________________________________________________________
 301 TF1* AliHighPtDeDxBase::GetStandardPhiCutHigh()
 302 {
 303   //  TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.55/x/x+pi/18.0+0.03", 0, 50);
 304   TF1* cutHigh = new TF1("StandardPhiCutHigh", "0.12/x+pi/18.0+0.035", 0, 50);
 305   return cutHigh;
 306 }
 307
 308
 309 //___________________________________________________________________________
 310 TCanvas* AliHighPtDeDxBase::DrawPhiCutHistograms()
 311 {
 312   gStyle->SetOptStat(0);
 313
 314   TCanvas* cPhiCut = FindCanvas("cPhiCut", 1200, 800);
 315   cPhiCut->SetTitle("phi cut histograms");
 316   cPhiCut->Clear();
 317   cPhiCut->Divide(3,2);
 318
 319   // cPhiCut->cd(1);
 320   // TH2D* hPhiVsPtBefore = (TH2D*)hNclVsPhiVsPtBefore->Project3D("yx");
 321   // hPhiVsPtBefore->SetName("hPhiVsPtBefore");
 322   // hPhiVsPtBefore->SetTitle("Phi vs p_{T} (before cuts)");
 323   // MakeNice2dHisto(hPhiVsPtBefore, gPad, kTRUE);
 324   // hPhiVsPtBefore->Draw("COL");
 325   // fPhiCutHigh->SetRange(2.0, 50.0);
 326   // fPhiCutHigh->Draw("SAME");
 327   // fPhiCutLow->SetRange(2.0, 50.0);
 328   // fPhiCutLow->Draw("SAME");
 329
 330   cPhiCut->cd(1);
 331   TProfile2D* hNclBefore = hNclVsPhiVsPtBefore->Project3DProfile("yx");
 332   hNclBefore->SetName("hNclBefore");
 333   hNclBefore->SetTitle("<Ncl> (before cuts); p_{T} [GeV/c]; #varphi'");
 334   MakeNice2dHisto(hNclBefore, gPad, kTRUE);
 335   DrawNice(hNclBefore, 0, 0, 0, "COLZ");
 336   fPhiCutHigh->Draw("SAME");
 337   fPhiCutLow->Draw("SAME");
 338
 339   cPhiCut->cd(2);
 340   gPad->SetGridy();
 341   TH2D* hNclVsPtBefore = (TH2D*)hNclVsPhiVsPtBefore->Project3D("zx");
 342   hNclVsPtBefore->SetName("hNclVsPtBefore");
 343   hNclVsPtBefore->SetTitle("; p_{T} [GeV/c]; Ncl;");
 344   MakeNice2dHisto(hNclVsPtBefore, gPad, kTRUE);
 345   hNclVsPtBefore->Draw("COL");
 346   TProfile* hNclVsPtBeforeProf = hNclVsPtBefore->ProfileX();
 347   hNclVsPtBeforeProf->SetMarkerStyle(29);
 348   hNclVsPtBeforeProf->SetMarkerColor(2);
 349   hNclVsPtBeforeProf->Draw("SAME");
 350
 351   // cPhiCut->cd(2);
 352   // TH2D* hPhiVsPtAfter = (TH2D*)hNclVsPhiVsPtAfter->Project3D("yx");
 353   // MakeNice2dHisto(hNclVsPtAfter, gPad, kTRUE);
 354   // hPhiVsPtAfter->SetName("hPhiVsPtAfter");
 355   // hPhiVsPtAfter->SetTitle("Phi vs p_{T} (after cuts)");
 356   // DrawNice(hPhiVsPtAfter, 0, 0, 0, "COL");
 357
 358   cPhiCut->cd(4);
 359   TProfile2D* hNclAfter = (TProfile2D*)hNclVsPhiVsPtAfter->Project3DProfile("yx");
 360   hNclAfter->SetName("hNclAfter");
 361   hNclAfter->SetTitle("<Ncl> (after cuts); p_{T} [GeV/c]; #varphi'");
 362   MakeNice2dHisto(hNclAfter, gPad, kTRUE);
 363   DrawNice(hNclAfter, 0, 0, 0, "COLZ");
 364
 365   cPhiCut->cd(5);
 366   gPad->SetGridy();
 367   TH2D* hNclVsPtAfter = (TH2D*)hNclVsPhiVsPtAfter->Project3D("zx");
 368   hNclVsPtAfter->SetName("hNclVsPtAfter");
 369   hNclVsPtAfter->SetTitle("; p_{T} [GeV/c]; Ncl;");
 370   MakeNice2dHisto(hNclVsPtAfter, gPad, kTRUE);
 371   hNclVsPtAfter->Draw("COL");
 372   TProfile* hNclVsPtAfterProf = hNclVsPtAfter->ProfileX();
 373   hNclVsPtAfterProf->SetMarkerStyle(29);
 374   hNclVsPtAfterProf->SetMarkerColor(2);
 375   hNclVsPtAfterProf->Draw("SAME");
 376
 377   cPhiCut->cd(3);
 378   gPad->SetGridy();
 379   TH1D* hEfficiency = (TH1D*)hNclVsPhiVsPtAfter->Project3D("x");
 380   hEfficiency->SetName("hEfficiency");
 381   hEfficiency->SetTitle("; p_{T} [GeV/c]; Cut efficiency;");
 382   TH1D* hHelp = (TH1D*)hNclVsPhiVsPtBefore->Project3D("x");
 383   hEfficiency->Divide(hEfficiency, hHelp, 1, 1, "B");
 384   delete hHelp;
 385   MakeNice1dHisto(hEfficiency, gPad);
 386   hEfficiency->SetMarkerStyle(20);
 387   hEfficiency->Draw("P");
 388
 389   return cPhiCut;
 390 }
 391
 392 /*
 393   These methods does not fit in very well here, but for now I put them
 394   here. Maybe they should go in some tool class.
 395  */
 396
 397 //_______________________________________________________________________
 398 void AliHighPtDeDxBase::MakeNice1dHisto(TH1* hist, TVirtualPad* c1)
 399 {
 400   if(c1) {
 401
 402     c1->SetLeftMargin(0.12);
 403     c1->SetRightMargin(0.02);
 404     c1->SetBottomMargin(0.12);
 405     c1->SetTopMargin(0.06);
 406   }
 407
 408   gStyle->SetTitleH(0.08);
 409   gStyle->SetTitleW(0.46);
 410   gStyle->SetTitleX(0.29);
 411   gStyle->SetTitleY(0.99);
 412   hist->GetXaxis()->SetTitleSize(0.07);
 413   hist->GetYaxis()->SetTitleSize(0.07);
 414   hist->GetXaxis()->SetTitleOffset(0.7);
 415   hist->GetYaxis()->SetTitleOffset(0.8);
 416 }
 417
 418 //_______________________________________________________________________
 419 void AliHighPtDeDxBase::MakeNice2dHisto(TH2* hist, TVirtualPad* c1, Bool_t colz)
 420 {
 421   gStyle->SetTitleH(0.08);
 422   gStyle->SetTitleW(0.76);
 423   gStyle->SetTitleX(0.14);
 424   gStyle->SetTitleY(0.98);
 425   if(colz)
 426     gStyle->SetTitleW(0.68);
 427
 428   hist->GetXaxis()->SetTitleSize(0.07);
 429   hist->GetYaxis()->SetTitleSize(0.07);
 430   hist->GetXaxis()->SetTitleOffset(0.7);
 431   hist->GetYaxis()->SetTitleOffset(0.8);
 432   if(c1) {
 433
 434     c1->SetLeftMargin(0.12);
 435     if(!colz)
 436       c1->SetRightMargin(0.02);
 437     else
 438       c1->SetRightMargin(0.12);
 439     c1->SetBottomMargin(0.12);
 440     c1->SetTopMargin(0.06);
 441   }
 442 }
 443
 444 //______________________________________________________________________
 445 TCanvas* AliHighPtDeDxBase::FindCanvas(const Char_t* canvasName,
 446                                        Int_t xwidth, Int_t ywidth)
 447 {
 448
 449   TCanvas *c1 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject(canvasName);
 450   if(!c1)
 451     c1 = new TCanvas(canvasName, canvasName, xwidth, ywidth);
 452   else
 453     c1->SetWindowSize(xwidth, ywidth);
 454
 455   return c1;
 456 }
 457
 458 //_______________________________________________________________________
 459 TCanvas* AliHighPtDeDxBase::DrawNice(TH1* hist, const Char_t* canvasName,
 460                                      Int_t xwidth, Int_t ywidth, const Char_t* option)
 461 {
 462   TCanvas* canv = 0;
 463
 464   if(canvasName) {
 465     canv = FindCanvas(canvasName, xwidth, ywidth);
 466     canv->Clear();
 467   }
 468
 469   Bool_t is2dHist = hist->InheritsFrom("TH2");
 470
 471   if(is2dHist) {
 472
 473     TString opt(option);
 474     opt.ToUpper();
 475
 476     if (opt.Contains("COL"))
 477       MakeNice2dHisto((TH2*)hist, canv, kTRUE);
 478     else
 479       MakeNice2dHisto((TH2*)hist, canv);
 480
 481   } else
 482     MakeNice1dHisto(hist, canv);
 483
 484   hist->Draw(option);
 485
 486   return canv;
 487 }