PWGLF/SPECTRA/Nuclei/B2/macros/RatioMult.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // particle ratios as a function of multiplicity
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #include <TROOT.h>
  20 #include <TFile.h>
  21 #include <TString.h>
  22 #include <TGraphErrors.h>
  23 #include <TF1.h>
  24 #include <TCanvas.h>
  25 #include <TMath.h>
  26 #include <TStyle.h>
  27
  28 #include "B2.h"
  29 #include "Config.h"
  30
  31 void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr);
  32 void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle);
  33
  34 void RatioMult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
  35                const TString& ptag     = "lhc10d",
  36                const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
  37                const TString& dtag     = "lhc10d")
  38 {
  39 //
  40 // particle ratios as a function of multiplicity (from fitting distributions)
  41 //
  42         using namespace B2mult;
  43
  44         const Int_t kNpart = 2;
  45
  46         const TString kProton[kNpart] = { "Proton", "AntiProton" };
  47         const TString kDeuteron[kNpart] = { "Deuteron", "AntiDeuteron" };
  48
  49         // open files
  50
  51         TFile* fproton = new TFile(pSpectra.Data());
  52         if(fproton->IsZombie()) exit(1);
  53
  54         TFile* fdeuteron = new TFile(dSpectra.Data());
  55         if(fdeuteron->IsZombie()) exit(1);
  56
  57         // particle ratios
  58         TGraphErrors*  grRatio[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
  59         TGraphErrors*  grMixRatio[kNpart] = { new TGraphErrors(), new TGraphErrors() };
  60
  61         // model parameters
  62         TGraphErrors*  grProtondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
  63         TGraphErrors*  grProtonN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
  64         TGraphErrors*  grProtonC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
  65         TGraphErrors*  grProtonChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
  66         //TGraphErrors*  grProtonProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
  67
  68         TGraphErrors*  grDeuterondNdy[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
  69         TGraphErrors*  grDeuteronN[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
  70         TGraphErrors*  grDeuteronC[kNpart]       = { new TGraphErrors(), new TGraphErrors() };
  71         TGraphErrors*  grDeuteronChi2Ndf[kNpart] = { new TGraphErrors(), new TGraphErrors() };
  72         //TGraphErrors*  grDeuteronProb[kNpart]    = { new TGraphErrors(), new TGraphErrors() };
  73
  74         for(Int_t i=0; i<kNmult; ++i)
  75         {
  76                 TF1* fncProton[kNpart];
  77                 TF1* fncDeuteron[kNpart];
  78
  79                 for(Int_t j=0; j<kNpart; ++j)
  80                 {
  81                         fncProton[j] = (TF1*)FindObj(fproton, ptag + "-" + kMultClass[i], kProton[j] + "_Fit_InvDiffYield_Pt");
  82                         fncDeuteron[j] = (TF1*)FindObj(fdeuteron, dtag + "-" + kMultClass[i], kDeuteron[j] + "_Fit_InvDiffYield_Pt");
  83                 }
  84
  85                 // integrated ratios
  86
  87                 Double_t dNdy[kNpart] = {fncProton[0]->GetParameter(0), fncDeuteron[0]->GetParameter(0)};
  88                 Double_t dNdyErr[kNpart] = {fncProton[0]->GetParError(0), fncDeuteron[0]->GetParError(0)};
  89
  90                 Double_t dNdyBar[kNpart] = {fncProton[1]->GetParameter(0), fncDeuteron[1]->GetParameter(0)};
  91                 Double_t dNdyBarErr[kNpart] = {fncProton[1]->GetParError(0), fncDeuteron[1]->GetParError(0)};
  92
  93                 // ratios
  94                 Double_t ratio[kNpart], ratioErr[kNpart];
  95                 for(Int_t j=0; j<kNpart; ++j)
  96                 {
  97                         GetRatio(dNdyBar[j], dNdy[j], dNdyBarErr[j], dNdyErr[j], ratio[j], ratioErr[j]);
  98                         grRatio[j]->SetPoint(i, kKNOmult[i], ratio[j]);
  99                         grRatio[j]->SetPointError(i, 0, ratioErr[j]);
 100                 }
 101
 102                 // mixed ratios
 103                 Double_t mixRatio[kNpart], mixRatioErr[kNpart];
 104
 105                 GetRatio(dNdy[1], dNdy[0], dNdyErr[1], dNdyErr[0], mixRatio[0], mixRatioErr[0]);
 106                 GetRatio(dNdyBar[1], dNdyBar[0], dNdyBarErr[1], dNdyBarErr[0], mixRatio[1], mixRatioErr[1]);
 107
 108                 for(Int_t j=0; j<kNpart; ++j)
 109                 {
 110                         grMixRatio[j]->SetPoint(i, kKNOmult[i], mixRatio[j]);
 111                         grMixRatio[j]->SetPointError(i, 0, mixRatioErr[j]);
 112                 }
 113
 114                 // model parameters
 115                 for(Int_t j=0; j<kNpart; ++j)
 116                 {
 117                         grProtondNdy[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(0));
 118                         grProtondNdy[j]->SetPointError(i, 0, fncProton[j]->GetParError(0));
 119
 120                         grProtonN[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(1));
 121                         grProtonN[j]->SetPointError(i, 0, fncProton[j]->GetParError(1));
 122
 123                         grProtonC[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetParameter(2));
 124                         grProtonC[j]->SetPointError(i, 0, fncProton[j]->GetParError(2));
 125
 126                         grProtonChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetChisquare()/fncProton[j]->GetNDF());
 127                         //grProtonProb[j]->SetPoint(i, kKNOmult[i], fncProton[j]->GetProb());
 128
 129                         // deuteron
 130                         grDeuterondNdy[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(0));
 131                         grDeuterondNdy[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(0));
 132
 133                         grDeuteronN[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(1));
 134                         grDeuteronN[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(1));
 135
 136                         grDeuteronC[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetParameter(2));
 137                         grDeuteronC[j]->SetPointError(i, 0, fncDeuteron[j]->GetParError(2));
 138
 139                         grDeuteronChi2Ndf[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetChisquare()/fncDeuteron[j]->GetNDF());
 140                         //grDeuteronProb[j]->SetPoint(i, kKNOmult[i], fncDeuteron[j]->GetProb());
 141                 }
 142         }
 143
 144         delete fproton;
 145         delete fdeuteron;
 146
 147         // draw
 148
 149         TStyle* st = new TStyle();
 150
 151         st->SetPadTickX(1);
 152         st->SetPadTickY(1);
 153         st->SetPadGridX(1);
 154         st->SetPadGridY(1);
 155         st->SetCanvasColor(0);
 156         st->SetFrameBorderMode(0);
 157         st->SetFrameFillColor(0);
 158         st->SetLabelFont(62,"XYZ");
 159         st->SetTitleFont(62,"XYZ");
 160
 161         st->cd();
 162         gROOT->ForceStyle();
 163
 164         const Int_t kNCol = 4;
 165
 166         TCanvas* c0 = new TCanvas("c0", "proton model parameters");
 167         c0->Divide(kNCol,2);
 168
 169         for(Int_t j=0; j<kNpart; ++j)
 170         {
 171                 c0->cd(kNCol*j+1);
 172                 Draw(grProtondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
 173
 174                 c0->cd(kNCol*j+2);
 175                 Draw(grProtonN[j], kFullCircle, kBlue, "z", "n");
 176
 177                 c0->cd(kNCol*j+3);
 178                 Draw(grProtonC[j], kFullCircle, kBlue, "z", "C");
 179
 180                 c0->cd(kNCol*j+4);
 181                 Draw(grProtonChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
 182
 183         /*      c0->cd(kNCol*j+5);
 184                 Draw(grProtonProb[j], kFullCircle, kBlue, "z", "Prob");
 185         */
 186         }
 187
 188         TCanvas* c1 = new TCanvas("c1", "deuteron model parameters");
 189         c1->Divide(kNCol,2);
 190
 191         for(Int_t j=0; j<kNpart; ++j)
 192         {
 193                 c1->cd(kNCol*j+1);
 194                 Draw(grDeuterondNdy[j], kFullCircle, kBlue, "z", "dN/dy");
 195
 196                 c1->cd(kNCol*j+2);
 197                 Draw(grDeuteronN[j], kFullCircle, kBlue, "z", "n");
 198
 199                 c1->cd(kNCol*j+3);
 200                 Draw(grDeuteronC[j], kFullCircle, kBlue, "z", "C");
 201
 202                 c1->cd(kNCol*j+4);
 203                 Draw(grDeuteronChi2Ndf[j], kFullCircle, kBlue, "z", "#chi^{2}/ndf");
 204
 205         /*      c1->cd(kNCol*j+5);
 206                 Draw(grDeuteronProb[j], kFullCircle, kBlue, "z", "Prob");
 207         */
 208         }
 209
 210         TCanvas* c2 = new TCanvas("c2","particle ratios");
 211         c2->Divide(2,2);
 212
 213         c2->cd(1);
 214         Draw(grRatio[0], kFullCircle, kRed, "z", "#bar{p}/p");
 215
 216         c2->cd(2);
 217         Draw(grRatio[1], kFullCircle, kRed, "z", "#bar{d}/d");
 218
 219         c2->cd(3);
 220         Draw(grMixRatio[0], kFullCircle, kRed, "z", "d/p");
 221
 222         c2->cd(4);
 223         Draw(grMixRatio[1], kFullCircle, kRed, "z", "#bar{d}/#bar{p}");
 224 }
 225
 226 void GetRatio(Double_t x, Double_t y, Double_t errx, Double_t erry, Double_t& r, Double_t& rerr)
 227 {
 228 //
 229 // get the ratio of x/y, its error and save in r and rerr
 230 //
 231         if(x == 0 || y == 0)
 232         {
 233                 r=0;
 234                 rerr=1e+16;
 235
 236                 return;
 237         }
 238
 239         r = x/y;
 240         rerr = r*TMath::Sqrt(TMath::Power(errx/x,2)+TMath::Power(erry/y,2));
 241 }
 242
 243 void Draw(TGraph* gr, Int_t marker, Int_t color, const TString& xtitle, const TString& ytitle)
 244 {
 245 //
 246 // draw a graph in current canvas
 247 //
 248         gr->SetMarkerStyle(marker);
 249         gr->SetMarkerColor(color);
 250         gr->SetLineColor(color);
 251         gr->GetXaxis()->SetTitle(xtitle.Data());
 252         gr->GetYaxis()->SetTitle(ytitle.Data());
 253         gr->Draw("AP");
 254 }