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   // KStar
  75   centrFile = "npart_PbPb.txt";
  76   //  centrFile = "dndeta_PbPb.txt";
  77   maxy = 50;
  78   systemAndEnergy = "Pb-Pb #sqrt{#it{s}}_{NN} = 2.76 TeV";
  79   const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080" , 0};
  80   const char * centrToExtrapolate = "V0M0010";
  81   Int_t pdg = 313;
  82   TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./PbPb_2760_Kstar892.txt");
  83   // Deuteron
  84   // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
  85   // const char * centralityToPlot[] = {   "V0M0010" ,  "V0M1020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080",0 };
  86   // const char * centrToExtrapolate = "V0M0010";
  87   // Int_t pdg = 1000010020;
  88   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
  89   // Deuteron pPb
  90   //  centrFile = "dndeta_pPb.txt";
  91   // const char * centralityToPlot[] = {   "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
  92   // //  const char * centrToExtrapolate = "V0A0005";
  93   // const char * centrToExtrapolate = "V0A6080";
  94   // Int_t pdg = 1000010020;
  95   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
  96   // maxy = 0.01;
  97   // systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
  98   // energy = 5020;
  99   // collSystem = 1;
 100   // Helium3
 101   // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");
 102   // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2080" ,0};
 103   // const char * centrToExtrapolate = "V0M0010";
 104   // Int_t pdg = 1000020030;
 105   // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
 106
 107
 108
 109
 110   //  gPad->GetCanvas()->SetTitle(Form("c%s", TDatabasePDG::Instance()->GetParticle(pdg)->GetName()));
 111
 112   ReadCentralityFromFile();
 113   TGraphErrors * grSyst = new TGraphErrors;
 114   TGraphErrors * grStat = new TGraphErrors;
 115   Double_t maxx = 1.1*npartPbPb[centrToExtrapolate];
 116   //Double_t maxx = 1.1*npartPbPb["V0A0005"];
 117   // Function
 118   // TF1 * f1 = new TF1 ("f1", "[0] + [1]*x", 0, maxx);
 119   // f1->SetParameters(1,1);
 120   TF1 * f1 = new TF1 ("f1", "[0] + [1]*x^[2]", 0, maxx);
 121   f1->SetParameters(0, 6.3266e-04, 8.99883e-01);
 122
 123
 124   gROOT->ProcessLine(Form(".x figTemplate.C(0,0.1,%f,%f)", maxx, maxy));
 125   gPad->GetCanvas()->SetTitle("cKstar");
 126
 127   //  f1->FixParameter(0, 0);
 128   Int_t icentr = 0;
 129   AliParticleYield * part = 0;
 130   while (centralityToPlot[icentr]) {
 131     part =  AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centralityToPlot[icentr]);
 132     if(part) {
 133       grSyst->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
 134       grSyst->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetSystError());
 135       grStat->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
 136       grStat->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetStatError());
 137
 138       //      part->Print();
 139     }
 140     icentr++;
 141   }
 142   grStat->Draw("P");
 143   grStat->SetMarkerStyle(24);
 144   grSyst->Draw("PE2");
 145   grSyst->SetFillStyle(0);
 146   Double_t yield          = FitShiftedGraphAndExtrapolate(grStat, kNoShift  , f1, centrToExtrapolate, kRed  );
 147
 148   // We need to cache the covariance matrix and the function here, otherwise things get messed up
 149   fForErrors = (TF1*) f1->Clone();
 150   fForErrors->SetLineColor(kRed);
 151   fForErrors->SetLineStyle(0);
 152   fForErrors->Draw("same");
 153   gMinuit->mnemat(&matrix[0][0],npar);
 154   // The statistical error is computed propagating the uncertainty on the parameters to the function via the covariance matrix
 155   TF1 * fShiftedPlus   = new TF1("fShiftedPlus"  , FuncPlusErr, 0, maxx,1); fShiftedPlus->SetParameter(0,+1);
 156   TF1 * fShiftedMinus  = new TF1("fShiftedMinus" , FuncPlusErr, 0, maxx,1); fShiftedMinus->SetParameter(0,-1);
 157   fShiftedMinus->SetLineStyle(kDotted);
 158   fShiftedPlus->SetLineStyle (kDotted);
 159   fShiftedPlus ->Draw("same");
 160   fShiftedMinus->Draw("same");
 161   TF1 * fError = new TF1("fError", ErrorFunction, 0,maxx, 0);
 162
 163   // The uncertainty on the systematics is computed shifting the graph up and down + refitting
 164   // Double_t errorSystPlus  = FitShiftedGraphAndExtrapolate(grSyst, kShiftUp  , f1, centrToExtrapolate, kRed)-yield;
 165   // Double_t errorSystMinus = FitShiftedGraphAndExtrapolate(grSyst, kShiftDown, f1, centrToExtrapolate, kRed)-yield;
 166   Double_t errorSystPlus  = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftHarder, f1, centrToExtrapolate, kRed)-yield);
 167   Double_t errorSystMinus = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftSofter, f1, centrToExtrapolate, kRed)-yield);
 168
 169   // Double_t errorStatPlus  = FitShiftedGraphAndExtrapolate(grStat, kShiftUp  , f1, centrToExtrapolate, kBlue) -yield;
 170   // Double_t errorStatMinus = FitShiftedGraphAndExtrapolate(grStat, kShiftDown, f1, centrToExtrapolate, kBlue) -yield;
 171
 172   Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);
 173
 174
 175   std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")"
 176             << yield
 177             << "   +" << errorStat << " "
 178             << "   +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
 179             << "   +" << fError->Eval(npartPbPb[centrToExtrapolate])  << std::endl;
 180   std::cout << yield << "    "
 181             << errorStat << "     "
 182             << (errorSystPlus+errorSystMinus)/2 << std::endl;
 183
 184
 185   TGraphErrors * gExtrap = new TGraphErrors();
 186   gExtrap->SetMarkerStyle(20);
 187   gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
 188   gExtrap->SetPointError(0, 0, (errorSystPlus+errorSystMinus)/2);
 189   gExtrap->Draw("P");
 190
 191   TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
 192   text->SetNDC();
 193   text->Draw();
 194   TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
 195   text2->SetNDC();
 196   text2->Draw();
 197
 198   TLatex * func = new TLatex(0.53,0.23,label);
 199   func->SetNDC();
 200   func->Draw();
 201
 202
 203 }
 204
 205 Double_t ErrorFunction (Double_t *xx, Double_t *p ) {
 206
 207
 208
 209   //  Double_t x = xx[0];
 210   Double_t func[npar];
 211   // func[0] = 1;
 212   // func[1] = TMath::Power(x, fForErrors->GetParameter(2));
 213   // func[2] = fForErrors->GetParameter(1)*TMath::Power(x, fForErrors->GetParameter(2))*TMath::Log(x);
 214   //  In general, one can compute the derivative numerically using root.
 215   for(Int_t ipar = 0; ipar < npar; ipar++){
 216         func[ipar] =  fForErrors->GradientPar(ipar, xx);
 217   }
 218
 219
 220
 221   Double_t variance = 0;
 222
 223   for(Int_t ipar = 0; ipar < npar; ipar++){
 224     for(Int_t jpar = 0; jpar < npar; jpar++){
 225       variance = variance + func[ipar]*func[jpar] * matrix[ipar][jpar];
 226     }
 227
 228   }
 229   Double_t error = TMath::Sqrt(variance);
 230   return error;
 231 }
 232
 233 Double_t FuncPlusErr (Double_t *xx, Double_t *p) {
 234
 235   Double_t value = fForErrors->Eval(xx[0]);
 236   Double_t error = ErrorFunction(xx, p);
 237
 238   if(p[0] > 0) {
 239     return value + error;
 240   }
 241   return value -error;
 242 }
 243
 244
 245 Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw) {
 246   // Shift graph
 247   // 1 = up
 248   // 2 = down
 249
 250   TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
 251   Int_t npoint = grLocal->GetN() ;
 252   if(shift != kNoShift) {
 253     for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
 254       Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
 255       Double_t value    = grLocal->GetY() [ipoint];
 256       if(shift == kShiftUp  ) value += grLocal->GetEY()[ipoint];
 257       if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];
 258       if(shift == kShiftHarder) {
 259         Double_t valueShift = -1. + 2.*Float_t(ipoint)/(npoint-1);
 260         value += valueShift*grLocal->GetEY()[ipoint];
 261         std::cout << "valueShift " << valueShift << std::endl;
 262       }
 263       if(shift == kShiftSofter) {
 264         Double_t valueShift = +1. - 2.*Float_t(ipoint)/(npoint-1);
 265         value += valueShift*grLocal->GetEY()[ipoint];
 266         std::cout << "valueShift " << valueShift << std::endl;
 267       }
 268       grLocal->SetPoint     (ipoint, grLocal->GetX() [ipoint], value);
 269       grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
 270     }
 271   }
 272   grLocal->Fit(f1, "", "Q");
 273
 274
 275   Double_t yield = f1->Eval(npartPbPb[centr]);
 276   TF1 * clone = (TF1*) f1->Clone();
 277   clone->SetLineWidth(1);
 278   clone->SetLineStyle(kDashed);
 279   clone->SetLineColor(color);
 280   clone->Draw("Same");
 281   if(draw) {
 282     TVirtualPad * oldPad = gPad;
 283     new TCanvas(Form("shifted_%d", shift), Form("shifted_%d", shift));
 284     grLocal->Draw("AP");
 285     clone->Draw("same");
 286     oldPad->cd();
 287   } else {
 288     delete grLocal;
 289   }
 290   return yield;
 291
 292 }
 293
 294 void ReadCentralityFromFile() {
 295   ifstream filein (centrFile);
 296   TString line;
 297   while (line.ReadLine((filein))) {
 298     TObjArray * tokens = line.Tokenize(" ");
 299     if(tokens->GetEntries() != 3) {
 300       delete tokens;
 301       continue;
 302     }
 303     TString str        = ((TObjString*)tokens->At(0))->String();
 304     Float_t npart     = ((TObjString*)tokens->At(1))->String().Atof();
 305     Float_t npartErr  = ((TObjString*)tokens->At(2))->String().Atof();
 306     //    std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
 307     if(npart) {
 308       npartPbPb   [str] = npart;
 309       npartPbPbErr[str]  = npartErr;
 310     }
 311     delete tokens;
 312   }
 313
 314 }