PWGLF/ThermalFits/FitNPartDependence.C

   1 #if !defined (__CINT__) || (defined(__MAKECINT__))
   2 #include <iostream>
   3 #include <map>
   4 #include <fstream>
   5 #include "TString.h"
   6 #include "TObjArray.h"
   7 #include "TObjString.h"
   8 #include "AliParticleYield.h"
   9 #include "TGraphErrors.h"
  10 #include "TF1.h"
  11 #include "TROOT.h"
  12 #include "TVirtualFitter.h"
  13 #include "TMath.h"
  14 #include "TFitResult.h"
  15 #include "TMinuit.h"
  16 #include "TLatex.h"
  17 #include "TDatabasePDG.h"
  18 #include "TPad.h"
  19 #include "TCanvas.h"
  20
  21 #endif
  22
  23
  24
  25
  26 std::map<TString,Float_t> npartPbPb;
  27 std::map<TString,Float_t> npartPbPbErr;
  28
  29 void ReadCentralityFromFile() ;
  30 Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color,  Bool_t draw = kFALSE) ;
  31 Double_t ErrorFunction (Double_t *xx, Double_t *p ) ;
  32 Double_t FuncPlusErr (Double_t *xx, Double_t *p) ;
  33
  34 enum {kNoShift, kShiftUp, kShiftDown, kShiftHarder, kShiftSofter};
  35
  36
  37 TF1 * fForErrors = 0;
  38 const int npar = 3;
  39 Double_t matrix[npar][npar];
  40 TString centrFile;
  41 TString systemAndEnergy;
  42 Double_t maxy = 0;
  43 //TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (#it{N}_{part})^{#it{c}}";
  44 TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (d#it{N}/d#it{#eta})^{#it{c}}";
  45 Int_t collSystem  = 2;
  46 Float_t energy = 2760;
  47
  48 void FitNPartDependence() {
  49
  50   // This macro can be used to fit the centrality dependence of
  51   // different particles in order to extrapolate or interpolate the
  52   // yield as a function of centrality.
  53   // The value of npart or dndeta is read from a text file, which must
  54   // be in the format "centrTag value error"
  55
  56   // The particle yields are taken from the machine readable files.
  57   // If you want to fit a different particle, you have to uncomment
  58   // the corresponding section below.  if you want to try a different
  59   // function, scroll down to the definition of f1
  60
  61   // For the estimate of the systematic uncertainties, the syst graph
  62   // is shifted up and down within its uncertainties. Few strategies
  63   // are possible (see the call to FitShiftedGraphAndExtrapolate below
  64   // to change the behavior):
  65   // 1. All points are coherently shifted up and down by 1 syst bar
  66   // (kShiftUp, kShiftDown)
  67   // 2. A worst case scenario is implemented, where the first point is
  68   // shifted by +- 1 sigma, the last one by -+ 1 sigma and all the
  69   // other points are shifted smoothly in between (kShiftHarder,
  70   // kShiftSofter)
  71
  72   //__________________________________________________________________//
  73
  74   // WARNING: check isSum
  75
  76   // KStar
  77   // centrFile = "npart_PbPb.txt";
  78   // //  centrFile = "dndeta_PbPb.txt";
  79   // maxy = 50;
  80   // systemAndEnergy = "Pb-Pb #sqrt{#it{s}}_{NN} = 2.76 TeV";
  81   // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080" , 0};
  82   // const char * centrToExtrapolate = "V0M0010";
  83   // Int_t pdg = 313;
  84   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./PbPb_2760_Kstar892.txt");
  85   // Deuteron
  86   // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
  87   // const char * centralityToPlot[] = {   "V0M0010" ,  "V0M1020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080",0 };
  88   // const char * centrToExtrapolate = "V0M0010";
  89   // Int_t pdg = 1000010020;
  90   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
  91   //  Deuteron pPb
  92   centrFile = "dndeta_pPb.txt";
  93   const char * centralityToPlot[] = {   "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
  94   const char * centrToExtrapolate = "V0A0005";
  95   //const char * centrToExtrapolate = "V0A6080";
  96   Int_t pdg = -1000010020;
  97   TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
  98   maxy = 0.01;
  99   systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
 100   energy = 5020;
 101   collSystem = 1;
 102   // K* pPb
 103   // centrFile =  "dndeta_pPb.txt";
 104   // //  centrFile = "dndeta_PbPb.txt";
 105   // maxy = 1;
 106   // systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
 107   // const char * centralityToPlot[] = {   "V0A0020" ,  "V0A2040" ,  "V0A4060" ,  "V0A6080", "V0A8000" , 0};
 108   // const char * centrToExtrapolate = "V0A0005";
 109   // Int_t pdg = 313;
 110   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./pPb_5020_Kstar.txt");
 111   // energy = 5020;
 112   // collSystem = 1;
 113
 114
 115   // Helium3
 116   // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
 117   // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2080" ,0};
 118   // const char * centrToExtrapolate = "V0M0010";
 119   // Int_t pdg = 1000020030;
 120   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
 121
 122
 123
 124
 125   //  gPad->GetCanvas()->SetTitle(Form("c%s", TDatabasePDG::Instance()->GetParticle(pdg)->GetName()));
 126
 127   ReadCentralityFromFile();
 128   TGraphErrors * grSyst = new TGraphErrors;
 129   TGraphErrors * grStat = new TGraphErrors;
 130   Double_t maxx = 1.1*npartPbPb[centrToExtrapolate];
 131   //  Double_t maxx = 1.1*npartPbPb["V0A0005"];
 132   // Function
 133   // TF1 * f1 = new TF1 ("f1", "[0] + [1]*x", 0, maxx);
 134   // f1->SetParameters(1,1);
 135   TF1 * f1 = new TF1 ("f1", "[0] + [1]*x^[2]", 0, maxx);
 136   f1->SetParameters(0, 6.3266e-04, 8.99883e-01);
 137
 138
 139   gROOT->ProcessLine(Form(".x figTemplate.C(0,0.1,%f,%f)", maxx, maxy));
 140   gPad->GetCanvas()->SetTitle("cKstar");
 141
 142   //  f1->FixParameter(0, 0);
 143   Int_t icentr = 0;
 144   AliParticleYield * part = 0;
 145   while (centralityToPlot[icentr]) {
 146     part =  AliParticleYield::FindParticle (arr, pdg, collSystem, energy, centralityToPlot[icentr]);
 147     if(part) {
 148       grSyst->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
 149       grSyst->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetSystError());
 150       grStat->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
 151       grStat->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetStatError());
 152
 153       //      part->Print();
 154     }
 155     icentr++;
 156   }
 157   grStat->Draw("P");
 158   grStat->SetMarkerStyle(24);
 159   grSyst->Draw("PE2");
 160   grSyst->SetFillStyle(0);
 161   Double_t yield          = FitShiftedGraphAndExtrapolate(grStat, kNoShift  , f1, centrToExtrapolate, kRed  );
 162
 163   // We need to cache the covariance matrix and the function here, otherwise things get messed up
 164   fForErrors = (TF1*) f1->Clone();
 165   fForErrors->SetLineColor(kRed);
 166   fForErrors->SetLineStyle(0);
 167   fForErrors->Draw("same");
 168   gMinuit->mnemat(&matrix[0][0],npar);
 169   // The statistical error is computed propagating the uncertainty on the parameters to the function via the covariance matrix
 170   TF1 * fShiftedPlus   = new TF1("fShiftedPlus"  , FuncPlusErr, 0, maxx,1); fShiftedPlus->SetParameter(0,+1);
 171   TF1 * fShiftedMinus  = new TF1("fShiftedMinus" , FuncPlusErr, 0, maxx,1); fShiftedMinus->SetParameter(0,-1);
 172   fShiftedMinus->SetLineStyle(kDotted);
 173   fShiftedPlus->SetLineStyle (kDotted);
 174   fShiftedPlus ->Draw("same");
 175   fShiftedMinus->Draw("same");
 176   TF1 * fError = new TF1("fError", ErrorFunction, 0,maxx, 0);
 177
 178   // The uncertainty on the systematics is computed shifting the graph up and down + refitting
 179   // Double_t errorSystPlus  = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftUp  , f1, centrToExtrapolate, kRed)-yield);
 180   // Double_t errorSystMinus = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftDown, f1, centrToExtrapolate, kRed)-yield);
 181   Double_t errorSystPlus  = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftHarder, f1, centrToExtrapolate, kRed)-yield);
 182   Double_t errorSystMinus = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftSofter, f1, centrToExtrapolate, kRed)-yield);
 183
 184
 185   Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);
 186
 187
 188   std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")"
 189             << yield
 190             << "   +" << errorStat << " "
 191             << "   +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
 192             << "   +" << fError->Eval(npartPbPb[centrToExtrapolate])  << std::endl;
 193   std::cout << yield << "    "
 194             << errorStat << "     "
 195             << (errorSystPlus+errorSystMinus)/2 << std::endl;
 196
 197
 198   TGraphErrors * gExtrap = new TGraphErrors();
 199   gExtrap->SetMarkerStyle(20);
 200   gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
 201   gExtrap->SetPointError(0, 0, (errorSystPlus+errorSystMinus)/2);
 202   gExtrap->Draw("P");
 203
 204   TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
 205   text->SetNDC();
 206   text->Draw();
 207   TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
 208   text2->SetNDC();
 209   text2->Draw();
 210
 211   TLatex * func = new TLatex(0.53,0.23,label);
 212   func->SetNDC();
 213   func->Draw();
 214
 215
 216 }
 217
 218 Double_t ErrorFunction (Double_t *xx, Double_t *p ) {
 219
 220
 221
 222   //  Double_t x = xx[0];
 223   Double_t func[npar];
 224   // func[0] = 1;
 225   // func[1] = TMath::Power(x, fForErrors->GetParameter(2));
 226   // func[2] = fForErrors->GetParameter(1)*TMath::Power(x, fForErrors->GetParameter(2))*TMath::Log(x);
 227   //  In general, one can compute the derivative numerically using root.
 228   for(Int_t ipar = 0; ipar < npar; ipar++){
 229         func[ipar] =  fForErrors->GradientPar(ipar, xx);
 230   }
 231
 232
 233
 234   Double_t variance = 0;
 235
 236   for(Int_t ipar = 0; ipar < npar; ipar++){
 237     for(Int_t jpar = 0; jpar < npar; jpar++){
 238       variance = variance + func[ipar]*func[jpar] * matrix[ipar][jpar];
 239     }
 240
 241   }
 242   Double_t error = TMath::Sqrt(variance);
 243   return error;
 244 }
 245
 246 Double_t FuncPlusErr (Double_t *xx, Double_t *p) {
 247
 248   Double_t value = fForErrors->Eval(xx[0]);
 249   Double_t error = ErrorFunction(xx, p);
 250
 251   if(p[0] > 0) {
 252     return value + error;
 253   }
 254   return value -error;
 255 }
 256
 257
 258 Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw) {
 259   // Shift graph
 260   // 1 = up
 261   // 2 = down
 262
 263   TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
 264   Int_t npoint = grLocal->GetN() ;
 265   if(shift != kNoShift) {
 266     for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
 267       Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
 268       Double_t value    = grLocal->GetY() [ipoint];
 269       if(shift == kShiftUp  ) value += grLocal->GetEY()[ipoint];
 270       if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];
 271       if(shift == kShiftHarder) {
 272         Double_t valueShift = -1. + 2.*Float_t(ipoint)/(npoint-1);
 273         value += valueShift*grLocal->GetEY()[ipoint];
 274         std::cout << "valueShift " << valueShift << std::endl;
 275       }
 276       if(shift == kShiftSofter) {
 277         Double_t valueShift = +1. - 2.*Float_t(ipoint)/(npoint-1);
 278         value += valueShift*grLocal->GetEY()[ipoint];
 279         std::cout << "valueShift " << valueShift << std::endl;
 280       }
 281       grLocal->SetPoint     (ipoint, grLocal->GetX() [ipoint], value);
 282       grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
 283     }
 284   }
 285   grLocal->Fit(f1, "", "Q");
 286
 287
 288   Double_t yield = f1->Eval(npartPbPb[centr]);
 289   TF1 * clone = (TF1*) f1->Clone();
 290   clone->SetLineWidth(1);
 291   clone->SetLineStyle(kDashed);
 292   clone->SetLineColor(color);
 293   clone->Draw("Same");
 294   if(draw) {
 295     TVirtualPad * oldPad = gPad;
 296     new TCanvas(Form("shifted_%d", shift), Form("shifted_%d", shift));
 297     grLocal->Draw("AP");
 298     clone->Draw("same");
 299     oldPad->cd();
 300   } else {
 301     delete grLocal;
 302   }
 303   return yield;
 304
 305 }
 306
 307 void ReadCentralityFromFile() {
 308   ifstream filein (centrFile);
 309   TString line;
 310   while (line.ReadLine((filein))) {
 311     TObjArray * tokens = line.Tokenize(" ");
 312     if(tokens->GetEntries() != 3) {
 313       delete tokens;
 314       continue;
 315     }
 316     TString str        = ((TObjString*)tokens->At(0))->String();
 317     Float_t npart     = ((TObjString*)tokens->At(1))->String().Atof();
 318     Float_t npartErr  = ((TObjString*)tokens->At(2))->String().Atof();
 319     //    std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
 320     if(npart) {
 321       npartPbPb   [str] = npart;
 322       npartPbPbErr[str]  = npartErr;
 323     }
 324     delete tokens;
 325   }
 326
 327 }