#if !defined (__CINT__) || (defined(__MAKECINT__)) #include #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 () { // MF COMMON LIBRARIES SHOULD NOT BE LOADED FOR THIS MACRO TO WORK LoadStuff(1); // 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"); // AverageAndExtrapolate("pions_pos_2040"); AverageAndExtrapolate("pions_neg_2040"); AverageAndExtrapolate("kaons_pos_2040"); AverageAndExtrapolate("kaons_neg_2040"); AverageAndExtrapolate("protons_pos_2040"); AverageAndExtrapolate("protons_neg_2040"); } 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(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 == "pions_pos_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyPion][kMyPos]); arrStat->Add(hSpectraCentr_stat[4][kMyPion][kMyPos]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyPion][kMyPos]); arrSyst->Add(hSpectraCentr_sys [4][kMyPion][kMyPos]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]); } if(what == "pions_neg_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyPion][kMyNeg]); arrStat->Add(hSpectraCentr_stat[4][kMyPion][kMyNeg]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyPion][kMyNeg]); arrSyst->Add(hSpectraCentr_sys [4][kMyPion][kMyNeg]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[0]); } if(what == "kaons_pos_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyKaon][kMyPos]); arrStat->Add(hSpectraCentr_stat[4][kMyKaon][kMyPos]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyKaon][kMyPos]); arrSyst->Add(hSpectraCentr_sys [4][kMyKaon][kMyPos]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]); } if(what == "kaons_neg_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyKaon][kMyNeg]); arrStat->Add(hSpectraCentr_stat[4][kMyKaon][kMyNeg]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyKaon][kMyNeg]); arrSyst->Add(hSpectraCentr_sys [4][kMyKaon][kMyNeg]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[1]); } if(what == "protons_pos_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyProton][kMyPos]); arrStat->Add(hSpectraCentr_stat[4][kMyProton][kMyPos]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyProton][kMyPos]); arrSyst->Add(hSpectraCentr_sys [4][kMyProton][kMyPos]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]); } if(what == "protons_neg_2040"){ TObjArray * arrStat = new TObjArray(); arrStat->Add(hSpectraCentr_stat[3][kMyProton][kMyNeg]); arrStat->Add(hSpectraCentr_stat[4][kMyProton][kMyNeg]); TObjArray * arrSyst = new TObjArray(); arrSyst->Add(hSpectraCentr_sys [3][kMyProton][kMyNeg]); arrSyst->Add(hSpectraCentr_sys [4][kMyProton][kMyNeg]); Double_t weights[] = {0.5, 0.5}; TH1 * hstat = ReweightSpectra(arrStat, weights, 0, "_2040"); TH1 * hsyst = ReweightSpectra(arrSyst, weights, 1, "_2040"); hyieldmean = PrintYieldAndError(hsyst, hstat, 0, mass[2]); } 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(); }