PWGLF/SPECTRA/PiKaPr/TOF/PbPb276/macros/TPCcalib.C

   1 Double_t tofReso = 85.;
   2
   3 /**************************************************************/
   4 /*** HISTOS AND BINNING ***************************************/
   5 /**************************************************************/
   6
   7 /**************************************************************/
   8 enum ECharge_t {
   9   kPositive,
  10   kNegative,
  11   kNCharges
  12 };
  13 const Char_t *chargeName[kNCharges] = {
  14   "positive",
  15   "negative",
  16 };
  17 /**************************************************************/
  18 enum EHistoParam_t {
  19   kCentrality,
  20   kPtpc,
  21   kTPCNcls,
  22   kTPCdelta,
  23   kTPCgain,
  24   kTOFsignal,
  25   kNHistoParams
  26 };
  27 /**************************************************************/
  28 const Int_t NcentralityBins = 10;
  29 Double_t centralityBin[NcentralityBins + 1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90.};
  30 /**************************************************************/
  31 const Int_t NpBins = 46;
  32 Double_t pBin[NpBins + 1] = {0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0};
  33 /**************************************************************/
  34 const Int_t NclsBins = 200;
  35 Double_t clsBin[NclsBins + 1];
  36 Double_t clsMin = 0., clsMax = 200., clsStep = (clsMax - clsMin) / NclsBins;
  37 /**************************************************************/
  38 const Int_t NtpcdeltaBins = 200;
  39 Double_t tpcdeltaBin[NtpcdeltaBins + 1];
  40 Double_t tpcdeltaMin = -100., tpcdeltaMax = 100., tpcdeltaStep = (tpcdeltaMax - tpcdeltaMin) / NtpcdeltaBins;
  41 /**************************************************************/
  42 const Int_t NtpcgainBins = 200;
  43 Double_t tpcgainBin[NtpcgainBins + 1];
  44 Double_t tpcgainMin = 0., tpcgainMax = 2., tpcgainStep = (tpcgainMax - tpcgainMin) / NtpcgainBins;
  45 /**************************************************************/
  46 const Int_t NsigmaBins = 20;
  47 Double_t sigmaBin[NsigmaBins + 1];
  48 Double_t sigmaMin = -5., sigmaMax = 5., sigmaStep = (sigmaMax - sigmaMin) / NsigmaBins;
  49 /**************************************************************/
  50 Int_t NparamsBins[kNHistoParams] = {NcentralityBins, NpBins, NclsBins, NtpcdeltaBins, NtpcgainBins, NsigmaBins};
  51 Double_t *paramBin[kNHistoParams] = {centralityBin, pBin, clsBin, tpcdeltaBin, tpcgainBin, sigmaBin};
  52 /**************************************************************/
  53
  54 /**************************************************************/
  55
  56
  57 TPCcalib(const Char_t *filename, Int_t evMax = kMaxInt, Int_t startEv = 0)
  58 {
  59
  60   /* include path for ACLic */
  61   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  62   gSystem->AddIncludePath("-I$ALICE_ROOT/TOF");
  63   /* load libraries */
  64   gSystem->Load("libANALYSIS");
  65   gSystem->Load("libANALYSISalice");
  66   /* build analysis task class */
  67   gROOT->LoadMacro("AliAnalysisParticle.cxx+g");
  68   gROOT->LoadMacro("AliAnalysisEvent.cxx+g");
  69   gROOT->LoadMacro("AliAnalysisTrack.cxx+g");
  70
  71   /* open file, get tree and connect */
  72   TFile *filein = TFile::Open(filename);
  73   TTree *treein = (TTree *)filein->Get("aodTree");
  74   printf("got \"aodTree\": %d entries\n", treein->GetEntries());
  75   AliAnalysisEvent *analysisEvent = new AliAnalysisEvent();
  76   TClonesArray *analysisTrackArray = new TClonesArray("AliAnalysisTrack");
  77   AliAnalysisTrack *analysisTrack = NULL;
  78   treein->SetBranchAddress("AnalysisEvent", &analysisEvent);
  79   treein->SetBranchAddress("AnalysisTrack", &analysisTrackArray);
  80
  81   /**************************************************************/
  82   /*** HISTOS ***************************************************/
  83   /**************************************************************/
  84
  85   /* run-time binning */
  86   for (Int_t ibin = 0; ibin < NclsBins + 1; ibin++)
  87     clsBin[ibin] = clsMin + ibin * clsStep;
  88   for (Int_t ibin = 0; ibin < NtpcdeltaBins + 1; ibin++)
  89     tpcdeltaBin[ibin] = tpcdeltaMin + ibin * tpcdeltaStep;
  90   for (Int_t ibin = 0; ibin < NtpcgainBins + 1; ibin++)
  91     tpcgainBin[ibin] = tpcgainMin + ibin * tpcgainStep;
  92   for (Int_t ibin = 0; ibin < NsigmaBins + 1; ibin++)
  93     sigmaBin[ibin] = sigmaMin + ibin * sigmaStep;
  94
  95   /* THnSparse */
  96   THnSparse *hTPCcalib[AliPID::kSPECIES][kNCharges];
  97   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
  98     for (Int_t icharge = 0; icharge< kNCharges; icharge++) {
  99       hTPCcalib[ipart][icharge] = new THnSparseF(Form("hTPCcalib_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]), "", kNHistoParams, NparamsBins);
 100       for (Int_t iparam = 0; iparam < kNHistoParams; iparam++) {
 101         hTPCcalib[ipart][icharge]->SetBinEdges(iparam, paramBin[iparam]);
 102       }
 103     }
 104   }
 105
 106   /**************************************************************/
 107   /**************************************************************/
 108   /**************************************************************/
 109
 110   /* TOF PID response */
 111   AliTOFPIDResponse tofResponse;
 112   tofResponse.SetTimeResolution(tofReso);
 113   /* TPC PID response */
 114   AliTPCPIDResponse *tpcResponse = AliAnalysisTrack::GetTPCResponse();
 115
 116   /* start stopwatch */
 117   TStopwatch timer;
 118   timer.Start();
 119
 120   /* loop over events */
 121   Bool_t hastofpid;
 122   Int_t charge, index;
 123   UShort_t dedxN;
 124   Double_t ptpc, dedx, bethe, deltadedx, dedx_sigma, tpcsignal;
 125   Double_t p, time, time_sigma, timezero, timezero_sigma, tof, tof_sigma, texp, texp_sigma, deltat, deltat_sigma, tofsignal;
 126   Double_t tpctofsignal;
 127   Double_t param[kNHistoParams];
 128   for (Int_t iev = startEv; iev < treein->GetEntries() && iev < evMax; iev++) {
 129     /* get event */
 130     treein->GetEvent(iev);
 131     if (iev % 100 == 0) printf("iev = %d\n", iev);
 132     /* check vertex */
 133     if (!analysisEvent->AcceptVertex()) continue;
 134     /* check collision candidate */
 135     if (!analysisEvent->IsCollisionCandidate()) continue;
 136     /* check centrality quality */
 137     if (analysisEvent->GetCentralityQuality() != 0.) continue;
 138
 139     /*** ACCEPTED EVENT ***/
 140
 141     /* apply time-zero TOF correction */
 142     analysisEvent->ApplyTimeZeroTOFCorrection();
 143
 144     /* get centrality */
 145     param[kCentrality] = analysisEvent->GetCentralityPercentile(AliAnalysisEvent::kCentEst_V0M);
 146
 147     /* loop over tracks */
 148     for (Int_t itrk = 0; itrk < analysisTrackArray->GetEntries(); itrk++) {
 149       /* get track */
 150       analysisTrack = (AliAnalysisTrack *)analysisTrackArray->At(itrk);
 151       if (!analysisTrack) continue;
 152       /* check accepted track */
 153       if (!analysisTrack->AcceptTrack()) continue;
 154       /* get charge */
 155       charge = analysisTrack->GetSign() > 0. ? kPositive : kNegative;
 156
 157       /*** ACCEPTED TRACK ***/
 158
 159       /* check TOF pid */
 160       if (!analysisTrack->HasTOFPID() || !analysisTrack->HasTPCPID()) continue;
 161
 162       /*** ACCEPTED TRACK WITH TPC+TOF PID ***/
 163
 164       /* apply expected time correction */
 165       analysisTrack->ApplyTOFExpectedTimeCorrection();
 166
 167       /* get track info */
 168       p = analysisTrack->GetP();
 169
 170       /* get TPC info */
 171       dedx = analysisTrack->GetTPCdEdx();
 172       dedxN = analysisTrack->GetTPCdEdxN();
 173       ptpc = analysisTrack->GetTPCmomentum();
 174       param[kPtpc] = ptpc;
 175       param[kTPCNcls] = dedxN;
 176
 177       /* get TOF info */
 178       time = analysisTrack->GetTOFTime();
 179       time_sigma = tofReso;
 180       timezero = analysisEvent->GetTimeZeroTOF(p);
 181       timezero_sigma = analysisEvent->GetTimeZeroTOFSigma(p);
 182       tof = time - timezero;
 183       tof_sigma = TMath::Sqrt(time_sigma * time_sigma + timezero_sigma * timezero_sigma);
 184
 185       /* loop over particle IDs */
 186       for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
 187
 188         /* check rapidity */
 189         if (TMath::Abs(analysisTrack->GetY(AliPID::ParticleMass(ipart))) > 0.5) continue;
 190
 191         /*** ACCEPTED TRACK WITHIN CORRECT RAPIDITY ***/
 192
 193         /* TOF expected time */
 194         texp = analysisTrack->GetTOFExpTime(ipart);
 195         texp_sigma = analysisTrack->GetTOFExpTimeSigma(ipart);
 196
 197         /* TOF signal */
 198         deltat = tof - texp;
 199         deltat_sigma = TMath::Sqrt(tof_sigma * tof_sigma + texp_sigma * texp_sigma);
 200         tofsignal = deltat / deltat_sigma;
 201         param[kTOFsignal] = tofsignal;
 202
 203         /* check TOF PID with 3-sigma cut */
 204         if (TMath::Abs(tofsignal) > 3.) continue;
 205
 206         /*** ACCEPTED TRACK WITH COMPATIBLE TOF PID ***/
 207
 208         /* TPC signal */
 209         bethe = tpcResponse->GetExpectedSignal(ptpc, ipart);
 210         deltadedx = dedx - bethe;
 211         param[kTPCdelta] = deltadedx;
 212         param[kTPCgain] = dedx / bethe;
 213
 214         /* fill histo */
 215         hTPCcalib[ipart][charge]->Fill(param);
 216
 217       } /* end of loop over particle IDs */
 218     } /* end of loop over tracks */
 219   } /* end of loop over events */
 220
 221   /* start stopwatch */
 222   timer.Stop();
 223   timer.Print();
 224
 225   /* output */
 226   TFile *fileout = TFile::Open(Form("TPCcalib.%s", filename), "RECREATE");
 227   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
 228     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 229       hTPCcalib[ipart][icharge]->Write();
 230     }
 231   }
 232   fileout->Close();
 233
 234 }
 235
 236 /**************************************************************/
 237
 238 TPCcalib_delta(const Char_t *filename, Int_t ipart, Int_t icharge, Float_t tofsigmaMin = -3., Float_t tofsigmaMax = 3.)
 239 {
 240
 241   Double_t ptMin[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.5};
 242   Double_t ptMax[AliPID::kSPECIES] = {0.5, 0., 1.5, 1.5, 2.0};
 243
 244   /* load HistoUtils */
 245   gROOT->LoadMacro("HistoUtils.C");
 246
 247   /* open data */
 248   TFile *filein = TFile::Open(filename);
 249   THnSparseF *hTPCcalib = (THnSparseF *)filein->Get(Form("hTPCcalib_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 250
 251   /* set momentum range */
 252   hTPCcalib->GetAxis(kPtpc)->SetRangeUser(ptMin[ipart] + 0.001, ptMax[ipart] - 0.001);
 253
 254   /* set TOF pid range */
 255   //  hTPCcalib->GetAxis(kTOFsignal)->SetRangeUser(tofsigmaMin + 0.001, tofsigmaMax - 0.001);
 256
 257   /* minimum-bias projection */
 258   hTPCcalib->GetAxis(kCentrality)->SetRange(1, NcentralityBins);
 259   TH2 *hDeltaPtpc = hTPCcalib->Projection(kTPCdelta, kPtpc);
 260   TObjArray *oaTPCdelta_mb = HistoUtils_FitPeak(NULL, hDeltaPtpc, 5, 3., 3., 100, Form("hTPCdelta_mb_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 261   delete hDeltaPtpc;
 262
 263   /* convert to bethe-blok */
 264   TH1D *hTPCdelta_mb = (TH1D *)oaTPCdelta_mb->At(1);
 265   TH1D *hTPCbethe_mb = TPCcalib_delta2bethe(hTPCdelta_mb, ipart, Form("hTPCbethe_mb_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 266
 267   /* centrality projections */
 268   TObjArray *oaTPCdelta_cent[NcentralityBins];
 269   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 270
 271     hTPCcalib->GetAxis(kCentrality)->SetRange(icent + 1, icent + 1);
 272     hDeltaPtpc = hTPCcalib->Projection(kTPCdelta, kPtpc);
 273     oaTPCdelta_cent[icent] = HistoUtils_FitPeak(NULL, hDeltaPtpc, 5, 3., 3., 100, Form("hTPCdelta_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 274     delete hDeltaPtpc;
 275   }
 276
 277 }
 278
 279 /**************************************************************/
 280
 281 TF1 *
 282 TPCcalib_BetheBlochAleph(Double_t min, Double_t max)
 283 {
 284
 285   TF1 *f = new TF1("fBetheBlochAleph", "[0] * AliExternalTrackParam::BetheBlochAleph(x, [1], [2], [3], [4], [5])", min, max);
 286   f->SetParameter(0, 50.);
 287   f->SetParameter(1, 0.0283086);
 288   f->SetParameter(2, 2.63394e+01);
 289   f->SetParameter(3, 5.04114e-11);
 290   f->SetParameter(4, 2.12543);
 291   f->SetParameter(5, 4.88663);
 292   return f;
 293
 294 };
 295
 296 /**************************************************************/
 297
 298 void
 299 TPCcalib_fitBetheBloch(TGraphErrors *g, TF1 *f, Int_t param = -1)
 300 {
 301
 302   if (param >= 0) {
 303     for (Int_t ipar = 0; ipar < 6; ipar++)
 304       if (ipar != param)
 305         f->FixParameter(ipar, f->GetParameter(ipar));
 306       else
 307         f->ReleaseParameter(ipar);
 308   }
 309   else {
 310     for (Int_t ipar = 0; ipar < 6; ipar++)
 311       f->ReleaseParameter(ipar);
 312   }
 313   g->Fit(f);
 314
 315 }
 316
 317 /**************************************************************/
 318
 319 TPCcalib_bethe(const Char_t *filename, Int_t icent = -1, Float_t tofsigmaMin = -3., Float_t tofsigmaMax = 3.)
 320 {
 321
 322   Double_t ptMin[AliPID::kSPECIES] = {0., 0., 0.5, 0.5, 0.5};
 323   Double_t ptMax[AliPID::kSPECIES] = {0., 0., 1.5, 1.5, 2.0};
 324
 325   /* set centrality name */
 326   Char_t centName[1024];
 327   if (icent < 0 || icent >= NcentralityBins)
 328     sprintf(centName, "cent0090");
 329   else
 330     sprintf(centName, "cent%02d%02d", , centralityBin[icent], centralityBin[icent + 1]);
 331
 332   /* load HistoUtils */
 333   gROOT->LoadMacro("HistoUtils.C");
 334
 335   /* open data */
 336   TFile *filein = TFile::Open(filename);
 337
 338   /* output */
 339   TFile *fileout = TFile::Open(Form("TPCcalib_bethe_%s.%s", centName, filename), "RECREATE");
 340
 341   /* loop over particles and charges */
 342   TH1D *hTPCbethe[AliPID::kSPECIES][kNCharges];
 343   TGraphErrors *gTPCbethe = new TGraphErrors();
 344   gTPCbethe->SetName("gTPCbethe");
 345   Int_t nPoints = 0;
 346   for (Int_t ipart = 0; ipart < AliPID::kSPECIES; ipart++) {
 347     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 348
 349       /* check momentum range */
 350       if (ptMax[ipart] == 0. || ptMax[ipart] < ptMin[ipart]) continue;
 351
 352       /* get histo */
 353       THnSparseF *hTPCcalib = (THnSparseF *)filein->Get(Form("hTPCcalib_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 354
 355       /* set centrality range */
 356       if (icent < 0 || icent >= NcentralityBins)
 357         hTPCcalib->GetAxis(kCentrality)->SetRange(1, NcentralityBins);
 358       else
 359         hTPCcalib->GetAxis(kCentrality)->SetRange(icent + 1, icent + 1);
 360
 361       /* set TOF pid range */
 362       //  hTPCcalib->GetAxis(kTOFsignal)->SetRangeUser(tofsigmaMin + 0.001, tofsigmaMax - 0.001);
 363
 364       /* set momentum range */
 365       hTPCcalib->GetAxis(kPtpc)->SetRangeUser(ptMin[ipart] + 0.001, ptMax[ipart] - 0.001);
 366
 367       /* minimum-bias projection */
 368       TH2 *hDeltaPtpc = hTPCcalib->Projection(kTPCdelta, kPtpc);
 369       TObjArray *oaTPCdelta = HistoUtils_FitPeak(NULL, hDeltaPtpc, 5, 3., 3., 100, Form("hTPCdelta_%s_%s_%s", centName, AliPID::ParticleName(ipart), chargeName[icharge]));
 370       delete hDeltaPtpc;
 371
 372       /* convert to bethe-blok */
 373       TH1D *hTPCdelta = (TH1D *)oaTPCdelta->At(1);
 374       hTPCbethe[ipart][icharge] = TPCcalib_delta2bethe(hTPCdelta, ipart, Form("hTPCbethe_%s_%s_%s", centName, AliPID::ParticleName(ipart), chargeName[icharge]));
 375
 376       /* delete array */
 377       delete oaTPCdelta;
 378
 379       /* add points to graph */
 380       for (Int_t ibin = 0; ibin < hTPCbethe[ipart][icharge]->GetNbinsX(); ibin++) {
 381         gTPCbethe->SetPoint(nPoints, hTPCbethe[ipart][icharge]->GetBinCenter(ibin + 1), hTPCbethe[ipart][icharge]->GetBinContent(ibin + 1));
 382         gTPCbethe->SetPointError(nPoints, 0.5 * hTPCbethe[ipart][icharge]->GetBinWidth(ibin + 1), hTPCbethe[ipart][icharge]->GetBinError(ibin + 1));
 383         nPoints++;
 384       }
 385
 386       /* write */
 387       fileout->cd();
 388       hTPCbethe[ipart][icharge]->Write();
 389     }
 390   }
 391
 392   TF1 *fBethe = TPCcalib_BetheBlochAleph(0.1, 1.e3);
 393   for (Int_t i = 0; i < 50; i++)
 394     for (Int_t ii = 0; ii < 6; ii++)
 395       TPCcalib_fitBetheBloch(gTPCbethe, fBethe, ii);
 396   TPCcalib_fitBetheBloch(gTPCbethe, fBethe, -1);
 397
 398   /* write */
 399   fileout->cd();
 400   gTPCbethe->Write("gTPCbethe");
 401   fBethe->Write("fBetheBlochAleph");
 402
 403   /* close file */
 404   fileout->Close();
 405
 406   return;
 407
 408 }
 409
 410 /**************************************************************/
 411
 412 TPCcalib_gain(const Char_t *filename, Int_t ipart, Int_t icharge, Float_t tofsigmaMin = -3., Float_t tofsigmaMax = 3.)
 413 {
 414
 415   /* load HistoUtils */
 416   gROOT->LoadMacro("HistoUtils.C");
 417
 418   /* open data */
 419   TFile *filein = TFile::Open(filename);
 420   THnSparseF *hTPCcalib = (THnSparseF *)filein->Get(Form("hTPCcalib_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 421
 422   /* set TOF pid range */
 423   //  hTPCcalib->GetAxis(kTOFsignal)->SetRangeUser(tofsigmaMin + 0.001, tofsigmaMax - 0.001);
 424
 425   /* minimum-bias projection */
 426   hTPCcalib->GetAxis(kCentrality)->SetRange(1, NcentralityBins);
 427   TH2 *hGainPtpc = hTPCcalib->Projection(kTPCgain, kPtpc);
 428   TObjArray *oaTPCgain_mb = HistoUtils_FitPeak(NULL, hGainPtpc, 5, 3., 3., 100, Form("hTPCgain_mb_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 429   delete hGainPtpc;
 430   TH1D *hmb = (TH1D *)oaTPCgain_mb->At(1);
 431
 432   /* centrality projections */
 433   TObjArray *oaTPCgain_cent[NcentralityBins];
 434   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 435
 436     hTPCcalib->GetAxis(kCentrality)->SetRange(icent + 1, icent + 1);
 437     hGainPtpc = hTPCcalib->Projection(kTPCgain, kPtpc);
 438     oaTPCgain_cent[icent] = HistoUtils_FitPeak(NULL, hGainPtpc, 5, 3., 3., 100, Form("hTPCgain_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 439     TH1D *hcent = (TH1D *)oaTPCgain_cent[icent]->At(1);
 440     hcent->Divide(hmb);
 441     delete hGainPtpc;
 442   }
 443
 444 }
 445
 446 /**************************************************************/
 447
 448 TH1D *
 449 TPCcalib_delta2bethe(TH1D *hDelta, Int_t ipart, const Char_t *name = "hTPCbethe")
 450 {
 451
 452   Int_t nbins = hDelta->GetNbinsX();
 453   Float_t *xbins = new Float_t[nbins + 1];
 454   for(Int_t ibin = 0; ibin <= nbins; ibin++)
 455     xbins[ibin] = hDelta->GetBinLowEdge(ibin + 1) / AliPID::ParticleMass(ipart);
 456
 457   TH1D *hBethe = new TH1D(name, "", nbins, xbins);
 458   AliTPCPIDResponse tpcResponse;
 459   Double_t ptpc, delta, bethe, deltae;
 460   for (Int_t ibin = 0; ibin < nbins; ibin++) {
 461     ptpc = hDelta->GetBinCenter(ibin + 1);
 462     delta = hDelta->GetBinContent(ibin + 1);
 463     deltae = hDelta->GetBinError(ibin + 1);
 464     bethe = tpcResponse.GetExpectedSignal(ptpc, ipart);
 465     hBethe->SetBinContent(ibin + 1, bethe + delta);
 466     hBethe->SetBinError(ibin + 1, deltae);
 467   }
 468
 469   return hBethe;
 470 }