--- /dev/null
+#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() ;
+void AverageAndExtrapolate(TString what);
+
+
+
+TH1 * hLambdaStat[7];
+TH1 * hLambdaSyst[7];
+TH1 * hK0SStat[7] ;
+TH1 * hK0SSyst[7] ;
+
+
+
+void AverageAndExtrapolate () {
+
+ LoadStuff();
+
+ // 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("kaons_pos_0010");
+ // AverageAndExtrapolate("kaons_neg_0010");
+ // AverageAndExtrapolate("protons_pos_0010");
+ // AverageAndExtrapolate("protons_neg_0010");
+ // AverageAndExtrapolate("lambda_0010");
+ // AverageAndExtrapolate("k0s_0010");
+
+}
+
+TH1* PrintYieldAndError(TH1 * hsyst, TH1* hstat, Int_t ifunc, Double_t massParticle) {
+
+ TF1 * f = 0;
+ if (ifunc == 0) {
+ f= BGBlastWave("fBlastWave", massParticle);
+ }
+ else if (ifunc == 1) {
+ f = fm.GetTsallis(massParticle, 0.2, 11, 800, "fTsallisPion");
+ }
+
+ TH1 * h = YieldMean(hstat, hsyst, f);
+ // 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() {
+ gROOT->LoadMacro("LoadSpectraPiKPPbPb.C");
+ LoadSpectraPiKPPbPb();
+ 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_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 == "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();
+
+
+}
+
+
--- /dev/null
+#if !defined (__CINT__) || (defined(__MAKECINT__))
+#include <iostream>
+#include "TClonesArray.h"
+#include "AliParticleYield.h"
+
+#endif
+TClonesArray * arr =0;
+
+void InterpolateRatios(Int_t pdg1, Int_t pdg2) ;
+void Interpolate0010(Int_t pdg) ;
+void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest);
+
+
+void InterpolateRatiosAndYields() {
+#if !(!defined (__CINT__) || (defined(__MAKECINT__)))
+ LoadLibs();
+#endif
+ // *************** pi, K, pr *****************
+ // arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt");
+ // Interpolate0010(211);
+ // Interpolate0010(-211);
+ // Interpolate0010(321);
+ // Interpolate0010(-321);
+ // Interpolate0010(2212);
+ // Interpolate0010(-2212);
+ // InterpolateRatios(2212,211);
+ // InterpolateRatios(321,211);
+ // *************** Lambda and K0 *****************
+ // arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt");
+ // Interpolate0010(3122);
+ // Interpolate0010(310);
+ // *************** Helium 3 *****************
+ // arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt");
+ // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
+ // ExtrapolateWithConstantRatioToPions(1000020030, "V0M0020", "V0M0010");
+ // *************** Kstar *****************
+ arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt");
+ arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt"));
+ ExtrapolateWithConstantRatioToPions(313, "V0M0020", "V0M0010");
+
+
+}
+
+void DUMP(){
+ AliParticleYield * p0 = AliParticleYield::FindParticle(arr, 211, 2, 2760., "V0M0005");
+ // p0 = AliParticleYield::Add(p0, AliParticleYield::FindParticle(arr, -211, 2, 2760., "V0M0005"));
+ AliParticleYield * p2 = AliParticleYield::FindParticle(arr, 2212, 2, 2760., "V0M0005");
+ // p2 = AliParticleYield::Add(p2,AliParticleYield::FindParticle(arr, -2212, 2, 2760., "V0M0005"));
+ AliParticleYield *pratio = AliParticleYield::Divide(p2,p0);
+ pratio->Print();
+ AliParticleYield::FindRatio(arr, 2212, 211, 2, 2760, "V0M0005")->Print();
+
+}
+
+void Interpolate0010(Int_t pdg) {
+
+ AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, 2, 2760., "V0M0005");
+ AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, 2, 2760., "V0M0510");
+ p0->Scale(0.5);
+ p1->Scale(0.5);
+ AliParticleYield * pav = AliParticleYield::Add(p0,p1);
+ std::cout << pav->GetYield() << ", " << pav->GetStatError() << ", "<< pav->GetSystError() << std::endl;
+
+
+}
+
+void InterpolateRatios(Int_t pdg1, Int_t pdg2) {
+
+ AliParticleYield * ratio[2];
+ ratio[0] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., "V0M0005", 1);
+ ratio[1] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., "V0M0510", 1);
+ AliParticleYield * average = AliParticleYield::Add(ratio[0], ratio[1]);
+ average->Scale(0.5);
+ AliParticleYield * pi[2];
+ pi[0] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., "V0M0005", 0);
+ pi[0] = AliParticleYield::Add(pi[0],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., "V0M0005", 0));
+ pi[1] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., "V0M0510", 0);
+ pi[1] = AliParticleYield::Add(pi[1],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., "V0M0510", 0));
+
+
+ // Scale to get protons with errors corresponding to the ratio
+ ratio[0]->Scale(pi[0]->GetYield()) ;
+ ratio[1]->Scale(pi[1]->GetYield()) ;
+
+ ratio[0]->Add(ratio[0], ratio[1]);
+ pi[0]->Add(pi[0],pi[1]);
+ pi[0]->SetNormError(0);
+ pi[0]->SetStatError(0);
+ pi[0]->SetSystError(0);
+
+ ratio[0]->Scale(1./pi[0]->GetYield());
+ ratio[0]->SetCentr("V0M0010");
+
+ ratio[0]->Print();
+ // average->Dump();
+ average->Print();
+
+
+}
+
+void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest) {
+
+ AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg, 2, 2760., centrOrigin);
+ AliParticleYield * pionOrigin = AliParticleYield::FindParticle(arr, 211, 2, 2760., centrOrigin);
+ AliParticleYield * pionDest = AliParticleYield::FindParticle(arr, 211, 2, 2760., centrDest);
+ if(!part || !pionOrigin | !pionDest) {
+ return;
+ }
+ // The pion yield is takes
+ part->Scale (1./pionOrigin->GetYield());
+ part->Scale (pionDest->GetYield());
+ part->SetCentr(centrDest);
+ part->SetMeasurementType(part->GetMeasurementType() | AliParticleYield::kTypeExtrPionRatio);
+ part->Print();
+ // std::cout << "1" << std::endl;
+ TClonesArray * arrOut = new TClonesArray("AliParticleYield");
+ // std::cout << "2" << std::endl;
+ new((*arrOut)[0]) AliParticleYield(*part) ;
+
+ // std::cout << "3" << std::endl;
+ // std::cout << "4" << std::endl;
+ AliParticleYield::SaveAsASCIIFile(arrOut, "temp.txt");
+}
--- /dev/null
+#if !defined (__CINT__) || (defined(__MAKECINT__))
+#include <iostream>
+#include "TClonesArray.h"
+#include "AliParticleYield.h"
+#include "TH1F.h"
+#include "TCanvas.h"
+#include <fstream>
+#endif
+
+enum MyParticles { kPDGPi = 211, kPDGK = 321, kPDGProton = 2212, kPDGKS0 = 310, kPDGLambda=3122, kPDGXi=3312,kPDGOmega=3334,kPDGPhi=333,kPDGKStar=313,kPDGDeuteron=1000010020,kPDGHE3 = 1000020030, kPDGHyperTriton = 1010010030, kPDGSigmaStarPlus=3224,kPDGSigmaStarMinus=3114,kPDGXiStar=3324};
+
+TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) ;
+void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
+
+// Plots ratios for QM and saves input files for thermal models
+
+TClonesArray * PlotRatiosForQM14() {
+ LoadLibs();
+ TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
+ arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
+
+ TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
+
+ // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
+ // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
+
+ TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
+ c1->SetLogy();
+ // GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010")->Draw();
+ GetHistoRatios(arr, AliParticleYield::kCSpp, 7000, ".*")->Draw("");
+ return arr;
+}
+
+TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) {
+
+ const Int_t nratio = 10;
+ Int_t num [nratio] = {kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar};
+ Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGProton , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK};
+ Int_t isSum[nratio] = {1 ,1 ,0 ,1 ,1 ,0 ,0 ,1 ,0 ,1};
+ TH1F * h = new TH1F("hRatio", "hRatio", nratio, 1, nratio+1);
+
+ // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
+ for(Int_t iratio = 1; iratio <= nratio; iratio++){
+ AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
+ std::cout << num[iratio-1] << " " << denum[iratio-1]<< " " ;
+ if(ratio)ratio->Print("short");
+
+
+ if(!ratio) {
+ // If the ratio is not found, try to build it!
+ AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality, isSum[iratio-1]);
+ if(!part1) {// Try with the !sum, if part 1 is not found
+ part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
+ }
+ AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
+ if(!part2) {// Try with the !sum, if part 2 is not found
+ part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
+ }
+ ratio = AliParticleYield::Divide(part1, part2);
+ if(ratio) {
+ std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
+ ratio->Print("short");
+ }
+ }
+ if(ratio){
+ h->SetBinContent(iratio, ratio->GetYield());
+ h->SetBinError (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
+ h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",ratio->GetCentr().Data() , ratio->GetLatexName()));
+ }
+ else {
+ h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
+
+ }
+ }
+
+ return h;
+
+
+
+}
+
+void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges) {
+ // If "Separate charges" is true, tries to dump both charges are dumped
+ TClonesArray * arrOut = new TClonesArray("AliParticleYield");
+ const Int_t npart = 12;
+ Int_t particles [npart] = {kPDGPi ,kPDGK ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3 };
+ Int_t isSum[npart] = {1 ,1 ,0 , 1 , 0 , 1 ,0 ,1 ,1 ,0 , 1 , 0 };
+
+ Int_t ipartOut = 0; // Index for the array
+ for(Int_t ipart = 0; ipart < npart; ipart++){
+ if(!separateCharges) {
+ AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, isSum[ipart]);
+ if(!part && isSum[ipart]) {
+ //Could not find the particle, but the sum was requested: build the sum!
+ part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
+ AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
+ if(part2 && part) part = AliParticleYield::Add(part, part2);
+ else part = 0;
+ }
+ if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
+ }
+ else {
+ // ignore isSum and try to find both particles
+ Bool_t notFound = 0;
+ AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
+ if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
+ else notFound=1;
+ // Try to find antiparticle (-pdg code)
+ part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
+ if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
+ else if (notFound) {
+ // If neither charge was found, check if we at least have the sum
+ part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality, 1);
+ if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
+ }
+
+ }
+ }
+ std::cout << "Particles for thermal model fits:" << std::endl;
+ arrOut->Print("short");
+ std::cout << "" << std::endl;
+ // Write GSI input file
+ TIter it(arrOut);
+ AliParticleYield * part = 0;
+ ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
+ while ((part = (AliParticleYield*) it.Next())){
+ fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
+ }
+ fout.close();
+ // Write thermus file
+ AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff