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

#if !defined (__CINT__) || (defined(__MAKECINT__))
#include <iostream>
#include "TClonesArray.h"
#include "TH1.h"
#include "AliPWGFunc.h"
//#include "LoadSpectraPiKPPbPb.C"
#include "AliParticleYield.h"
#include "TFile.h"
#include "TDatabasePDG.h"

#endif


void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale = 1. , Int_t isLinear = 0) ;
TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc = 0, Double_t massParticle=0.139) ;
TH1* ReweightSpectra(TClonesArray * histos, Double_t * weights, Int_t isLinear = 0, TString suffix = "") ;
void LoadStuff(Bool_t loadPbPb = kTRUE) ;
void AverageAndExtrapolate(TString what);


TH1 * hLambdaStat[7];
TH1 * hLambdaSyst[7];
TH1 * hK0SStat[7]   ;
TH1 * hK0SSyst[7]   ;


void AverageAndExtrapolate () {

  LoadStuff(0); // True PbPb, False pPb

  //  PrintYieldAndError(hSpectraCentr_sys[0][kMyProton][0], hSpectraCentr_stat[0][kMyProton][0], 0, mass[2]);

  //  PrintYieldAndError(hLambdaSyst[0], hLambdaStat[0], 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
  // icentr, ipart, icharge

  //  AverageAndExtrapolate("pions_pos_0020");
  //  AverageAndExtrapolate("pions_pos_0010");
  // AverageAndExtrapolate("pions_neg_0010");
  AverageAndExtrapolate("pions_sum_0010");
  // AverageAndExtrapolate("kaons_pos_0010");
  // AverageAndExtrapolate("kaons_neg_0010");
  // AverageAndExtrapolate("protons_pos_0010");
  // AverageAndExtrapolate("protons_neg_0010");
  //  AverageAndExtrapolate("lambda_0010");
  //  AverageAndExtrapolate("k0s_0010");
  // AverageAndExtrapolate("pions_pos_6080");
  // AverageAndExtrapolate("pions_neg_6080");
  // AverageAndExtrapolate("kaons_pos_6080");
  // AverageAndExtrapolate("kaons_neg_6080");
  // AverageAndExtrapolate("protons_pos_6080");
  // AverageAndExtrapolate("protons_neg_6080");
  
}

TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc, Double_t massParticle) {

  TF1 * f = 0;
  if (ifunc == 0) {
    f= BGBlastWave("fBlastWave",  massParticle);
    f->SetParameters(massParticle, 0.640, 0.097, 0.73, 100); // FIXME
  }
  else if (ifunc == 1) {
    f  = fm.GetTsallis(massParticle, 0.2, 11, 800, "fTsallisPion");
  }
  
  //  TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100, 0.01, 0.1, "");
  TH1 * h =  YieldMean(hstat, hsyst, f, 0.0, 100);
  //  std::cout << "H " << h << std::endl;
  std::cout << "" << std::endl;
  std::cout << h->GetBinContent(1) << ", " << h->GetBinContent(2) << ", " << h->GetBinContent(3) << std::endl;
  std::cout << "" << std::endl;
  return h;


}


TH1* ReweightSpectra(TObjArray * histos, Double_t * weights, Int_t isLinear, TString suffix) {

  // sums a number of spectra with a given weight. Used to combine
  // several centrality bins.  The weights should add up to 1 and be
  // proportional to the width of the centrality bin for each
  // histogram
  // Suffix is added to the name of the histogram.
  // if linear = 1 errors are summed linearly rather than in quadrature


  TIter it(histos);
  TH1 * h = 0;
  TH1 * hsum = 0;
  Int_t ihisto = 0;
  while ((h = dynamic_cast<TH1*>(it.Next()))) {
    if(!h) {
      std::cout << "ERROR cannot get one of the histos!" << std::endl;
      return 0;
    }

    if (!hsum) {
      // First histogram, clone it
      hsum =  (TH1D*) h->Clone(TString(h->GetName())+suffix);
      hsum->Scale(weights[ihisto++]);
    }else {
      AddHistograms(hsum, h, weights[ihisto++], isLinear);
    }
  }
  return hsum;
}

void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale, Int_t isLinear) {
  // THis method assumes that the 2 histos are consistent!!
  TH1 * hsourceLoc = (TH1*) hsource->Clone("hsourceLoc");
  hsourceLoc->Scale(scale);
  if(!isLinear) {
    hdest->Add(hsourceLoc);
  } 
  else {
    Int_t nbin = hsourceLoc->GetNbinsX();
    for(Int_t ibin = 0; ibin < nbin; ibin++){
      Double_t content = hdest->GetBinContent(ibin) + hsourceLoc->GetBinContent(ibin);
      Double_t error   = hdest->GetBinError  (ibin) + hsourceLoc->GetBinError  (ibin);
      hdest->SetBinContent(ibin,content);
      hdest->SetBinError(ibin, error);
    }
    

  }

  delete hsourceLoc;
  
  
}


void LoadStuff(Bool_t loadPbPb) {
  if(loadPbPb){
    gROOT->LoadMacro("LoadSpectraPiKPPbPb.C");
    LoadSpectraPiKPPbPb();
  }
  else {
    gROOT->LoadMacro("LoadSpectraPiKPProtonLead.C");
    LoadSpectraPiKPProtonLead();
  }
  LoadLibs();
  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/YieldMean.C");
  gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/SpectraUtils.C");

  // Load Lambdas and K0s
  TFile * f = new TFile("k0s_lambda_final_spectra.root");
  const char * multTags[] = {  "0005",  "0510",  "1020",  "2040",  "4060",  "6080", "8090"};
  for (Int_t icentr = 0; icentr<7; icentr++) {
    hLambdaStat[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
    hLambdaSyst[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
    hK0SStat[icentr]    = (TH1*) f->Get(Form("systonly_cent%s_K0s",multTags[icentr]));
    hK0SSyst[icentr]    = (TH1*) f->Get(Form("statonly_cent%s_K0s",multTags[icentr]));

    // The bin 300-400 MeV was not used in the analysis
    hK0SStat[icentr]->SetBinContent(4,0);   
    hK0SSyst[icentr]->SetBinContent(4,0);
    hK0SStat[icentr]->SetBinError(4,0);   
    hK0SSyst[icentr]->SetBinError(4,0);
  }

}

void AverageAndExtrapolate(TString what) {

  TH1 * hstat =0;
  TH1 * hsyst =0;
  TH1 * hyieldmean =0;

  if(what == "pions_pos_0020"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[2][kMyPion][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [2][kMyPion][kMyPos]);
    
    Double_t weights[] = {0.25, 0.25, 0.5};
    //Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0020");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0020");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "pions_pos_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "pions_sum_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5, 0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "pions_neg_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "protons_pos_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);

    
  }
  if(what == "protons_neg_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);

    
  }
    if(what == "kaons_pos_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);

    
  }
  if(what == "kaons_neg_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);

    
  }

  if(what == "pions_pos_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "pions_neg_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);

    
  }
  if(what == "protons_pos_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);

    
  }
  if(what == "protons_neg_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);

    
  }
    if(what == "kaons_pos_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyPos]);
    arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyPos]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyPos]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyPos]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);

    
  }
  if(what == "kaons_neg_6080"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyNeg]);
    arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyNeg]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyNeg]);
    arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyNeg]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);

    
  }


  if(what == "lambda_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hLambdaStat[0]);
    arrStat->Add(hLambdaStat[1]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hLambdaSyst[0]);
    arrSyst->Add(hLambdaSyst[1]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());

    
  }

  if(what == "k0s_0010"){ 

    TObjArray * arrStat = new TObjArray();
    arrStat->Add(hK0SStat[0]);
    arrStat->Add(hK0SStat[1]);
    
    TObjArray * arrSyst = new TObjArray();
    arrSyst->Add(hK0SSyst[0]);
    arrSyst->Add(hK0SSyst[1]);
    
    Double_t weights[] = {0.5, 0.5};
    
    TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
    TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");

    hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("K_S0")->Mass());

    
  }


  TCanvas * c1 = new TCanvas("Averaging", "Averaging");
  c1->Divide(2,1);
  c1->cd(1);
  hstat->Draw();
  hsyst->Draw("same,e3");
  c1->cd(2);
  hyieldmean->Draw();
  

}
Commit	Line	Data
9d84731e	1	#if !defined (__CINT__) \|\| (defined(__MAKECINT__))
	2	#include <iostream>
	3	#include "TClonesArray.h"
	4	#include "TH1.h"
	5	#include "AliPWGFunc.h"
	6	//#include "LoadSpectraPiKPPbPb.C"
	7	#include "AliParticleYield.h"
	8	#include "TFile.h"
	9	#include "TDatabasePDG.h"
	10
	11	#endif
	12
	13
	14	void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale = 1. , Int_t isLinear = 0) ;
	15	TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc = 0, Double_t massParticle=0.139) ;
	16	TH1* ReweightSpectra(TClonesArray * histos, Double_t * weights, Int_t isLinear = 0, TString suffix = "") ;
4e849078	17	void LoadStuff(Bool_t loadPbPb = kTRUE) ;
9d84731e	18	void AverageAndExtrapolate(TString what);
	19
	20
	21
	22	TH1 * hLambdaStat[7];
	23	TH1 * hLambdaSyst[7];
	24	TH1 * hK0SStat[7] ;
	25	TH1 * hK0SSyst[7] ;
	26
	27
	28
	29	void AverageAndExtrapolate () {
	30
4e849078	31	LoadStuff(0); // True PbPb, False pPb
9d84731e	32
	33	// PrintYieldAndError(hSpectraCentr_sys[0][kMyProton][0], hSpectraCentr_stat[0][kMyProton][0], 0, mass[2]);
	34
	35	// PrintYieldAndError(hLambdaSyst[0], hLambdaStat[0], 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
	36	// icentr, ipart, icharge
	37
	38	// AverageAndExtrapolate("pions_pos_0020");
4e849078	39	// AverageAndExtrapolate("pions_pos_0010");
9d84731e	40	// AverageAndExtrapolate("pions_neg_0010");
4e849078	41	AverageAndExtrapolate("pions_sum_0010");
9d84731e	42	// AverageAndExtrapolate("kaons_pos_0010");
	43	// AverageAndExtrapolate("kaons_neg_0010");
	44	// AverageAndExtrapolate("protons_pos_0010");
	45	// AverageAndExtrapolate("protons_neg_0010");
	46	// AverageAndExtrapolate("lambda_0010");
	47	// AverageAndExtrapolate("k0s_0010");
ad431d97	48	// AverageAndExtrapolate("pions_pos_6080");
	49	// AverageAndExtrapolate("pions_neg_6080");
	50	// AverageAndExtrapolate("kaons_pos_6080");
	51	// AverageAndExtrapolate("kaons_neg_6080");
	52	// AverageAndExtrapolate("protons_pos_6080");
	53	// AverageAndExtrapolate("protons_neg_6080");
4e849078	54
9d84731e	55	}
	56
	57	TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc, Double_t massParticle) {
	58
	59	TF1 * f = 0;
	60	if (ifunc == 0) {
	61	f= BGBlastWave("fBlastWave", massParticle);
4e849078	62	f->SetParameters(massParticle, 0.640, 0.097, 0.73, 100); // FIXME
9d84731e	63	}
	64	else if (ifunc == 1) {
	65	f = fm.GetTsallis(massParticle, 0.2, 11, 800, "fTsallisPion");
	66	}
	67
4e849078	68	// TH1 * h = YieldMean(hstat, hsyst, f, 0.0, 100, 0.01, 0.1, "");
4e849078	69	TH1 * h = YieldMean(hstat, hsyst, f, 0.0, 100);
9d84731e	70	// std::cout << "H " << h << std::endl;
	71	std::cout << "" << std::endl;
	72	std::cout << h->GetBinContent(1) << ", " << h->GetBinContent(2) << ", " << h->GetBinContent(3) << std::endl;
	73	std::cout << "" << std::endl;
	74	return h;
	75
	76
	77	}
	78
	79
	80	TH1* ReweightSpectra(TObjArray * histos, Double_t * weights, Int_t isLinear, TString suffix) {
	81
	82	// sums a number of spectra with a given weight. Used to combine
	83	// several centrality bins. The weights should add up to 1 and be
	84	// proportional to the width of the centrality bin for each
	85	// histogram
	86	// Suffix is added to the name of the histogram.
	87	// if linear = 1 errors are summed linearly rather than in quadrature
	88
	89
	90	TIter it(histos);
	91	TH1 * h = 0;
	92	TH1 * hsum = 0;
	93	Int_t ihisto = 0;
	94	while ((h = dynamic_cast<TH1*>(it.Next()))) {
	95	if(!h) {
	96	std::cout << "ERROR cannot get one of the histos!" << std::endl;
	97	return 0;
	98	}
	99
	100	if (!hsum) {
	101	// First histogram, clone it
	102	hsum = (TH1D*) h->Clone(TString(h->GetName())+suffix);
	103	hsum->Scale(weights[ihisto++]);
	104	}else {
	105	AddHistograms(hsum, h, weights[ihisto++], isLinear);
	106	}
	107	}
	108	return hsum;
	109	}
	110
	111	void AddHistograms (TH1 * hdest, TH1 * hsource, Float_t scale, Int_t isLinear) {
	112	// THis method assumes that the 2 histos are consistent!!
	113	TH1 * hsourceLoc = (TH1*) hsource->Clone("hsourceLoc");
	114	hsourceLoc->Scale(scale);
	115	if(!isLinear) {
	116	hdest->Add(hsourceLoc);
	117	}
	118	else {
	119	Int_t nbin = hsourceLoc->GetNbinsX();
	120	for(Int_t ibin = 0; ibin < nbin; ibin++){
	121	Double_t content = hdest->GetBinContent(ibin) + hsourceLoc->GetBinContent(ibin);
	122	Double_t error = hdest->GetBinError (ibin) + hsourceLoc->GetBinError (ibin);
	123	hdest->SetBinContent(ibin,content);
	124	hdest->SetBinError(ibin, error);
	125	}
	126
	127
	128	}
	129
	130	delete hsourceLoc;
	131
	132
	133
134	}
135
136
4e849078	137	void LoadStuff(Bool_t loadPbPb) {
	138	if(loadPbPb){
	139	gROOT->LoadMacro("LoadSpectraPiKPPbPb.C");
	140	LoadSpectraPiKPPbPb();
	141	}
	142	else {
	143	gROOT->LoadMacro("LoadSpectraPiKPProtonLead.C");
	144	LoadSpectraPiKPProtonLead();
	145	}
9d84731e	146	LoadLibs();
	147	gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/YieldMean.C");
	148	gROOT->LoadMacro("$ALICE_ROOT/PWGLF/SPECTRA/UTILS/SpectraUtils.C");
	149
	150	// Load Lambdas and K0s
	151	TFile * f = new TFile("k0s_lambda_final_spectra.root");
	152	const char * multTags[] = { "0005", "0510", "1020", "2040", "4060", "6080", "8090"};
	153	for (Int_t icentr = 0; icentr<7; icentr++) {
	154	hLambdaStat[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
	155	hLambdaSyst[icentr] = (TH1*) f->Get(Form("statonly_cent%s_Lambda",multTags[icentr]));
	156	hK0SStat[icentr] = (TH1*) f->Get(Form("systonly_cent%s_K0s",multTags[icentr]));
	157	hK0SSyst[icentr] = (TH1*) f->Get(Form("statonly_cent%s_K0s",multTags[icentr]));
	158
	159	// The bin 300-400 MeV was not used in the analysis
	160	hK0SStat[icentr]->SetBinContent(4,0);
	161	hK0SSyst[icentr]->SetBinContent(4,0);
	162	hK0SStat[icentr]->SetBinError(4,0);
	163	hK0SSyst[icentr]->SetBinError(4,0);
	164	}
	165
	166	}
	167
	168	void AverageAndExtrapolate(TString what) {
	169
	170	TH1 * hstat =0;
	171	TH1 * hsyst =0;
	172	TH1 * hyieldmean =0;
	173
	174	if(what == "pions_pos_0020"){
	175
	176	TObjArray * arrStat = new TObjArray();
	177	arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
	178	arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
	179	arrStat->Add(hSpectraCentr_stat[2][kMyPion][kMyPos]);
	180
	181	TObjArray * arrSyst = new TObjArray();
	182	arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
	183	arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
	184	arrSyst->Add(hSpectraCentr_sys [2][kMyPion][kMyPos]);
	185
	186	Double_t weights[] = {0.25, 0.25, 0.5};
	187	//Double_t weights[] = {0.5, 0.5};
	188
	189	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0020");
	190	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0020");
	191
	192	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	193
	194
	195	}
	196	if(what == "pions_pos_0010"){
	197
	198	TObjArray * arrStat = new TObjArray();
	199	arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
	200	arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
	201
	202	TObjArray * arrSyst = new TObjArray();
	203	arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
	204	arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
	205
	206	Double_t weights[] = {0.5, 0.5};
	207
	208	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
4e849078	209	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	210
	211	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	212
	213
	214	}
	215	if(what == "pions_sum_0010"){
	216
	217	TObjArray * arrStat = new TObjArray();
	218	arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyPos]);
	219	arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyPos]);
	220	arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
	221	arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
	222
	223	TObjArray * arrSyst = new TObjArray();
	224	arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyPos]);
	225	arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyPos]);
	226	arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
	227	arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
	228
	229	Double_t weights[] = {0.5, 0.5, 0.5, 0.5};
	230
	231	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
9d84731e	232	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	233
	234	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	235
	236
	237	}
	238	if(what == "pions_neg_0010"){
	239
	240	TObjArray * arrStat = new TObjArray();
	241	arrStat->Add(hSpectraCentr_stat[0][kMyPion][kMyNeg]);
	242	arrStat->Add(hSpectraCentr_stat[1][kMyPion][kMyNeg]);
	243
	244	TObjArray * arrSyst = new TObjArray();
	245	arrSyst->Add(hSpectraCentr_sys [0][kMyPion][kMyNeg]);
	246	arrSyst->Add(hSpectraCentr_sys [1][kMyPion][kMyNeg]);
	247
	248	Double_t weights[] = {0.5, 0.5};
	249
	250	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
	251	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	252
	253	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	254
	255
	256	}
	257	if(what == "protons_pos_0010"){
	258
	259	TObjArray * arrStat = new TObjArray();
	260	arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyPos]);
	261	arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyPos]);
	262
	263	TObjArray * arrSyst = new TObjArray();
	264	arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyPos]);
	265	arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyPos]);
	266
	267	Double_t weights[] = {0.5, 0.5};
	268
	269	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
	270	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	271
	272	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
	273
	274
	275	}
	276	if(what == "protons_neg_0010"){
	277
	278	TObjArray * arrStat = new TObjArray();
	279	arrStat->Add(hSpectraCentr_stat[0][kMyProton][kMyNeg]);
	280	arrStat->Add(hSpectraCentr_stat[1][kMyProton][kMyNeg]);
	281
	282	TObjArray * arrSyst = new TObjArray();
	283	arrSyst->Add(hSpectraCentr_sys [0][kMyProton][kMyNeg]);
	284	arrSyst->Add(hSpectraCentr_sys [1][kMyProton][kMyNeg]);
	285
	286	Double_t weights[] = {0.5, 0.5};
	287
	288	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
	289	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	290
	291	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
	292
	293
	294	}
	295	if(what == "kaons_pos_0010"){
296
297	TObjArray * arrStat = new TObjArray();
298	arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyPos]);
299	arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyPos]);
300
301	TObjArray * arrSyst = new TObjArray();
302	arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyPos]);
303	arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyPos]);
304
305	Double_t weights[] = {0.5, 0.5};
306
307	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
308	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
309
310	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
311
312
313	}
314	if(what == "kaons_neg_0010"){
315
316	TObjArray * arrStat = new TObjArray();
317	arrStat->Add(hSpectraCentr_stat[0][kMyKaon][kMyNeg]);
318	arrStat->Add(hSpectraCentr_stat[1][kMyKaon][kMyNeg]);
319
320	TObjArray * arrSyst = new TObjArray();
321	arrSyst->Add(hSpectraCentr_sys [0][kMyKaon][kMyNeg]);
322	arrSyst->Add(hSpectraCentr_sys [1][kMyKaon][kMyNeg]);
323
324	Double_t weights[] = {0.5, 0.5};
325
326	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
327	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
328
329	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
330
331
332	}
ad431d97	333
	334	if(what == "pions_pos_6080"){
	335
	336	TObjArray * arrStat = new TObjArray();
	337	arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyPos]);
	338	arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyPos]);
	339
	340	TObjArray * arrSyst = new TObjArray();
	341	arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyPos]);
	342	arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyPos]);
	343
	344	Double_t weights[] = {0.5, 0.5};
	345
	346	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
	347	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
	348
	349	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	350
	351
	352	}
	353	if(what == "pions_neg_6080"){
	354
	355	TObjArray * arrStat = new TObjArray();
	356	arrStat->Add(hSpectraCentr_stat[7][kMyPion][kMyNeg]);
	357	arrStat->Add(hSpectraCentr_stat[8][kMyPion][kMyNeg]);
	358
	359	TObjArray * arrSyst = new TObjArray();
	360	arrSyst->Add(hSpectraCentr_sys [7][kMyPion][kMyNeg]);
	361	arrSyst->Add(hSpectraCentr_sys [8][kMyPion][kMyNeg]);
	362
	363	Double_t weights[] = {0.5, 0.5};
	364
	365	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
	366	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
	367
	368	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]);
	369
	370
	371	}
	372	if(what == "protons_pos_6080"){
	373
	374	TObjArray * arrStat = new TObjArray();
	375	arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyPos]);
	376	arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyPos]);
	377
	378	TObjArray * arrSyst = new TObjArray();
	379	arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyPos]);
	380	arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyPos]);
	381
	382	Double_t weights[] = {0.5, 0.5};
	383
	384	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
	385	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
	386
	387	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
	388
	389
	390	}
	391	if(what == "protons_neg_6080"){
	392
	393	TObjArray * arrStat = new TObjArray();
	394	arrStat->Add(hSpectraCentr_stat[7][kMyProton][kMyNeg]);
	395	arrStat->Add(hSpectraCentr_stat[8][kMyProton][kMyNeg]);
	396
397	TObjArray * arrSyst = new TObjArray();
398	arrSyst->Add(hSpectraCentr_sys [7][kMyProton][kMyNeg]);
399	arrSyst->Add(hSpectraCentr_sys [8][kMyProton][kMyNeg]);
400
401	Double_t weights[] = {0.5, 0.5};
402
403	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
404	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
405
406	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]);
407
408
409	}
410	if(what == "kaons_pos_6080"){
411
412	TObjArray * arrStat = new TObjArray();
413	arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyPos]);
414	arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyPos]);
415
416	TObjArray * arrSyst = new TObjArray();
417	arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyPos]);
418	arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyPos]);
419
420	Double_t weights[] = {0.5, 0.5};
421
422	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
423	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
424
425	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
426
427
428	}
429	if(what == "kaons_neg_6080"){
430
431	TObjArray * arrStat = new TObjArray();
432	arrStat->Add(hSpectraCentr_stat[7][kMyKaon][kMyNeg]);
433	arrStat->Add(hSpectraCentr_stat[8][kMyKaon][kMyNeg]);
434
435	TObjArray * arrSyst = new TObjArray();
436	arrSyst->Add(hSpectraCentr_sys [7][kMyKaon][kMyNeg]);
437	arrSyst->Add(hSpectraCentr_sys [8][kMyKaon][kMyNeg]);
438
439	Double_t weights[] = {0.5, 0.5};
440
441	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_6080");
442	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_6080");
443
444	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]);
445
446
447	}
448
449
9d84731e	450	if(what == "lambda_0010"){
	451
	452	TObjArray * arrStat = new TObjArray();
	453	arrStat->Add(hLambdaStat[0]);
	454	arrStat->Add(hLambdaStat[1]);
	455
	456	TObjArray * arrSyst = new TObjArray();
	457	arrSyst->Add(hLambdaSyst[0]);
	458	arrSyst->Add(hLambdaSyst[1]);
	459
	460	Double_t weights[] = {0.5, 0.5};
	461
	462	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
	463	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	464
	465	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("Lambda0")->Mass());
	466
	467
	468	}
	469
	470	if(what == "k0s_0010"){
	471
	472	TObjArray * arrStat = new TObjArray();
	473	arrStat->Add(hK0SStat[0]);
	474	arrStat->Add(hK0SStat[1]);
	475
	476	TObjArray * arrSyst = new TObjArray();
	477	arrSyst->Add(hK0SSyst[0]);
	478	arrSyst->Add(hK0SSyst[1]);
	479
	480	Double_t weights[] = {0.5, 0.5};
	481
	482	TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_0010");
	483	TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_0010");
	484
	485	hyieldmean = PrintYieldAndError(hsyst, hstat, 0, TDatabasePDG::Instance()->GetParticle("K_S0")->Mass());
	486
	487
	488	}
	489
	490
	491
	492
	493	TCanvas * c1 = new TCanvas("Averaging", "Averaging");
	494	c1->Divide(2,1);
	495	c1->cd(1);
	496	hstat->Draw();
	497	hsyst->Draw("same,e3");
	498	c1->cd(2);
	499	hyieldmean->Draw();
	500
	501
	502	}
	503
	504