[u/mrichter/AliRoot.git] / PWGLF / ThermalFits / FitNPartDependence.C

#if !defined (__CINT__) || (defined(__MAKECINT__))
#include <iostream>
#include <map>
#include <fstream>
#include "TString.h"
#include "TObjArray.h"
#include "TObjString.h"
#include "AliParticleYield.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TROOT.h"
#include "TVirtualFitter.h"
#include "TMath.h"
#include "TFitResult.h"
#include "TMinuit.h"
#include "TLatex.h"
#include "TDatabasePDG.h"
#include "TPad.h"
#include "TCanvas.h"

#endif


std::map<TString,Float_t> npartPbPb;
std::map<TString,Float_t> npartPbPbErr;

void ReadCentralityFromFile() ;
Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color,  Bool_t draw = kFALSE) ;
Double_t ErrorFunction (Double_t *xx, Double_t *p ) ;
Double_t FuncPlusErr (Double_t *xx, Double_t *p) ;

enum {kNoShift, kShiftUp, kShiftDown, kShiftHarder, kShiftSofter};


TF1 * fForErrors = 0;
const int npar = 3;
Double_t matrix[npar][npar];
TString centrFile;
TString systemAndEnergy;
Double_t maxy = 0;
//TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (#it{N}_{part})^{#it{c}}";
TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (d#it{N}/d#it{#eta})^{#it{c}}";
Int_t collSystem  = 2;
Float_t energy = 2760;

void FitNPartDependence() {

  // This macro can be used to fit the centrality dependence of
  // different particles in order to extrapolate or interpolate the
  // yield as a function of centrality.
  // The value of npart or dndeta is read from a text file, which must
  // be in the format "centrTag value error"

  // The particle yields are taken from the machine readable files.
  // If you want to fit a different particle, you have to uncomment
  // the corresponding section below.  if you want to try a different
  // function, scroll down to the definition of f1

  // For the estimate of the systematic uncertainties, the syst graph
  // is shifted up and down within its uncertainties. Few strategies
  // are possible (see the call to FitShiftedGraphAndExtrapolate below
  // to change the behavior):
  // 1. All points are coherently shifted up and down by 1 syst bar
  // (kShiftUp, kShiftDown)
  // 2. A worst case scenario is implemented, where the first point is
  // shifted by +- 1 sigma, the last one by -+ 1 sigma and all the
  // other points are shifted smoothly in between (kShiftHarder,
  // kShiftSofter)

  //__________________________________________________________________//

  // WARNING: check isSum

  // KStar
  // centrFile = "npart_PbPb.txt";
  // //  centrFile = "dndeta_PbPb.txt";
  // maxy = 50;
  // systemAndEnergy = "Pb-Pb #sqrt{#it{s}}_{NN} = 2.76 TeV";
  // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080" , 0};
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 313;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./PbPb_2760_Kstar892.txt");  
  // Deuteron
  // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");  
  // const char * centralityToPlot[] = {   "V0M0010" ,  "V0M1020" ,  "V0M2040" ,  "V0M4060" ,  "V0M6080",0 };
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 1000010020;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
  //  Deuteron pPb
  centrFile = "dndeta_pPb.txt";
  const char * centralityToPlot[] = {   "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
  const char * centrToExtrapolate = "V0A0005";
  //const char * centrToExtrapolate = "V0A6080";
  Int_t pdg = -1000010020;
  TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt");
  maxy = 0.01;
  systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
  energy = 5020;
  collSystem = 1;
  // K* pPb
  // centrFile =  "dndeta_pPb.txt";
  // //  centrFile = "dndeta_PbPb.txt";
  // maxy = 1;
  // systemAndEnergy = "p-Pb #sqrt{#it{s}}_{NN} = 5.02 TeV";
  // const char * centralityToPlot[] = {   "V0A0020" ,  "V0A2040" ,  "V0A4060" ,  "V0A6080", "V0A8000" , 0};
  // const char * centrToExtrapolate = "V0A0005";
  // Int_t pdg = 313;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("./pPb_5020_Kstar.txt");  
  // energy = 5020;
  // collSystem = 1;


  // Helium3
  // gROOT->ProcessLine(".x figTemplate.C(0,0.00001,400,0.2)");  
  // const char * centralityToPlot[] = {   "V0M0020" ,  "V0M2080" ,0};
  // const char * centrToExtrapolate = "V0M0010";
  // Int_t pdg = 1000020030;
  // TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");


  //  gPad->GetCanvas()->SetTitle(Form("c%s", TDatabasePDG::Instance()->GetParticle(pdg)->GetName()));

  ReadCentralityFromFile();
  TGraphErrors * grSyst = new TGraphErrors; 
  TGraphErrors * grStat = new TGraphErrors; 
  Double_t maxx = 1.1*npartPbPb[centrToExtrapolate];
  //  Double_t maxx = 1.1*npartPbPb["V0A0005"];
  // Function
  // TF1 * f1 = new TF1 ("f1", "[0] + [1]*x", 0, maxx);
  // f1->SetParameters(1,1);
  TF1 * f1 = new TF1 ("f1", "[0] + [1]*x^[2]", 0, maxx);
  f1->SetParameters(0, 6.3266e-04, 8.99883e-01);

  
  gROOT->ProcessLine(Form(".x figTemplate.C(0,0.1,%f,%f)", maxx, maxy));
  gPad->GetCanvas()->SetTitle("cKstar");

  //  f1->FixParameter(0, 0);
  Int_t icentr = 0;
  AliParticleYield * part = 0;
  while (centralityToPlot[icentr]) {
    part =  AliParticleYield::FindParticle (arr, pdg, collSystem, energy, centralityToPlot[icentr]);
    if(part) {
      grSyst->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
      grSyst->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetSystError());
      grStat->SetPoint     (icentr , npartPbPb[centralityToPlot[icentr]]    , part->GetYield());
      grStat->SetPointError(icentr , npartPbPbErr[centralityToPlot[icentr]] , part->GetStatError());
      
      //      part->Print();
    }
    icentr++;
  }
  grStat->Draw("P");
  grStat->SetMarkerStyle(24);
  grSyst->Draw("PE2");
  grSyst->SetFillStyle(0);
  Double_t yield          = FitShiftedGraphAndExtrapolate(grStat, kNoShift  , f1, centrToExtrapolate, kRed  );

  // We need to cache the covariance matrix and the function here, otherwise things get messed up
  fForErrors = (TF1*) f1->Clone();
  fForErrors->SetLineColor(kRed);
  fForErrors->SetLineStyle(0);
  fForErrors->Draw("same");
  gMinuit->mnemat(&matrix[0][0],npar);
  // The statistical error is computed propagating the uncertainty on the parameters to the function via the covariance matrix
  TF1 * fShiftedPlus   = new TF1("fShiftedPlus"  , FuncPlusErr, 0, maxx,1); fShiftedPlus->SetParameter(0,+1);
  TF1 * fShiftedMinus  = new TF1("fShiftedMinus" , FuncPlusErr, 0, maxx,1); fShiftedMinus->SetParameter(0,-1);
  fShiftedMinus->SetLineStyle(kDotted);
  fShiftedPlus->SetLineStyle (kDotted);
  fShiftedPlus ->Draw("same");
  fShiftedMinus->Draw("same");
  TF1 * fError = new TF1("fError", ErrorFunction, 0,maxx, 0);

  // The uncertainty on the systematics is computed shifting the graph up and down + refitting
  // Double_t errorSystPlus  = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftUp  , f1, centrToExtrapolate, kRed)-yield);
  // Double_t errorSystMinus = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftDown, f1, centrToExtrapolate, kRed)-yield);
  Double_t errorSystPlus  = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftHarder, f1, centrToExtrapolate, kRed)-yield);
  Double_t errorSystMinus = TMath::Abs(FitShiftedGraphAndExtrapolate(grSyst, kShiftSofter, f1, centrToExtrapolate, kRed)-yield);


  Double_t errorStat = fError->Eval(npartPbPb[centrToExtrapolate]);


  std::cout << "Yield from fit: ("<<centrToExtrapolate<< "="<<npartPbPb[centrToExtrapolate]<<")" 
            << yield 
            << "   +" << errorStat << " " 
            << "   +" << errorSystPlus << " " << errorSystMinus//<< std::endl;
            << "   +" << fError->Eval(npartPbPb[centrToExtrapolate])  << std::endl;
  std::cout << yield << "    " 
            << errorStat << "     " 
            << (errorSystPlus+errorSystMinus)/2 << std::endl;
  

  TGraphErrors * gExtrap = new TGraphErrors();
  gExtrap->SetMarkerStyle(20);
  gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
  gExtrap->SetPointError(0, 0, (errorSystPlus+errorSystMinus)/2);
  gExtrap->Draw("P");

  TLatex * text = new TLatex (0.2,0.81,systemAndEnergy);
  text->SetNDC();
  text->Draw();
  TLatex * text2 = new TLatex (0.2,0.72,part->GetLatexName());
  text2->SetNDC();
  text2->Draw(); 

  TLatex * func = new TLatex(0.53,0.23,label);
  func->SetNDC();
  func->Draw();


}

Double_t ErrorFunction (Double_t *xx, Double_t *p ) {


  //  Double_t x = xx[0];
  Double_t func[npar];
  // func[0] = 1;
  // func[1] = TMath::Power(x, fForErrors->GetParameter(2));
  // func[2] = fForErrors->GetParameter(1)*TMath::Power(x, fForErrors->GetParameter(2))*TMath::Log(x);
  //  In general, one can compute the derivative numerically using root.
  for(Int_t ipar = 0; ipar < npar; ipar++){
        func[ipar] =  fForErrors->GradientPar(ipar, xx);
  }

  
  Double_t variance = 0;
  
  for(Int_t ipar = 0; ipar < npar; ipar++){
    for(Int_t jpar = 0; jpar < npar; jpar++){
      variance = variance + func[ipar]*func[jpar] * matrix[ipar][jpar];
    }
    
  }
  Double_t error = TMath::Sqrt(variance);
  return error;
}

Double_t FuncPlusErr (Double_t *xx, Double_t *p) {

  Double_t value = fForErrors->Eval(xx[0]);
  Double_t error = ErrorFunction(xx, p);

  if(p[0] > 0) {
    return value + error;
  }
  return value -error;
}


Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw) {
  // Shift graph
  // 1 = up
  // 2 = down

  TGraphErrors * grLocal = (TGraphErrors*) gr->Clone();
  Int_t npoint = grLocal->GetN() ;
  if(shift != kNoShift) {
    for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
      Double_t relError = grLocal->GetEY()[ipoint]/grLocal->GetY()[ipoint];
      Double_t value    = grLocal->GetY() [ipoint];
      if(shift == kShiftUp  ) value += grLocal->GetEY()[ipoint];        
      if(shift == kShiftDown) value -= grLocal->GetEY()[ipoint];   
      if(shift == kShiftHarder) {        
        Double_t valueShift = -1. + 2.*Float_t(ipoint)/(npoint-1);
        value += valueShift*grLocal->GetEY()[ipoint];
        std::cout << "valueShift " << valueShift << std::endl;
      }
      if(shift == kShiftSofter) {        
        Double_t valueShift = +1. - 2.*Float_t(ipoint)/(npoint-1);
        value += valueShift*grLocal->GetEY()[ipoint];
        std::cout << "valueShift " << valueShift << std::endl;
      }
      grLocal->SetPoint     (ipoint, grLocal->GetX() [ipoint], value);
      grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
    }
  }
  grLocal->Fit(f1, "", "Q");


  Double_t yield = f1->Eval(npartPbPb[centr]);
  TF1 * clone = (TF1*) f1->Clone();
  clone->SetLineWidth(1);
  clone->SetLineStyle(kDashed);
  clone->SetLineColor(color);
  clone->Draw("Same");
  if(draw) {
    TVirtualPad * oldPad = gPad;
    new TCanvas(Form("shifted_%d", shift), Form("shifted_%d", shift));
    grLocal->Draw("AP");
    clone->Draw("same");
    oldPad->cd();
  } else {
    delete grLocal;
  }
  return yield;

}

void ReadCentralityFromFile() {
  ifstream filein (centrFile);
  TString line;
  while (line.ReadLine((filein))) {    
    TObjArray * tokens = line.Tokenize(" ");
    if(tokens->GetEntries() != 3) {
      delete tokens;
      continue;
    }
    TString str        = ((TObjString*)tokens->At(0))->String();
    Float_t npart     = ((TObjString*)tokens->At(1))->String().Atof();
    Float_t npartErr  = ((TObjString*)tokens->At(2))->String().Atof();
    //    std::cout << "["<<str.Data() <<"]" << npart << " " << npartErr << std::endl;
    if(npart) {
      npartPbPb   [str] = npart;
      npartPbPbErr[str]  = npartErr;
    }
    delete tokens;
  }

}
Commit	Line	Data
368b4921	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() ;
6a734c5e	30	Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw = kFALSE) ;
368b4921	31	Double_t ErrorFunction (Double_t xx, Double_t p ) ;
	32	Double_t FuncPlusErr (Double_t xx, Double_t p) ;
	33
6a734c5e	34	enum {kNoShift, kShiftUp, kShiftDown, kShiftHarder, kShiftSofter};
368b4921	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;
4c7bd9d1	43	//TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (#it{N}_{part})^{#it{c}}";
4c7bd9d1	44	TString label = "#frac{d#it{N}}{d#it{y}} = #it{a} + #it{b} #times (d#it{N}/d#it{#eta})^{#it{c}}";
368b4921	45	Int_t collSystem = 2;
	46	Float_t energy = 2760;
	47
	48	void FitNPartDependence() {
	49
6a734c5e	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
4c7bd9d1	74	// WARNING: check isSum
4c7bd9d1	75
368b4921	76	// KStar
53369598	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");
368b4921	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");
0bf9d3e4	91	// Deuteron pPb
	92	centrFile = "dndeta_pPb.txt";
	93	const char * centralityToPlot[] = { "V0A0010", "V0A1020", "V0A2040", "V0A4060", "V0A6000" ,0};
53369598	94	const char * centrToExtrapolate = "V0A0005";
0bf9d3e4	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";
53369598	100	energy = 5020;
53369598	101	collSystem = 1;
0bf9d3e4	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;
53369598	113
53369598	114
368b4921	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];
4c7bd9d1	131	// Double_t maxx = 1.1*npartPbPb["V0A0005"];
368b4921	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]) {
4c7bd9d1	146	part = AliParticleYield::FindParticle (arr, pdg, collSystem, energy, centralityToPlot[icentr]);
368b4921	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
53369598	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);
368b4921	183
368b4921	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 << " "
6a734c5e	195	<< (errorSystPlus+errorSystMinus)/2 << std::endl;
368b4921	196
	197
	198	TGraphErrors * gExtrap = new TGraphErrors();
	199	gExtrap->SetMarkerStyle(20);
	200	gExtrap->SetPoint(0, npartPbPb[centrToExtrapolate], yield);
6a734c5e	201	gExtrap->SetPointError(0, 0, (errorSystPlus+errorSystMinus)/2);
368b4921	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
6a734c5e	258	Double_t FitShiftedGraphAndExtrapolate(TGraphErrors * gr, Int_t shift, TF1 * f1, const char * centr, Color_t color, Bool_t draw) {
368b4921	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];
6a734c5e	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	}
368b4921	281	grLocal->SetPoint (ipoint, grLocal->GetX() [ipoint], value);
	282	grLocal->SetPointError(ipoint, grLocal->GetEX()[ipoint], relError*value);
	283	}
	284	}
368b4921	285	grLocal->Fit(f1, "", "Q");
368b4921	286
6a734c5e	287
368b4921	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");
6a734c5e	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	}
368b4921	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	}