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   // Double_t errorStatPlus  = FitShiftedGraphAndExtrapolate(grStat, kShiftUp  , f1, centrToExtrapolate, kBlue) -yield;
 185   // Double_t errorStatMinus = FitShiftedGraphAndExtrapolate(grStat, kShiftDown, f1, centrToExtrapolate, kBlue) -yield;
 186
 187   Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);
 188
 189
 190   std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")"
 191             << yield
 192             << "   +" << errorStat << " "
 193             << "   +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
 194             << "   +" << fError->Eval(npartPbPb[centrToExtrapolate])  << std::endl;
 195   std::cout << yield << "    "
 196             << errorStat << "     "
 197             << (errorSystPlus+errorSystMinus)/2 << std::endl;
 198
 199
 200   TGraphErrors * gExtrap = new TGraphErrors();
 201   gExtrap->SetMarkerStyle(20);
 202   gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
 203   gExtrap->SetPointError(0, 0, (errorSystPlus+errorSystMinus)/2);
 204   gExtrap->Draw("P");
 205
 206   TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
 207   text->SetNDC();
 208   text->Draw();
 209   TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
 210   text2->SetNDC();
 211   text2->Draw();
 212
 213   TLatex * func = new TLatex(0.53,0.23,label);
 214   func->SetNDC();
 215   func->Draw();
 216
 217
 218 }
 219
 220 Double_t ErrorFunction (Double_t *xx, Double_t *p ) {
 221
 222
 223
 224   //  Double_t x = xx[0];
 225   Double_t func[npar];
 226   // func[0] = 1;
 227   // func[1] = TMath::Power(x, fForErrors->GetParameter(2));
 228   // func[2] = fForErrors->GetParameter(1)*TMath::Power(x, fForErrors->GetParameter(2))*TMath::Log(x);
 229   //  In general, one can compute the derivative numerically using root.
 230   for(Int_t ipar = 0; ipar < npar; ipar++){
 231         func[ipar] =  fForErrors->GradientPar(ipar, xx);
 232   }
 233
 234
 235
 236   Double_t variance = 0;
 237
 238   for(Int_t ipar = 0; ipar < npar; ipar++){
 239     for(Int_t jpar = 0; jpar < npar; jpar++){
 240       variance = variance + func[ipar]*func[jpar] * matrix[ipar][jpar];
 241     }
 242
 243   }
 244   Double_t error = TMath::Sqrt(variance);
 245   return error;
 246 }
 247
 248 Double_t FuncPlusErr (Double_t *xx, Double_t *p) {
 249
 250   Double_t value = fForErrors->Eval(xx[0]);
 251   Double_t error = ErrorFunction(xx, p);
 252
 253   if(p[0] > 0) {
 254     return value + error;
 255   }
 256   return value -error;
 257 }
 258
 259
 260 Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw) {
 261   // Shift graph
 262   // 1 = up
 263   // 2 = down
 264
 265   TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
 266   Int_t npoint = grLocal->GetN() ;
 267   if(shift != kNoShift) {
 268     for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
 269       Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
 270       Double_t value    = grLocal->GetY() [ipoint];
 271       if(shift == kShiftUp  ) value += grLocal->GetEY()[ipoint];
 272       if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];
 273       if(shift == kShiftHarder) {
 274         Double_t valueShift = -1. + 2.*Float_t(ipoint)/(npoint-1);
 275         value += valueShift*grLocal->GetEY()[ipoint];
 276         std::cout << "valueShift " << valueShift << std::endl;
 277       }
 278       if(shift == kShiftSofter) {
 279         Double_t valueShift = +1. - 2.*Float_t(ipoint)/(npoint-1);
 280         value += valueShift*grLocal->GetEY()[ipoint];
 281         std::cout << "valueShift " << valueShift << std::endl;
 282       }
 283       grLocal->SetPoint     (ipoint, grLocal->GetX() [ipoint], value);
 284       grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
 285     }
 286   }
 287   grLocal->Fit(f1, "", "Q");
 288
 289
 290   Double_t yield = f1->Eval(npartPbPb[centr]);
 291   TF1 * clone = (TF1*) f1->Clone();
 292   clone->SetLineWidth(1);
 293   clone->SetLineStyle(kDashed);
 294   clone->SetLineColor(color);
 295   clone->Draw("Same");
 296   if(draw) {
 297     TVirtualPad * oldPad = gPad;
 298     new TCanvas(Form("shifted_%d", shift), Form("shifted_%d", shift));
 299     grLocal->Draw("AP");
 300     clone->Draw("same");
 301     oldPad->cd();
 302   } else {
 303     delete grLocal;
 304   }
 305   return yield;
 306
 307 }
 308
 309 void ReadCentralityFromFile() {
 310   ifstream filein (centrFile);
 311   TString line;
 312   while (line.ReadLine((filein))) {
 313     TObjArray * tokens = line.Tokenize(" ");
 314     if(tokens->GetEntries() != 3) {
 315       delete tokens;
 316       continue;
 317     }
 318     TString str        = ((TObjString*)tokens->At(0))->String();
 319     Float_t npart     = ((TObjString*)tokens->At(1))->String().Atof();
 320     Float_t npartErr  = ((TObjString*)tokens->At(2))->String().Atof();
 321     //    std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
 322     if(npart) {
 323       npartPbPb   [str] = npart;
 324       npartPbPbErr[str]  = npartErr;
 325     }
 326     delete tokens;
 327   }
 328
 329 }