From: mfloris Date: Mon, 14 Apr 2014 07:23:27 +0000 (+0200) Subject: Macros improvements X-Git-Url: http://git.uio.no/git/?p=u%2Fmrichter%2FAliRoot.git;a=commitdiff_plain;h=ad431d97288dcac92b6709fbb7e36dcc5755a906 Macros improvements Added further scenarios to the Macros used to extrapolate and plot ratios and yields --- diff --git a/PWGLF/ThermalFits/AverageAndExtrapolate.C b/PWGLF/ThermalFits/AverageAndExtrapolate.C index d16c421fa7b..894fbd80a76 100644 --- a/PWGLF/ThermalFits/AverageAndExtrapolate.C +++ b/PWGLF/ThermalFits/AverageAndExtrapolate.C @@ -44,6 +44,12 @@ void AverageAndExtrapolate () { // 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"); } @@ -292,6 +298,123 @@ void AverageAndExtrapolate(TString what) { } + + 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(); diff --git a/PWGLF/ThermalFits/InterpolateRatiosAndYields.C b/PWGLF/ThermalFits/InterpolateRatiosAndYields.C index 3996bdbe6e3..63d5581555d 100644 --- a/PWGLF/ThermalFits/InterpolateRatiosAndYields.C +++ b/PWGLF/ThermalFits/InterpolateRatiosAndYields.C @@ -6,17 +6,21 @@ #endif TClonesArray * arr =0; -void InterpolateRatios(Int_t pdg1, Int_t pdg2) ; +void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1="V0M0005", TString centr2="V0M0510", TString centrfinal="V0M0010") ; void Interpolate0010(Int_t pdg) ; -void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest); +void Interpolate6080(Int_t pdg) ; +void ExtrapolateWithConstantRatioToPions(Int_t pdg, TString centrOrigin, TString centrDest); +Int_t collSystem = -1; +Float_t energy = -1; void InterpolateRatiosAndYields() { #if !(!defined (__CINT__) || (defined(__MAKECINT__))) LoadLibs(); #endif + collSystem = 2; energy =2760; // *************** pi, K, pr ***************** - // arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"); + arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"); // Interpolate0010(211); // Interpolate0010(-211); // Interpolate0010(321); @@ -25,6 +29,15 @@ void InterpolateRatiosAndYields() { // Interpolate0010(-2212); // InterpolateRatios(2212,211); // InterpolateRatios(321,211); + // Interpolate6080(211); + // Interpolate6080(-211); + // Interpolate6080(2212); + // Interpolate6080(-2212); + Interpolate6080(321); + Interpolate6080(-321); + // InterpolateRatios(2212,211, "V0M6070", "V0M7080", "V0M6080"); + // InterpolateRatios(321,211, "V0M6070", "V0M7080", "V0M6080"); + // *************** Lambda and K0 ***************** // arr= AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"); // Interpolate0010(3122); @@ -34,9 +47,21 @@ void InterpolateRatiosAndYields() { // 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"); + // arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"); + // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("./PbPb_2760_AveragedNumbers.txt")); + // ExtrapolateWithConstantRatioToPions(313, "V0M0020", "V0M0010"); + + // *************** pPb, deuteron ********************* + // collSystem = 1; energy = 5020; + // 1. Average pions + //arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"); + // Interpolate0010(211); + // Interpolate0010(-211); + // 2. Extrapolate the deuteron + // arr = AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"); + // arr->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt")); + // ExtrapolateWithConstantRatioToPions(1000010020, "V0A0010", "V0A0005"); + // ExtrapolateWithConstantRatioToPions(-1000010020, "V0A0010", "V0A0005"); } @@ -54,8 +79,23 @@ void DUMP(){ void Interpolate0010(Int_t pdg) { - AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, 2, 2760., "V0M0005"); - AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, 2, 2760., "V0M0510"); + TString centrPrefix = collSystem == 2 ? "V0M" : "V0A"; + + AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0005"); + AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"0510"); + 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 Interpolate6080(Int_t pdg) { + + TString centrPrefix = collSystem == 2 ? "V0M" : "V0A"; + + AliParticleYield * p0 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"6070"); + AliParticleYield * p1 = AliParticleYield::FindParticle(arr, pdg, collSystem, energy, centrPrefix+"7080"); p0->Scale(0.5); p1->Scale(0.5); AliParticleYield * pav = AliParticleYield::Add(p0,p1); @@ -64,18 +104,18 @@ void Interpolate0010(Int_t pdg) { } -void InterpolateRatios(Int_t pdg1, Int_t pdg2) { +void InterpolateRatios(Int_t pdg1, Int_t pdg2, TString centr1, TString centr2, TString centrfinal) { 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); + ratio[0] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr1, 1); + ratio[1] = AliParticleYield::FindRatio(arr, pdg1, pdg2, 2, 2760., centr2, 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)); + pi[0] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr1, 0); + pi[0] = AliParticleYield::Add(pi[0],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr1, 0)); + pi[1] = AliParticleYield::FindParticle(arr, pdg2, 2, 2760., centr2, 0); + pi[1] = AliParticleYield::Add(pi[1],AliParticleYield::FindParticle(arr, -pdg2, 2, 2760., centr2, 0)); // Scale to get protons with errors corresponding to the ratio @@ -89,7 +129,7 @@ void InterpolateRatios(Int_t pdg1, Int_t pdg2) { pi[0]->SetSystError(0); ratio[0]->Scale(1./pi[0]->GetYield()); - ratio[0]->SetCentr("V0M0010"); + ratio[0]->SetCentr(centrfinal); ratio[0]->Print(); // average->Dump(); @@ -100,9 +140,9 @@ void InterpolateRatios(Int_t pdg1, Int_t pdg2) { 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); + AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg,collSystem, energy, centrOrigin); + AliParticleYield * pionOrigin = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrOrigin); + AliParticleYield * pionDest = AliParticleYield::FindParticle(arr, 211,collSystem, energy, centrDest); if(!part || !pionOrigin | !pionDest) { return; } diff --git a/PWGLF/ThermalFits/PlotRatiosForQM14.C b/PWGLF/ThermalFits/PlotRatiosForQM14.C index a16b40a42ab..0d506b0ca7e 100644 --- a/PWGLF/ThermalFits/PlotRatiosForQM14.C +++ b/PWGLF/ThermalFits/PlotRatiosForQM14.C @@ -9,13 +9,16 @@ 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) ; +TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ; +TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ; 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(); +#if !(!defined (__CINT__) || (defined(__MAKECINT__))) + LoadLibs(); +#endif TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt"); arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt")); arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt")); @@ -25,25 +28,66 @@ TClonesArray * PlotRatiosForQM14() { arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt")); arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt")); - TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt"); + TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt"); + + TClonesArray * arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt"); + TClonesArray * arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt"); + + TClonesArray * arrpPb = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt"); + arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt")); + arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt")); + + TClonesArray * arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt"); - // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0); - // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1); + // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges0); + // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges1); + // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0); + // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*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(""); + c1->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368); + + // c1->SetLogy(); + // CENTRAL + TH1 * h = GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%"); + h->GetYaxis()->SetRangeUser(0, 0.4); + h->Draw(); + // //GetHistoRatios(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same"); + // GetHistoRatios(arrpp7, AliParticleYield::kCSpp, 7000, "" , "pp #sqrt{s} = 7 TeV" )->Draw("same"); + // GetHistoRatios(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same"); + // // GetHistoRatios(arrpp276, AliParticleYield::kCSpp, 2760, "" , "pp #sqrt{s} = 2.76 TeV" )->Draw("same"); + // // GetHistoRatios(arrpp900, AliParticleYield::kCSpp, 900, "" , "pp #sqrt{s} = 0.9 TeV" )->Draw("same"); + // NewLegend("", "lp",0,1,1); + + // Peripheral + GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M6080", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 60-80%")->Draw("same"); + //GetHistoRatios(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same"); + // GetHistoRatios(arrpp7, AliParticleYield::kCSpp, 7000, "" , "pp #sqrt{s} = 7 TeV" )->Draw("same"); + GetHistoRatios(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same"); + // GetHistoRatios(arrpp276, AliParticleYield::kCSpp, 2760, "" , "pp #sqrt{s} = 2.76 TeV" )->Draw("same"); + // GetHistoRatios(arrpp900, AliParticleYield::kCSpp, 900, "" , "pp #sqrt{s} = 0.9 TeV" )->Draw("same"); + NewLegend("", "lp",0,1,1); + + //return; + TCanvas * c2 = new TCanvas("Yields", "Yields", 1400, 600); + c2->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368); + + GetHistoYields(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%")->Draw(); + GetHistoYields(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same"); + NewLegend("", "lp",0,1,1); return arr; } -TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) { +TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) { + // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY 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); + Int_t num [nratio] = {kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar}; + // Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGPi , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK}; + 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 }; + Double_t scale[nratio] = {1 ,3 ,1 ,30 ,250 ,50 ,100 ,4e5 ,2 ,2 }; + TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, 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++){ @@ -55,23 +99,49 @@ TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString 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 + // Try with the !sum, if part 1 is not found + if(!part1) { part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]); + // If the sum was requested, try to recover it! + if(isSum[iratio-1]) { + AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0); + if(part1 && part1_bar) { + std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl; + + part1 = AliParticleYield::Add(part1, part1_bar); + + } + } else if(part1) { + // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable + part1->Scale(0.5); + } + } 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]); + if(isSum[iratio-1]) { + AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0); + if(part2 && part2_bar) part2 = AliParticleYield::Add(part2, part2_bar); + } else if(part2){ + // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable + part2->Scale(0.5); + } + } - ratio = AliParticleYield::Divide(part1, part2); + ratio = AliParticleYield::Divide(part1, part2, 0, "YQ"); if(ratio) { + std::cout << "" << std::endl; std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<Print("short"); + ratio->Print(""); + std::cout << "" << std::endl; } } if(ratio){ + ratio->Scale(scale[iratio-1]); 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())); + h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f", scale[iratio-1]), ratio->GetLatexName())); } else { h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1])); @@ -136,3 +206,27 @@ void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t // Write thermus file AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges)); } + + +TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) { + + const Int_t npart = 11; + Int_t pdg [npart] = {kPDGPi, kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar}; + // Int_t isSum[npart] = {1 ,1 ,1 ,0 ,1 ,1 ,0 ,0 ,1 ,0 ,1 }; + Int_t isSum[npart] = {0,0,0,0,0,0,0,0,0,0,1}; + // Double_t scale[npart] = {1 ,1 ,3 ,1 ,30 ,250 ,50 ,10 ,4e5 ,2 ,2 }; + Double_t scale[npart] = {1,5,30,30,200,1000,4000,2e6,2e6,20,20,}; + TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1, npart+1); + + // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored + for(Int_t ipart = 1; ipart <= npart; ipart++){ + AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg[ipart-1], system, energy, centrality,isSum[ipart-1]); + if(!part) continue; + part->Scale(scale[ipart-1]); + h->SetBinContent(ipart, part->GetYield()); + h->SetBinError (ipart, part->GetTotalError(0/* 0 = no normalization error */)); + h->GetXaxis()->SetBinLabel(ipart, Form("#splitline{%s}{%s}",Form("#times%2.0g", scale[ipart-1]), part->GetLatexName())); + + } + return h; +}