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

   1 /**************************************************************/
   2 enum ECharge_t {
   3   kPositive,
   4   kNegative,
   5   kNCharges
   6 };
   7 const Char_t *chargeName[kNCharges] = {
   8   "positive",
   9   "negative",
  10 };
  11 /**************************************************************/
  12 const Int_t NcentralityBins = 10;
  13 Double_t centralityBin[NcentralityBins + 1] = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80., 90.};
  14 /**************************************************************/
  15 const Int_t NptBins = 46;
  16 Double_t ptBin[NptBins + 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};
  17 /**************************************************************/
  18 Char_t *partLatex[AliPID::kSPECIES][2] = {
  19   "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
  20 };
  21  /**************************************************************/
  22
  23 FinalSpectra(const Char_t *rawfilename = "TOF_rawSpectra.root", const Char_t *matchfilename = "TOF_matchingEfficiency.root", const Char_t *trackfilename = "TOF_trackingEfficiency.root", const Char_t *primaryfilename = "TOF_primaryFraction.root", const Char_t *electronfilename = "TOF_electronCorrection.root")
  24  {
  25
  26    Int_t marker[2] = {20, 25};
  27    Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
  28
  29    TFile *rawfile = TFile::Open(rawfilename);
  30    TFile *matchfile = TFile::Open(matchfilename);
  31    TFile *trackfile = TFile::Open(trackfilename);
  32    TFile *primaryfile = TFile::Open(primaryfilename);
  33    TFile *electronfile = TFile::Open(electronfilename);
  34
  35    TFile *fileout = TFile::Open("TOF_finalSpectra.root", "RECREATE");
  36    TH1D *hRaw, *hMatch, *hTrack, *hPrim, *hElectron;
  37    Char_t title[1024];
  38    for (Int_t icent = 0; icent < NcentralityBins; icent++)
  39      for (Int_t icharge = 0; icharge < kNCharges; icharge++)
  40        for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
  41          hRaw = (TH1D *)rawfile->Get(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  42          if (!hRaw) {
  43            printf("cannot find %s in %s\n", Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), rawfilename);
  44            continue;
  45          }
  46          hMatch = (TH1D *)matchfile->Get(Form("hMatchEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  47          if (!hMatch) {
  48            printf("cannot find %s in %s\n", Form("hMatchEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), matchfilename);
  49            continue;
  50          }
  51          hTrack = (TH1D *)trackfile->Get(Form("hTrackingEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  52          if (!hTrack) {
  53            printf("cannot find %s in %s\n", Form("hTrackingEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), trackfilename);
  54            continue;
  55          }
  56          hPrim = (TH1D *)primaryfile->Get(Form("hPrimaryFrac_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  57          if (ipart != 3 && !hPrim) {
  58            printf("cannot find %s in %s\n", Form("hPrimaryFrac_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), primaryfilename);
  59            continue;
  60          }
  61          hElectron = (TH1D *)electronfile->Get("hElectronCorr_average");
  62          if (!hElectron) {
  63            printf("cannot find hElectronCorr_average in %s\n", electronfilename);
  64            continue;
  65          }
  66          hRaw->Divide(hMatch);
  67          hRaw->Divide(hTrack);
  68          if (hPrim) hRaw->Multiply(hPrim);
  69          if (ipart == 2) hRaw->Multiply(hElectron);
  70          sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partLatex[ipart][icharge], centralityBin[icent], centralityBin[icent + 1]);
  71          hRaw->SetName(Form("hFinal_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  72          hRaw->SetTitle(title);
  73          hRaw->SetMarkerStyle(marker[icharge]);
  74          hRaw->SetMarkerColor(color[ipart]);
  75          fileout->cd();
  76          hRaw->Write();
  77        }
  78
  79    fileout->Close();
  80
  81    FinalRatios();
  82  }
  83
  84
  85 /**************************************************************/
  86
  87 FinalRatios(const Char_t *filename = "TOF_finalSpectra.root")
  88 {
  89
  90   printf("WARNING: skipping ratios\n");
  91   return;
  92
  93   Int_t marker[2] = {20, 25};
  94   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
  95
  96   /* open data */
  97   TFile *filein = TFile::Open(filename);
  98   /* get spectra */
  99   TH1D *hSpectrum[NcentralityBins][AliPID::kSPECIES][kNCharges];
 100   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 101     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 102       for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 103         hSpectrum[icent][ipart][icharge] = (TH1D *)filein->Get(Form("hFinal_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 104       }
 105     }
 106   }
 107
 108   /* output */
 109   TFile *fileout = TFile::Open("TOF_finalRatios.root", "RECREATE");
 110
 111   /* particle/anti-particle ratios */
 112   TH1D *hRatio;
 113   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 114     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 115       if (!hSpectrum[icent][ipart][kNegative] ||
 116           !hSpectrum[icent][ipart][kPositive]) continue;
 117       hRatio = new TH1D(*hSpectrum[icent][ipart][kNegative]);
 118       hRatio->Divide(hSpectrum[icent][ipart][kPositive]);
 119       hRatio->SetName(Form("hRatio_cent%d_%s_negative_%s_positive", icent, AliPID::ParticleName(ipart), AliPID::ParticleName(ipart)));
 120       hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partLatex[ipart][kNegative], partLatex[ipart][kPositive], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partLatex[ipart][kNegative], partLatex[ipart][kPositive]));
 121       hRatio->SetMarkerStyle(20);
 122       hRatio->SetMarkerColor(color[ipart]);
 123       fileout->cd();
 124       hRatio->Write();
 125     }
 126   }
 127
 128   /* kaon/pion ratios */
 129   TH1D *hSum1, *hSum2;
 130   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 131     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 132       if (!hSpectrum[icent][ipart][kNegative] ||
 133           !hSpectrum[icent][AliPID::kPion][icharge]) continue;
 134       hRatio = new TH1D(*hSpectrum[icent][AliPID::kKaon][icharge]);
 135       hRatio->Divide(hSpectrum[icent][AliPID::kPion][icharge]);
 136       hRatio->SetName(Form("hRatio_cent%d_kaon_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 137       hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partLatex[AliPID::kKaon][icharge], partLatex[AliPID::kPion][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partLatex[AliPID::kKaon][icharge], partLatex[AliPID::kPion][icharge]));
 138       hRatio->SetMarkerStyle(marker[icharge]);
 139       hRatio->SetMarkerColor(color[AliPID::kKaon]);
 140       fileout->cd();
 141       hRatio->Write();
 142     }
 143     if (!hSpectrum[icent][AliPID::kKaon][kPositive] ||
 144         !hSpectrum[icent][AliPID::kKaon][kNegative] ||
 145         !hSpectrum[icent][AliPID::kPion][kPositive] ||
 146         !hSpectrum[icent][AliPID::kPion][kNegative]) continue;
 147     hSum1 = new TH1D(*hSpectrum[icent][AliPID::kKaon][kPositive]);
 148     hSum1->Add(hSpectrum[icent][AliPID::kKaon][kNegative]);
 149     hSum2 = new TH1D(*hSpectrum[icent][AliPID::kPion][kPositive]);
 150     hSum2->Add(hSpectrum[icent][AliPID::kPion][kNegative]);
 151     hRatio = new TH1D(*hSum1);
 152     hRatio->Divide(hSum2);
 153     hRatio->SetName(Form("hRatio_cent%d_kaon_pion", icent));
 154     hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (%d-%d%%);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partLatex[AliPID::kKaon][kPositive], partLatex[AliPID::kKaon][kNegative], partLatex[AliPID::kPion][kPositive], partLatex[AliPID::kPion][kNegative], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partLatex[AliPID::kKaon][kPositive], partLatex[AliPID::kKaon][kNegative], partLatex[AliPID::kPion][kPositive], partLatex[AliPID::kPion][kNegative]));
 155     hRatio->SetMarkerStyle(20);
 156     hRatio->SetMarkerColor(color[AliPID::kKaon]);
 157     fileout->cd();
 158     hRatio->Write();
 159   }
 160
 161   /* proton/pion ratios */
 162   TH1D *hSum1, *hSum2;
 163   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 164     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 165       if (!hSpectrum[icent][AliPID::kProton][icharge] ||
 166           !hSpectrum[icent][AliPID::kPion][icharge]) continue;
 167       hRatio = new TH1D(*hSpectrum[icent][AliPID::kProton][icharge]);
 168       hRatio->Divide(hSpectrum[icent][AliPID::kPion][icharge]);
 169       hRatio->SetName(Form("hRatio_cent%d_proton_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 170       hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partLatex[AliPID::kProton][icharge], partLatex[AliPID::kPion][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partLatex[AliPID::kProton][icharge], partLatex[AliPID::kPion][icharge]));
 171       hRatio->SetMarkerStyle(marker[icharge]);
 172       hRatio->SetMarkerColor(color[AliPID::kProton]);
 173       fileout->cd();
 174       hRatio->Write();
 175     }
 176     if (!hSpectrum[icent][AliPID::kProton][kPositive] ||
 177         !hSpectrum[icent][AliPID::kProton][kNegative] ||
 178         !hSpectrum[icent][AliPID::kPion][kPositive] ||
 179         !hSpectrum[icent][AliPID::kPion][kNegative]) continue;
 180     hSum1 = new TH1D(*hSpectrum[icent][AliPID::kProton][kPositive]);
 181     hSum1->Add(hSpectrum[icent][AliPID::kProton][kNegative]);
 182     hSum2 = new TH1D(*hSpectrum[icent][AliPID::kPion][kPositive]);
 183     hSum2->Add(hSpectrum[icent][AliPID::kPion][kNegative]);
 184     hRatio = new TH1D(*hSum1);
 185     hRatio->Divide(hSum2);
 186     hRatio->SetName(Form("hRatio_cent%d_proton_pion", icent));
 187     hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (%d-%d%%);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partLatex[AliPID::kProton][kPositive], partLatex[AliPID::kProton][kNegative], partLatex[AliPID::kPion][kPositive], partLatex[AliPID::kPion][kNegative], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partLatex[AliPID::kProton][kPositive], partLatex[AliPID::kProton][kNegative], partLatex[AliPID::kPion][kPositive], partLatex[AliPID::kPion][kNegative]));
 188     hRatio->SetMarkerStyle(20);
 189     hRatio->SetMarkerColor(color[AliPID::kProton]);
 190     fileout->cd();
 191     hRatio->Write();
 192   }
 193
 194
 195   fileout->Close();
 196 }
 197