Macros improvements
authormfloris <michele.floris@cern.ch>
Mon, 14 Apr 2014 07:23:27 +0000 (09:23 +0200)
committermfloris <michele.floris@cern.ch>
Mon, 14 Apr 2014 07:24:45 +0000 (09:24 +0200)
Added further scenarios to the Macros used to extrapolate and plot ratios and yields

PWGLF/ThermalFits/AverageAndExtrapolate.C
PWGLF/ThermalFits/InterpolateRatiosAndYields.C
PWGLF/ThermalFits/PlotRatiosForQM14.C

index d16c421..894fbd8 100644 (file)
@@ -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();
index 3996bdb..63d5581 100644 (file)
@@ -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;
   }
index a16b40a..0d506b0 100644 (file)
@@ -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] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
-        ratio->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;
+}