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

   1 #include "CommonDefs.C"
   2
   3 /**************************************************************/
   4 Char_t *partChargeName[AliPID::kSPECIES][2] = {
   5   "", "", "", "", "#pi^{+}", "#pi^{-}", "K^{+}", "K^{-}", "p", "#bar{p}"
   6 };
   7  /**************************************************************/
   8
   9 FinalSpectra()
  10 {
  11
  12   doFinalSpectra(1., "TOF_rawSpectra.root", "TOF_matchingEfficiency.root", "TOF_trackingEfficiency.root", "TOF_primaryFraction.root", "TOF_electronCorrection.root", "TOF_finalSpectra.root");
  13   doFinalSpectra(1., "TOF_rawSpectra.root", "TOF_matchingEfficiency.root", "TOF_trackingEfficiency.root", "TOF_primaryFraction.root", "TOF_electronCorrection.root", "TOF_finalSpectra_centCorr.root", kFALSE);
  14
  15   doFinalRatios("TOF_finalSpectra.root", "TOF_finalRatios.root");
  16
  17   FinalSpectra_mismatchCorrected();
  18 }
  19
  20 FinalSpectra_mismatchCorrected()
  21 {
  22   doFinalSpectra(1., "TOF_rawSpectra_mismatchCorrected.root", "TOF_matchingEfficiency_mismatchCorrected.root", "TOF_trackingEfficiency.root", "TOF_primaryFraction.root", "TOF_electronCorrection.root", "TOF_finalSpectra_mismatchCorrected.root");
  23   doFinalSpectra(1., "TOF_rawSpectra_mismatchCorrected.root", "TOF_matchingEfficiency_mismatchCorrected.root", "TOF_trackingEfficiency.root", "TOF_primaryFraction.root", "TOF_electronCorrection.root", "TOF_finalSpectra_mismatchCorrected_centCorr.root", kFALSE);
  24
  25   doFinalRatios("TOF_finalSpectra_mismatchCorrected.root", "TOF_finalRatios_mismatchCorrected.root");
  26 }
  27
  28 doFinalSpectra(Double_t scaleFact = 1., 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", const Char_t *outputfilename = "TOF_finalSpectra.root", Bool_t useMB = kTRUE)
  29  {
  30
  31    Int_t marker[2] = {20, 25};
  32    Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
  33
  34    TFile *rawfile = TFile::Open(rawfilename);
  35    TFile *matchfile = TFile::Open(matchfilename);
  36    TFile *trackfile = TFile::Open(trackfilename);
  37    TFile *primaryfile = TFile::Open(primaryfilename);
  38    TFile *electronfile = TFile::Open(electronfilename);
  39
  40    TFile *fileout = TFile::Open(outputfilename, "RECREATE");
  41    TH1D *hRaw, *hMatch, *hTrack, *hPrim, *hElectron;
  42    TH1D *hSpectrum[NcentralityBins+1][5][2];
  43    Char_t title[1024];
  44    for (Int_t icent = -1; icent < NcentralityBins; icent++)
  45      for (Int_t icharge = 0; icharge < kNCharges; icharge++)
  46        for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
  47          if (icent == -1)
  48            hRaw = (TH1D *)rawfile->Get(Form("hRaw_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
  49          else
  50            hRaw = (TH1D *)rawfile->Get(Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  51          if (!hRaw) {
  52            printf("cannot find %s in %s\n", Form("hRaw_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), rawfilename);
  53            continue;
  54          }
  55          if (0)
  56            hMatch = (TH1D *)matchfile->Get(Form("hMatchEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  57          else
  58            hMatch = (TH1D *)matchfile->Get(Form("hMatchEff_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
  59          if (!hMatch) {
  60            printf("cannot find %s in %s\n", Form("hMatchEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), matchfilename);
  61            continue;
  62          }
  63          if (0)
  64            hTrack = (TH1D *)trackfile->Get(Form("hTrackingEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  65          else
  66            hTrack = (TH1D *)trackfile->Get(Form("hTrackingEff_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
  67          if (!hTrack) {
  68            printf("cannot find %s in %s\n", Form("hTrackingEff_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), trackfilename);
  69            continue;
  70          }
  71          if (useMB || icent == -1) {
  72            hPrim = (TH1D *)primaryfile->Get(Form("hPrimaryFrac_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
  73            if (ipart != 3 && !hPrim) {
  74              printf("cannot find %s in %s\n", Form("hPrimaryFrac_MB_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), primaryfilename);
  75              continue;
  76            }
  77          }
  78          else {
  79            hPrim = (TH1D *)primaryfile->Get(Form("hPrimaryFrac_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
  80            if (ipart != 3 && !hPrim) {
  81              printf("cannot find %s in %s\n", Form("hPrimaryFrac_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]), primaryfilename);
  82              continue;
  83            }
  84          }
  85          hElectron = (TH1D *)electronfile->Get("hElectronCorr_average");
  86          if (!hElectron) {
  87            printf("cannot find hElectronCorr_average in %s\n", electronfilename);
  88            continue;
  89          }
  90          hRaw->Divide(hMatch);
  91          hRaw->Divide(hTrack);
  92          if (hPrim) hRaw->Multiply(hPrim);
  93          if (ipart == 2) hRaw->Multiply(hElectron);
  94          if (icent == -1)
  95            sprintf(title, "%s (MB);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partChargeName[ipart][icharge]);
  96          else
  97            sprintf(title, "%s (%d-%d%%);p_{T} (GeV/c);#frac{d^{2}N}{dy dp_{T}} (c/GeV);", partChargeName[ipart][icharge], centralityBin[icent], centralityBin[icent + 1]);
  98          if (icent == -1)
  99            hRaw->SetName(Form("hFinal_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 100          else
 101            hRaw->SetName(Form("hFinal_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 102          hRaw->SetTitle(title);
 103          hRaw->SetMarkerStyle(marker[icharge]);
 104          hRaw->SetMarkerColor(color[ipart]);
 105          fileout->cd();
 106          hRaw->Scale(scaleFact);
 107          hRaw->Write();
 108          if (icent == -1)
 109            hSpectrum[NcentralityBins][ipart][icharge] = hRaw;
 110          else
 111            hSpectrum[icent][ipart][icharge] = hRaw;
 112        }
 113
 114    TH1 *hRatioMB;
 115    for (Int_t icent = 0; icent < NcentralityBins; icent++)
 116      for (Int_t icharge = 0; icharge < kNCharges; icharge++)
 117        for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 118          hRatioMB = (TH1 *)hSpectrum[icent][ipart][icharge]->Clone(Form("hRatioMB_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 119          hRatioMB->Divide(hRatioMB, hSpectrum[NcentralityBins][ipart][icharge], 1., 1., "");
 120          fileout->cd();
 121          hRatioMB->Write();
 122        }
 123
 124    fileout->Close();
 125
 126  }
 127
 128
 129 /**************************************************************/
 130
 131 doFinalRatios(const Char_t *filename = "TOF_finalSpectra.root", const Char_t *outputfilename = "TOF_finalRatios.root")
 132 {
 133
 134   //  printf("WARNING: skipping ratios\n");
 135   //  return;
 136
 137   Int_t marker[2] = {20, 25};
 138   Int_t color[AliPID::kSPECIES] = {1, 1, 4, 8, 2};
 139
 140   /* open data */
 141   TFile *filein = TFile::Open(filename);
 142   /* get spectra */
 143   TH1D *hSpectrum[NcentralityBins+1][AliPID::kSPECIES][kNCharges];
 144   for (Int_t icent = 0; icent < NcentralityBins+1; icent++) {
 145     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 146       for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 147         if (icent == NcentralityBins)
 148           hSpectrum[icent][ipart][icharge] = (TH1D *)filein->Get(Form("hFinal_MB_%s_%s", AliPID::ParticleName(ipart), chargeName[icharge]));
 149         else
 150           hSpectrum[icent][ipart][icharge] = (TH1D *)filein->Get(Form("hFinal_cent%d_%s_%s", icent, AliPID::ParticleName(ipart), chargeName[icharge]));
 151       }
 152     }
 153   }
 154
 155   /* output */
 156   TFile *fileout = TFile::Open(outputfilename, "RECREATE");
 157
 158   /* particle/anti-particle ratios */
 159   TH1D *hRatio;
 160   TH1D *hPosNegRatio[NcentralityBins+1][5];
 161   for (Int_t icent = 0; icent < NcentralityBins+1; icent++) {
 162     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 163       if (!hSpectrum[icent][ipart][kNegative] ||
 164           !hSpectrum[icent][ipart][kPositive]) continue;
 165       hRatio = new TH1D(*hSpectrum[icent][ipart][kNegative]);
 166       hRatio->Divide(hSpectrum[icent][ipart][kPositive]);
 167       if (icent == NcentralityBins) {
 168         hRatio->SetName(Form("hRatio_MB_%s_negative_%s_positive", AliPID::ParticleName(ipart), AliPID::ParticleName(ipart)));
 169         hRatio->SetTitle(Form("%s/%s (MB);p_{T} (GeV/c);%s/%s;", partChargeName[ipart][kNegative], partChargeName[ipart][kPositive], partChargeName[ipart][kNegative], partChargeName[ipart][kPositive]));
 170       }
 171       else {
 172         hRatio->SetName(Form("hRatio_cent%d_%s_negative_%s_positive", icent, AliPID::ParticleName(ipart), AliPID::ParticleName(ipart)));
 173         hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partChargeName[ipart][kNegative], partChargeName[ipart][kPositive], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partChargeName[ipart][kNegative], partChargeName[ipart][kPositive]));
 174       }
 175       hRatio->SetMarkerStyle(20);
 176       hRatio->SetMarkerColor(color[ipart]);
 177       fileout->cd();
 178       hRatio->Write();
 179       hPosNegRatio[icent][ipart] = hRatio;
 180     }
 181   }
 182
 183   TH1 *hRatioMB;
 184   for (Int_t icent = 0; icent < NcentralityBins; icent++)
 185     for (Int_t ipart = 2; ipart < AliPID::kSPECIES; ipart++) {
 186       hRatioMB = (TH1 *)hPosNegRatio[icent][ipart]->Clone(Form("hDoubleRatioMB_cent%d_%s_negative_%s_positive", icent, AliPID::ParticleName(ipart), AliPID::ParticleName(ipart)));
 187       hRatioMB->Divide(hRatioMB, hPosNegRatio[NcentralityBins][ipart], 1., 1., "");
 188       fileout->cd();
 189       hRatioMB->Write();
 190     }
 191
 192   /* kaon/pion ratios */
 193   TH1D *hSum1, *hSum2;
 194   TH1D *hKaToPi[NcentralityBins+1][3];
 195   for (Int_t icent = 0; icent < NcentralityBins+1; icent++) {
 196     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 197       if (!hSpectrum[icent][AliPID::kKaon][icharge] ||
 198           !hSpectrum[icent][AliPID::kPion][icharge]) continue;
 199       hRatio = new TH1D(*hSpectrum[icent][AliPID::kKaon][icharge]);
 200       hRatio->Divide(hSpectrum[icent][AliPID::kPion][icharge]);
 201       if (icent == NcentralityBins) {
 202         hRatio->SetName(Form("hRatio_MB_kaon_%s_pion_%s", chargeName[icharge], chargeName[icharge]));
 203         hRatio->SetTitle(Form("%s/%s (MB);p_{T} (GeV/c);%s/%s;", partChargeName[AliPID::kKaon][icharge], partChargeName[AliPID::kPion][icharge], partChargeName[AliPID::kKaon][icharge], partChargeName[AliPID::kPion][icharge]));
 204       }
 205       else {
 206         hRatio->SetName(Form("hRatio_cent%d_kaon_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 207         hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partChargeName[AliPID::kKaon][icharge], partChargeName[AliPID::kPion][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partChargeName[AliPID::kKaon][icharge], partChargeName[AliPID::kPion][icharge]));
 208       }
 209       hRatio->SetMarkerStyle(marker[icharge]);
 210       hRatio->SetMarkerColor(color[AliPID::kKaon]);
 211       fileout->cd();
 212       hRatio->Write();
 213       hKaToPi[icent][icharge] = hRatio;
 214     }
 215     if (!hSpectrum[icent][AliPID::kKaon][kPositive] ||
 216         !hSpectrum[icent][AliPID::kKaon][kNegative] ||
 217         !hSpectrum[icent][AliPID::kPion][kPositive] ||
 218         !hSpectrum[icent][AliPID::kPion][kNegative]) continue;
 219     hSum1 = new TH1D(*hSpectrum[icent][AliPID::kKaon][kPositive]);
 220     hSum1->Add(hSpectrum[icent][AliPID::kKaon][kNegative]);
 221     hSum2 = new TH1D(*hSpectrum[icent][AliPID::kPion][kPositive]);
 222     hSum2->Add(hSpectrum[icent][AliPID::kPion][kNegative]);
 223     hRatio = new TH1D(*hSum1);
 224     hRatio->Divide(hSum2);
 225     if (icent == NcentralityBins) {
 226       hRatio->SetName(Form("hRatio_MB_kaon_pion"));
 227       hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (MB);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partChargeName[AliPID::kKaon][kPositive], partChargeName[AliPID::kKaon][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative], partChargeName[AliPID::kKaon][kPositive], partChargeName[AliPID::kKaon][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative]));
 228     }
 229     else {
 230       hRatio->SetName(Form("hRatio_cent%d_kaon_pion", icent));
 231       hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (%d-%d%%);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partChargeName[AliPID::kKaon][kPositive], partChargeName[AliPID::kKaon][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partChargeName[AliPID::kKaon][kPositive], partChargeName[AliPID::kKaon][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative]));
 232     }
 233     hRatio->SetMarkerStyle(20);
 234     hRatio->SetMarkerColor(color[AliPID::kKaon]);
 235     fileout->cd();
 236     hRatio->Write();
 237     hKaToPi[icent][2] = hRatio;
 238   }
 239
 240   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 241     for (Int_t icharge = 0; icharge < 2; icharge++) {
 242       hRatioMB = (TH1 *)hKaToPi[icent][icharge]->Clone(Form("hDoubleRatioMB_cent%d_kaon_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 243       hRatioMB->Divide(hRatioMB, hKaToPi[NcentralityBins][icharge], 1., 1., "");
 244       fileout->cd();
 245       hRatioMB->Write();
 246     }
 247     hRatioMB = (TH1 *)hKaToPi[icent][2]->Clone(Form("hDoubleRatioMB_cent%d_kaon_pion", icent));
 248     hRatioMB->Divide(hRatioMB, hKaToPi[NcentralityBins][2], 1., 1., "");
 249     fileout->cd();
 250     hRatioMB->Write();
 251   }
 252
 253
 254   /* proton/pion ratios */
 255   TH1D *hSum1, *hSum2;
 256   TH1D *hPrToPi[NcentralityBins+1][3];
 257   for (Int_t icent = 0; icent < NcentralityBins+1; icent++) {
 258     for (Int_t icharge = 0; icharge < kNCharges; icharge++) {
 259       if (!hSpectrum[icent][AliPID::kProton][icharge] ||
 260           !hSpectrum[icent][AliPID::kPion][icharge]) continue;
 261       hRatio = new TH1D(*hSpectrum[icent][AliPID::kProton][icharge]);
 262       hRatio->Divide(hSpectrum[icent][AliPID::kPion][icharge]);
 263       if (icent == NcentralityBins) {
 264         hRatio->SetName(Form("hRatio_MB_proton_%s_pion_%s", chargeName[icharge], chargeName[icharge]));
 265         hRatio->SetTitle(Form("%s/%s (MB);p_{T} (GeV/c);%s/%s;", partChargeName[AliPID::kProton][icharge], partChargeName[AliPID::kPion][icharge], partChargeName[AliPID::kProton][icharge], partChargeName[AliPID::kPion][icharge]));
 266       }
 267       else {
 268         hRatio->SetName(Form("hRatio_cent%d_proton_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 269         hRatio->SetTitle(Form("%s/%s (%d-%d%%);p_{T} (GeV/c);%s/%s;", partChargeName[AliPID::kProton][icharge], partChargeName[AliPID::kPion][icharge], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partChargeName[AliPID::kProton][icharge], partChargeName[AliPID::kPion][icharge]));
 270       }
 271       hRatio->SetMarkerStyle(marker[icharge]);
 272       hRatio->SetMarkerColor(color[AliPID::kProton]);
 273       fileout->cd();
 274       hRatio->Write();
 275       hPrToPi[icent][icharge] = hRatio;
 276     }
 277     if (!hSpectrum[icent][AliPID::kProton][kPositive] ||
 278         !hSpectrum[icent][AliPID::kProton][kNegative] ||
 279         !hSpectrum[icent][AliPID::kPion][kPositive] ||
 280         !hSpectrum[icent][AliPID::kPion][kNegative]) continue;
 281     hSum1 = new TH1D(*hSpectrum[icent][AliPID::kProton][kPositive]);
 282     hSum1->Add(hSpectrum[icent][AliPID::kProton][kNegative]);
 283     hSum2 = new TH1D(*hSpectrum[icent][AliPID::kPion][kPositive]);
 284     hSum2->Add(hSpectrum[icent][AliPID::kPion][kNegative]);
 285     hRatio = new TH1D(*hSum1);
 286     hRatio->Divide(hSum2);
 287     if (icent == NcentralityBins) {
 288       hRatio->SetName(Form("hRatio_MB_proton_pion"));
 289       hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (MB);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partChargeName[AliPID::kProton][kPositive], partChargeName[AliPID::kProton][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative], partChargeName[AliPID::kProton][kPositive], partChargeName[AliPID::kProton][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative]));
 290     }
 291     else {
 292       hRatio->SetName(Form("hRatio_cent%d_proton_pion", icent));
 293       hRatio->SetTitle(Form("(%s+%s)/(%s+%s) (%d-%d%%);p_{T} (GeV/c);(%s+%s)/(%s+%s);", partChargeName[AliPID::kProton][kPositive], partChargeName[AliPID::kProton][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative], (Int_t)centralityBin[icent], (Int_t)centralityBin[icent + 1], partChargeName[AliPID::kProton][kPositive], partChargeName[AliPID::kProton][kNegative], partChargeName[AliPID::kPion][kPositive], partChargeName[AliPID::kPion][kNegative]));
 294     }
 295     hRatio->SetMarkerStyle(20);
 296     hRatio->SetMarkerColor(color[AliPID::kProton]);
 297     fileout->cd();
 298     hRatio->Write();
 299     hPrToPi[icent][2] = hRatio;
 300   }
 301
 302   for (Int_t icent = 0; icent < NcentralityBins; icent++) {
 303     for (Int_t icharge = 0; icharge < 2; icharge++) {
 304       hRatioMB = (TH1 *)hPrToPi[icent][icharge]->Clone(Form("hDoubleRatioMB_cent%d_proton_%s_pion_%s", icent, chargeName[icharge], chargeName[icharge]));
 305       hRatioMB->Divide(hRatioMB, hPrToPi[NcentralityBins][icharge], 1., 1., "");
 306       fileout->cd();
 307       hRatioMB->Write();
 308     }
 309     hRatioMB = (TH1 *)hPrToPi[icent][2]->Clone(Form("hDoubleRatioMB_cent%d_proton_pion", icent));
 310     hRatioMB->Divide(hRatioMB, hPrToPi[NcentralityBins][2], 1., 1., "");
 311     fileout->cd();
 312     hRatioMB->Write();
 313   }
 314
 315
 316   fileout->Close();
 317 }
 318