PWGLF/ThermalFits/PlotRatiosForQM14.C

   1 #if !defined (__CINT__) || (defined(__MAKECINT__))
   2 #include <iostream>
   3 #include "TClonesArray.h"
   4 #include "AliParticleYield.h"
   5 #include "TH1F.h"
   6 #include "TCanvas.h"
   7 #include "TStyle.h"
   8 #include <fstream>
   9 #include "TLatex.h"
  10 #include "TLegend.h"
  11 #include "TList.h"
  12 #include "TF1.h"
  13 #include "AliPWGHistoTools.h"
  14 #include "TGraphErrors.h"
  15 #include "TMath.h"
  16 #include "TDatabasePDG.h"
  17 #include "TH2F.h"
  18 #include "TSystem.h"
  19 #include "TPaveText.h"
  20 #include <map>
  21 #include <fstream>
  22 #include <istream>
  23 #include "TMarker.h"
  24 #include "TObjString.h"
  25 #include "TLegendEntry.h"
  26
  27 #endif
  28
  29 // Plots ratios for QM and saves input files for thermal models
  30
  31 #if 1  // DUMMY IFDEF USED TO HIDE PREAMBLE in EMACS
  32
  33 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};
  34
  35 typedef enum {kStatError, kSystError, kTotalError} myerror_t;
  36
  37 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle, Int_t icolor, Int_t imarker = kOpenSquare, Int_t errorsType = kTotalError, Float_t shift = 0) ;
  38 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle, Int_t icolor, Int_t imarker = kOpenSquare, Int_t errorsType = kTotalError, Float_t shift = 0) ;
  39 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
  40 void SetStyle(Bool_t graypalette=0) ;
  41 TLegend * NewLegendQM(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Bool_t isYield = 0) ;
  42 void DrawRatio(TString what, Bool_t isYield = kFALSE, Double_t shiftloc=0.);
  43
  44 void DrawFrame(Bool_t yields = 0) ;
  45 void DrawExtrapolatedSymbolsAndLegendPbPb0010() ;
  46 void DrawExtrapolatedSymbolsAndLegendpPb0005() ;
  47 void DrawMarkerKStarNoFit(Bool_t plotLegend = 0) ;
  48 void DrawMarkerNucleiNoFit() ;
  49 void DrawExtrapNotInFitpPb0005(Bool_t drawExtrap = 1) ;
  50
  51 void DrawExtrapolatedSymbolsYieldsPbPb0010(Double_t x1=0.144578, Double_t y1=0.408249, Double_t x2=0.351406, Double_t y2=0.542403, Bool_t plotExtraploatedLegend=1);
  52 Float_t shiftRatioDataModel = 0;
  53 Double_t GetGraphRatioAndStdDev(TGraphErrors * gModel, TGraphErrors * &gRatio, TGraphErrors *&gStdDev) ;
  54 TString particlesToExcludeFromChi2 = "";// The above method recomputes the chi2. This string is used to store values to be excluded from this calculation. Each PDG cocde has to be enclosed in [..]. This is obviously not efficient as it involves many string conversions. But efficiency is not an issue here.
  55
  56
  57 void LoadArrays() ;
  58 //void AddLabel(Float_t x, Float_t y, TString text);
  59 void myLatexDraw(TLatex *currentLatex, Float_t currentSize=0.5, Int_t currentColor=1);
  60 void myPaveSetup(float rRatio=0, float rRange3=0, float rRange5=0,
  61                  int rFillColor=0);
  62 void myPadSetUp(TPad *currentPad);
  63 TGraphErrors*  PlotThermusYields(const char * filename, Int_t color, Int_t lineStyle,
  64                                  const char * tag);
  65 TGraphErrors * PlotGSIYields(const char * fileName, Int_t color=kBlack, Int_t lineStyle = kSolid,
  66                              const char * tag ="", Bool_t isPbPb = 1);
  67 TGraphErrors*  PlotFlorenceYields(const char * filename, Int_t color, Int_t lineStyle,
  68                                   const char * tag) ;
  69
  70 void AddLineToThermalLegend(TObject * obj, TString line, const char * optFirst = "L");
  71
  72 void SaveCanvas(const char * name) ;
  73
  74
  75 // Ratios to be draw. Remember to change the labels in DrawFrame if you change this
  76 const Int_t nratio = 10;
  77 Int_t num  [nratio]            = {kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
  78 Int_t denum[nratio]            = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGProton   , kPDGDeuteron , kPDGPi          , kPDGK   , kPDGK};
  79 Int_t isSum[nratio]            = {1      , 1          ,  1         ,   1     , 1         , 1            , 0            , 1               , 1       , 1      };
  80
  81 // const Int_t nratio = 10;
  82 // Int_t num  [nratio]            = {kPDGK  , kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron , kPDGHE3      , kPDGHyperTriton , kPDGPhi , kPDGKStar};
  83 // Int_t denum[nratio]            = {kPDGPi , kPDGPi     , kPDGKS0    ,  kPDGPi , kPDGPi    , kPDGProton   , kPDGDeuteron , kPDGPi          , kPDGK   , kPDGK};
  84 // Int_t isSum[nratio]            = {1      , 1          ,  1         ,   1     , 1         , 1            , 0            , 1               , 1       , 1      };
  85
  86 const char * ratiosLabels[]          = {"#frac{K^{+}+K^{-}}{#pi^{+}+#pi^{-}}",
  87                                         "#frac{p+#bar{p}}{#pi^{+}+#pi^{-}}",
  88                                         "#frac{2#Lambda}{K_{S}^{0}}",
  89                                         "#frac{#Xi^{-}+#Xi^{+}}{#pi^{+}+#pi^{-}}",
  90                                         "#frac{#Omega^{-}+#Omega^{+}}{#pi^{+}+#pi^{-}}",
  91                                         "#frac{d}{p+#bar{p}}",
  92                                         "#frac{{}^{3}He }{d}",
  93                                         "#frac{{}^{3}_{#Lambda}H+{}^{3}_{#Lambda}#bar{H} }{#pi^{+}+#pi^{-}}",
  94                                         "#frac{#phi}{K^{+}+K^{-}}",
  95                                         "#frac{K*+#bar{K}*}{K^{+}+K^{-}}",};
  96 //static const Double_t scale[]  = {1      , 3          ,  0.5       ,  30     ,  250      , 50           , 100          , 4e5             , 2       , 1      };
  97 static const Double_t scale[]  = {1      , 3          ,  0.5       ,  80     ,  1000      , 50           , 100          , 4e5             , 2       , 1      };
  98 //static const Double_t scale[]  = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,};
  99 const Int_t npart = 12;
 100 Int_t particleYields  [npart] = {kPDGPi ,kPDGK   ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi  , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3    };
 101 Int_t isSumYields[npart]      = {1      ,1       ,0      , 1        , 0      , 1           ,1            ,1        ,1          ,0           , 1              , 0          };
 102 //Int_t isSumInputFiles[npart]  = {1      ,1       ,0      , 1        , 0      , 1           ,1            ,1        ,1          ,0           , 1              , 0          };
 103 const char * yieldsLabels[]          = {"#frac{#pi^{+}+#pi^{-}}{2}",
 104                                         "#frac{K^{+}+K^{-}}{2}",
 105                                         "K_{S}^{0}",
 106                                         "#frac{K*+K*}{2}",
 107                                         "#phi",
 108                                         "#frac{p+#bar{p}}{2}",
 109                                         "#Lambda",
 110                                         "#frac{#Xi^{-}+#Xi^{+}}{2}",
 111                                         "#frac{#Omega^{-}+#Omega^{+}}{2}",
 112                                         "d",
 113                                         "#frac{{}^{3}_{#Lambda}H+{}^{3}_{#Lambda}#bar{H}}{2}",
 114                                         "{}^{3}He",
 115 };
 116
 117 //
 118 Int_t markerNoFit  = 28;
 119 Int_t markerExtrap = 27;
 120 Double_t maxy = 0.5;
 121
 122
 123 // Data arrays;
 124 TClonesArray *arrPbPb=0, *arrpp7=0, *arrpPb=0, * arrpp276=0, * arrpp900=0, * arrThermus=0;
 125 TClonesArray *arrSTARPbPb=0, *arrPHENIXPbPb=0, *arrBRAHMSPbPb=0;
 126 TClonesArray *arrSTARpp  =0, *arrPHENIXpp=0;
 127
 128 //const Double_t *scaleRatios = 0;
 129 Double_t *correlatedUnc = 0;
 130 //Double_t correlatedUncLocalPbPb[14] = {0.0424 , 0.0424     ,  0.041     ,  0      , 0         , 0.0424       , 0.0424       , 0               , 0.05    , 0.05   };
 131 Double_t correlatedUncLocalPbPbOnlyKStarPhi[14] = {0. , 0.     ,  0.     ,  0      , 0         , 0.       , 0.       , 0               , 0.05    , 0.05   };
 132 Double_t correlatedUncLocalPP  [14] = {0.0424 , 0.0424     ,  0.041     ,  0      , 0         , 0.0424       , 0.0424       , 0               , 0.0424    , 0.0424   };
 133 Double_t correlatedUncZero[14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0};
 134
 135 TCanvas *myCan = 0;
 136 TPad    *myPadStdDev =0;
 137 TPad    *myPadRatio  =0;
 138 TPad    *myPadHisto  =0;
 139 TPad    *myPadLabel  =0;
 140 TLegend * legThermal = 0;
 141 #endif
 142
 143 Bool_t saveCanvas = 1; // if true, the canvas is saved and copied to the analysis note folder
 144
 145 TClonesArray * PlotRatiosForQM14() {
 146 #if !(!defined (__CINT__) || (defined(__MAKECINT__)))
 147   LoadLibs();
 148 #endif
 149
 150   //
 151   LoadArrays();
 152
 153   // Uncomment stuff in this section to save the inputs for thermal models
 154   //#define SAVE_INPUT_THERMAL_MODEL
 155 #ifdef SAVE_INPUT_THERMAL_MODEL
 156     PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
 157   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/1);
 158   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", /*separateCharges*/1);
 159   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A2040", /*separateCharges*/1);
 160   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A6080", /*separateCharges*/1);
 161   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M6080", /* separateCharges*/1);
 162   //  PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M2040", /*separateCharges*/1);
 163
 164   PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
 165   // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0);
 166   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", /*separateCharges*/0);
 167   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A2040", /*separateCharges*/0);
 168   // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A6080", /*separateCharges*/0);
 169   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M6080", /*separateCharges*/0);
 170   // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M2040", /*separateCharges*/0);
 171
 172   return 0;
 173 #endif
 174
 175
 176   SetStyle();
 177
 178   DrawRatio("allpp");
 179   //  DrawRatio("allppWithRHIC");
 180   //  DrawRatio("PbPbWithPP7TeV");
 181   //DrawRatio("allsyst");
 182   //DrawRatio("PbPb6080andpPb0005");
 183   //  DrawRatio("pp_vsRHIC");
 184   //  DrawRatio("PbPb_vsRHIC");
 185   //DrawRatio("aliceall");
 186
 187
 188   // Yields and FITS
 189   //  maxy=2000;
 190
 191   // DrawRatio("fit_ReferenceFit_PbPb0010", 1);
 192   //  DrawRatio("fit_ReferenceFit_GSIONLY_PbPb0010", 1);
 193   //  DrawRatio("fit_ReferenceFit_GSITHERMUS_PbPb0010",1);
 194   // DrawRatio("fitSHARE_NoPionsNoProtons_PbPb0010",1);
 195   //  DrawRatio("fitGSI_NoPionsNoProtons_PbPb0010", 1);
 196   // DrawRatio("fitShare_All_PbPb0010", 1);
 197
 198   //  DrawRatio("fitShareWithWithoutNuclei_PbPb0010", 1);
 199   // maxy=200;
 200   // DrawRatio("fitGSI_PbPb6080", 1);
 201   //   maxy=150;
 202   //   DrawRatio("fitGSI_PbPb2040", 1);
 203   //  maxy = 60;
 204   // DrawRatio("fitThermus_GammaSFree_pPb0005");
 205   //  DrawRatio("fitShare_pPb0005");
 206   //  DrawRatio("fitShare_pPb0005_NoOmega", 1);
 207   //  maxy=20;
 208   //  DrawRatio("fitThermus_GammaSFree_pPb2040");
 209   //  maxy=9;
 210   //  DrawRatio("fitThermus_GammaSFree_pPb6080");
 211   //  maxy=9;
 212   //    DrawRatio("fitGSI_pp");
 213   //  DrawRatio("fitFlorence_pp");
 214
 215   //  NewLegendQM();
 216   return arrPbPb;
 217 }
 218
 219 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle, Int_t icolor, Int_t imarker, Int_t errorType, Float_t shift) {
 220   // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY
 221
 222
 223   //  scaleRatios = scale;
 224   TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s_%d",system,energy,centrality.Data(),errorType), histotitle, nratio, 1+shift, nratio+1+shift);
 225
 226   TClonesArray * arrRatios = new TClonesArray ("AliParticleYield");// We save to disk the precomputed ratios
 227   Int_t iratioArr = 0;// Index used only to add particles to the array above
 228
 229   //  Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
 230   for(Int_t iratio = 1; iratio <= nratio; iratio++){
 231     std::cout << "******** " << num[iratio-1] << "/" <<  denum[iratio-1]<< " ("<<isSum[iratio-1]<<")*******" << std::endl ;
 232
 233     AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
 234     if(ratio)
 235       {
 236         ratio = new AliParticleYield(*ratio); // We need to clone it to avoid a mess if we need to use this particle again later (e.g. double scalings)
 237         std::cout << "  " ;
 238         ratio->Print("short");
 239       }
 240
 241     if(!ratio) {
 242       // If the ratio is not found, try to build it!
 243       std::cout << "  Looking for " <<  num[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
 244       AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,  isSum[iratio-1]);
 245       if(part1) {
 246         part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
 247         if(isSum[iratio-1] && part1->IsTypeAverage()) {
 248           std::cout << "Sum requested, found average, scaling x2" << std::endl;
 249           part1->Scale(2.);
 250         }
 251       }
 252       // Try with the !sum, if part 1 is not found
 253       if(!part1) {
 254         std::cout << "  Looking for " <<  num[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
 255         part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 256         if(part1)
 257           {
 258             part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
 259             // If the sum was requested, try to recover it!
 260             if(isSum[iratio-1]) {
 261               std::cout << "  Looking for " <<  -num[iratio-1] <<std::endl;
 262               // If it's lambda and ALICE, use 2L instead of L + Lbar // FIXME: do the same for deuterons?
 263               if((num[iratio-1] == kPDGLambda || num[iratio-1] == kPDGDeuteron) && energy > 300) {
 264                 std::cout << "   It's Lambda or deuteron ALICE: Scaling x2 " << std::endl;
 265                 part1->Print();
 266                 part1->Scale(2.);
 267               } else {
 268                 AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0);
 269                 if(part1 && part1_bar) {
 270                   std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl;
 271                   part1 = AliParticleYield::Add(part1, part1_bar);
 272
 273                 } else if (TDatabasePDG::Instance()->GetParticle(-num[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
 274                   std::cout << "WARNING: Sum requested but not found, scaling x2 " << part1->GetName() << std::endl;
 275                   part1->Scale(2);
 276                 }
 277               }
 278             } else if(part1) { // if we are here, it means the sum was *not* requested (isSum=0), but we found something with (!isSum) = 1
 279               // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 280               std::cout << "WARNING: Using sum/2 for " << part1->GetName() << std::endl;
 281
 282               part1->Scale(0.5);
 283             }
 284         }
 285
 286       }
 287
 288
 289       std::cout << "  Looking for " <<  denum[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
 290       AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
 291       if(part2) {
 292         part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
 293         if(isSum[iratio-1] && part2->IsTypeAverage()) {
 294           std::cout << "Sum requested, found average, scaling x2" << std::endl;
 295           part2->Scale(2.);
 296         }
 297       }
 298       if(!part2) {// Try with the !sum, if part 2 is not found
 299         std::cout << "  Looking for " <<  denum[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
 300         part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
 301         if(part2)
 302           {
 303             part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
 304             if(isSum[iratio-1]) {
 305               std::cout << "  Looking for " <<  -denum[iratio-1] << std::endl;
 306               AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0);
 307               if(part2 && part2_bar){
 308                 std::cout << "Adding " << part2_bar->GetPartName() << " " << part2->GetPartName() << std::endl;
 309                 part2 = AliParticleYield::Add(part2, part2_bar);
 310               } else if (TDatabasePDG::Instance()->GetParticle(-denum[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
 311                 std::cout << "WARNING: Sum requested but not found, scaling x2 " << part2->GetName() << std::endl;
 312                 part2->Scale(2);
 313               }
 314             } else if(part2){
 315               // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
 316               std::cout << "WARNING: Using sum/2 for " << part2->GetName() << std::endl;
 317               part2->Scale(0.5);
 318             }
 319           }
 320
 321       }
 322       ratio = AliParticleYield::Divide(part1, part2, correlatedUnc[iratio-1], "YQ"); // Assume by that the systematics of part1 and part2 are uncorrelated.
 323
 324       if(ratio) {
 325         std::cout << "" << std::endl;
 326         std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
 327         ratio->Print("");
 328         std::cout << "" << std::endl;
 329       }
 330     }
 331     if(ratio){
 332       ratio->Scale(scale[iratio-1]);
 333       h->SetBinContent(iratio, ratio->GetYield());
 334       if(errorType == kTotalError) {
 335         h->SetBinError  (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
 336       } else if (errorType == kStatError) {
 337         h->SetBinError  (iratio, ratio->GetStatError());
 338       } else if (errorType == kSystError) {
 339         h->SetBinError  (iratio, ratio->GetSystError());
 340       } else {
 341         std::cout <<  "ERROR: Unknown Error Type " << errorType << std::endl;
 342       }
 343
 344       //      h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f",  scale[iratio-1]), ratio->GetLatexName()));
 345       h->GetXaxis()->SetBinLabel(iratio, ratio->GetLatexName());
 346       new ((*arrRatios)[iratioArr++]) AliParticleYield(*ratio);
 347     }
 348     else {
 349       h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
 350
 351     }
 352     std::cout << "*** END OF " << num[iratio-1] << "/" <<  denum[iratio-1]<< " *******" << std::endl ;
 353
 354   }
 355
 356   h->GetYaxis()->SetRangeUser(0, maxy);
 357   h->SetLineColor(icolor);
 358   h->SetMarkerColor(icolor);
 359   h->SetMarkerStyle(imarker);
 360   // the "if" avoids saving twice the same ratios
 361   if(errorType == kSystError) AliParticleYield::SaveAsASCIIFile(arrRatios, TString("ratios_")+h->GetName());
 362   return h;
 363
 364
 365
 366 }
 367
 368 void   PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges)  {
 369   // If "Separate charges" is true, tries to dump both charges are dumped
 370   TClonesArray * arrOut = new TClonesArray("AliParticleYield");
 371   TClonesArray * arrOutGSI = new TClonesArray("AliParticleYield"); // We add dummy lines to the GSI output file if needed!
 372 ;
 373
 374   Int_t ipartOut = 0; // Index for the array
 375   Int_t ipartOutGSI = 0; // Index for the array
 376
 377   for(Int_t ipart = 0; ipart < npart; ipart++){
 378     if(!separateCharges) {
 379       AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality,  isSumYields[ipart]);
 380       if(!part && isSumYields[ipart]) {
 381         //Could not find the particle, but the sum was requested: build the sum!
 382         part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality,  0);
 383         AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particleYields[ipart], system, energy, centrality,  0);
 384         if(part2 && part) part = AliParticleYield::Add(part, part2);
 385         else if(part) part->Scale(2.); // If we only found a particle, we can scale it by a factor 2.
 386         else part = 0;
 387       }
 388       // We want to save the average of particle and antiparticle in this case
 389       if(part) {
 390         if(isSumYields[ipart] && !part->IsTypeAverage()) part->Scale(0.5); // If it's not already an average, but just a sum, divide by 2
 391         new((*arrOut   )[ipartOut++]) AliParticleYield(*part);
 392         new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
 393       } else { // Add dummy particle to the GSI list
 394         new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
 395       }
 396     }
 397     else {
 398       // ignore isSumYields and try to find both particleYields
 399       Bool_t notFound = 0;
 400       AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality,  0);
 401       if(part) {
 402         new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 403         new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
 404       }
 405       else {
 406         // std::cout << "ADDING DUMMY part " << particleYields[ipart] << std::endl;
 407         // new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
 408         notFound=1;
 409       }
 410       // Try to find antiparticle (-pdg code)
 411       part = AliParticleYield::FindParticle(arr, -particleYields[ipart], system, energy, centrality,  0);
 412       if(part) {
 413         new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 414         new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
 415       }
 416       else if (notFound) {
 417         // If neither charge was found, check if we at least have the sum
 418         part = AliParticleYield::FindParticle(arr, abs(particleYields[ipart]), system, energy, centrality,  1);
 419         if (part) {
 420           if(!part->IsTypeAverage()) part->Scale(0.5); // If it's a sum (not an average) divide by 2
 421           new((*arrOut)[ipartOut++]) AliParticleYield(*part);
 422           new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
 423           if(TDatabasePDG::Instance()->GetParticle(-particleYields[ipart]) &&
 424              (particleYields[ipart] != kPDGLambda) &&
 425              (particleYields[ipart] != kPDGKStar)
 426              ){// if only the sum was found, add a dummy entry to the
 427                // GSI file, so that anton always has the same # of
 428                // lines. However, we don't do this for the Lambda and
 429                // KStar (always one of the charges or the average)
 430             new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(-particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
 431           }
 432         }
 433         else {
 434           std::cout << "ADDING DUMMY sum " << particleYields[ipart] << std::endl;
 435           new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
 436           if (particleYields[ipart]==kPDGHyperTriton) {
 437             // If is the 3LH, add another one for the antiparticle
 438             AliParticleYield(-particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
 439           }
 440         }
 441       }
 442
 443     }
 444   }
 445   std::cout << "Particles for thermal model fits:" << std::endl;
 446   arrOut->Print("short");
 447   //  arrOut->Print("");
 448   std::cout << "" << std::endl;
 449   // Write GSI input file
 450   TIter it(arrOutGSI);
 451   AliParticleYield * part = 0;
 452   ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 453   while ((part = (AliParticleYield*) it.Next())){
 454     fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
 455   }
 456   fout.close();
 457   // Write thermus file
 458   AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
 459 }
 460
 461
 462 TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle,
 463                       Int_t icolor, Int_t imarker, Int_t errorsType, Float_t shift) {
 464
 465   TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1+shift, npart+1+shift);
 466
 467   for(Int_t ipart = 1; ipart <= npart; ipart++){
 468     std::cout << "----- Searching " << particleYields[ipart-1] << " -------"  << std::endl;
 469
 470     AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart-1], system, energy, centrality,isSumYields[ipart-1]);
 471     if(part) {
 472       std::cout << "found" << std::endl;
 473       part->Print();
 474     }
 475     if(!part && isSumYields[ipart-1]) {
 476       //Could not find the particle, but the sum was requested: build the sum!
 477       part = AliParticleYield::FindParticle(arr, particleYields[ipart-1], system, energy, centrality,  0);
 478       AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particleYields[ipart-1], system, energy, centrality,  0);
 479       if(part2 && part) {
 480         std::cout << " Building sum" << std::endl;
 481         part->Print();
 482         part2->Print();
 483         part = AliParticleYield::Add(part, part2);
 484       }
 485       else if(part) {
 486         std::cout << "Scaling part" << std::endl;
 487         part->Print();
 488         part = new AliParticleYield(*part); // Always clone before scaling
 489         part->Scale(2.); // If we only found a particle, we can scale it by a factor 2.
 490       }
 491       else part = 0;
 492     }
 493     if(!part){
 494       std::cout << "Cannot find " << particleYields[ipart-1] << std::endl;
 495       continue;
 496     }
 497     if(isSumYields[ipart-1] && !part->IsTypeAverage()) {
 498       std::cout << " scaling /2" << std::endl;
 499       part = new AliParticleYield(*part); // Always clone before scaling
 500       part->Scale(0.5); // take average
 501     }
 502     std::cout << " Plotting " << particleYields[ipart-1] << std::endl;
 503     part->Print();
 504     //    part->Scale(scale[ipart-1]);
 505     h->SetBinContent(ipart, part->GetYield());
 506     if(errorsType == kTotalError) {
 507       h->SetBinError  (ipart, part->GetTotalError(0/* 0 = no normalization error */));
 508     } else if (errorsType == kSystError) {
 509       h->SetBinError  (ipart, part->GetSystError());
 510     } else if (errorsType == kStatError) {
 511       h->SetBinError  (ipart, part->GetStatError());
 512     }
 513     h->GetXaxis()->SetBinLabel(ipart, part->GetLatexName());
 514
 515   }
 516   h->SetMarkerStyle(imarker);
 517   h->SetMarkerColor(icolor);
 518   h->SetLineColor(icolor);
 519   h->SetMarkerSize(1.4);
 520
 521   return h;
 522 }
 523
 524 void LoadArrays() {
 525   arrPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
 526   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
 527   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
 528   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
 529   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
 530   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
 531   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
 532   arrPbPb->AbsorbObjects(  AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
 533
 534   arrpp7   = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
 535
 536   arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt");
 537   arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt");
 538
 539   arrpPb   = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt");
 540   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
 541   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"));
 542   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_AveragedNumbers.txt"));
 543   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_phi.txt"));
 544   arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_Kstar.txt"));
 545   arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt");
 546
 547   // RHIC data
 548   arrSTARPbPb   = AliParticleYield::ReadFromASCIIFile("PbPb_200_STAR-AntonQM12.txt");
 549   arrPHENIXPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_PHENIX-AntonQM12.txt");
 550   arrBRAHMSPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_BRAHMS-AntonQM12.txt");
 551   arrSTARpp     = AliParticleYield::ReadFromASCIIFile("pp_200_STAR.txt");
 552   arrPHENIXpp   = AliParticleYield::ReadFromASCIIFile("pp_200_PHENIX.txt");
 553
 554 }
 555
 556 void SetStyle(Bool_t graypalette) {
 557   std::cout << "Setting style!" << std::endl;
 558
 559   gStyle->Reset("Plain");
 560   gStyle->SetOptTitle(0);
 561   gStyle->SetOptStat(0);
 562   if(graypalette) gStyle->SetPalette(8,0);
 563   else gStyle->SetPalette(1);
 564   gStyle->SetDrawBorder(0);
 565   gStyle->SetCanvasColor(10);
 566   gStyle->SetCanvasBorderMode(0);
 567   gStyle->SetPadBorderMode(0);
 568   gStyle->SetFrameLineWidth(1);
 569   gStyle->SetFrameFillColor(kWhite);
 570   gStyle->SetPadColor(10);
 571   gStyle->SetPadTickX(1);
 572   gStyle->SetPadTickY(1);
 573   gStyle->SetPadBottomMargin(0.15);
 574   gStyle->SetPadLeftMargin(0.15);
 575   gStyle->SetHistLineWidth(1);
 576   gStyle->SetHistLineColor(kRed);
 577   gStyle->SetFuncWidth(2);
 578   gStyle->SetFuncColor(kGreen);
 579   gStyle->SetLineWidth(1);
 580   gStyle->SetLabelSize(0.045,"yz");
 581   gStyle->SetLabelSize(0.06,"x");
 582   gStyle->SetLabelOffset(0.01,"y");
 583   gStyle->SetLabelOffset(0.01,"x");
 584   gStyle->SetLabelColor(kBlack,"xyz");
 585   gStyle->SetTitleSize(0.05,"xyz");
 586   gStyle->SetTitleOffset(1.25,"y");
 587   gStyle->SetTitleOffset(1.2,"x");
 588   gStyle->SetTitleFillColor(kWhite);
 589   gStyle->SetTextSizePixels(26);
 590   gStyle->SetTextFont(42);
 591   gStyle->SetTickLength(0.012,"X");  gStyle->SetTickLength(0.012,"Y");
 592
 593   gStyle->SetLegendBorderSize(0);
 594   gStyle->SetLegendFillColor(kWhite);
 595   //  gStyle->SetFillColor(kWhite);
 596   gStyle->SetLegendFont(42);
 597
 598   gStyle->SetErrorX(0);
 599   gStyle->SetEndErrorSize(5);
 600 }
 601
 602 TLegend * NewLegendQM(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Bool_t isYield) {
 603
 604   const char * style = "p";
 605   const char ** labels=0;
 606   //  Bool_t beautify=kFALSE;
 607   Bool_t useTitle=kTRUE;
 608
 609   TLegend * l = new TLegend(x1, y1, x2, y2);
 610   l->SetFillColor(kWhite);
 611   l->SetTextFont(43);
 612   l->SetTextSize(25);
 613   l->SetBorderSize(1);
 614   l->SetLineWidth(1);
 615   l->SetMargin(0.1);
 616   // const Int_t markers[] = {20,24,21,25,23,28,33,20,24,21,25,23,28,33};
 617   // const Int_t colors[]  = {1,2,3,4,6,7,8,9,10,11,1,2,3,4,6,7,8,9,10};
 618
 619   TList * list = gPad->GetListOfPrimitives();
 620   TIterator * iter = list->MakeIterator();
 621   TObject * obj = 0;
 622   Int_t ilabel = -1;
 623   while ((obj = (TObject*) iter->Next())){
 624     if (obj->InheritsFrom("TH1") || obj->InheritsFrom("TGraph") || obj->InheritsFrom("TF1")) {
 625       if( (TString(obj->GetName()) == "hframe" ) ) continue;
 626       ilabel++;
 627       if (labels != 0)
 628         l->AddEntry(obj, labels[ilabel], style);
 629       else{
 630         if (useTitle)  {
 631           if(TString(obj->GetTitle()).Contains("NoLegend")) continue;
 632           TString title = obj->GetTitle();
 633           if(title.Contains("p-Pb")) {
 634             TObjArray * tokens = title.Tokenize(",");
 635             if(tokens) {
 636               TString system = ((TObjString*) tokens->At(0))->String().Strip(TString::kBoth, ' ');
 637               TString centr  = ((TObjString*) tokens->At(1))->String().Strip(TString::kBoth, ' ');
 638               l->AddEntry(obj, system, style);
 639               l->AddEntry(obj, centr, "0");
 640               delete tokens;
 641             }
 642
 643           } else {
 644             l->AddEntry(obj, obj->GetTitle(), style);
 645           }
 646
 647         }
 648         else
 649           l->AddEntry(obj, obj->GetName(), style);
 650       }
 651       // if(beautify) {
 652
 653       //   if(!obj->InheritsFrom("TF1")){
 654       //     ((TH1*)obj)->SetLineColor(colors[ilabel]);
 655       //     ((TH1*)obj)->SetMarkerColor(colors[ilabel]);
 656       //     ((TH1*)obj)->SetMarkerStyle(markers[ilabel]);
 657       //   } else {
 658       //     ((TF1*)obj)->SetLineColor(colors[ilabel]);
 659       //   }
 660
 661     }
 662   }
 663   if(isYield) {
 664     // Add some details on excluded stuff
 665     l->SetTextSize(22);
 666     l->SetBorderSize(0);
 667     l->GetEntry()->SetOption("0");
 668     l->SetMargin(0.01);
 669   }
 670     l->Draw();
 671   return l;
 672 }
 673
 674
 675
 676 void DrawRatio(TString what, Bool_t isYield, Double_t shift) {
 677   // This is used to simplify the code above
 678   // In order to draw syst error bars, we need to convert to graphs the syst errors histos
 679   // if isYield == true plots yields rather than ratios
 680
 681   // Sample colors
 682   //  const Int_t colors[]     = {kBlack, kRed+1 , kBlue+1, kGreen+3, kMagenta+1, kOrange-1,kCyan+2,kYellow+2  , kWhite};
 683
 684   TClonesArray * array = 0;
 685   Int_t system,  color, marker;
 686   Float_t energy = 0, shiftloc = shift;
 687   TString centrality, label;
 688   // FIXME: move this in the different sections below
 689   correlatedUnc = correlatedUncZero;
 690   std::cout << "Plotting " << what.Data() << std::endl;
 691
 692
 693   if (what == "frame" ) {
 694     correlatedUnc = correlatedUncZero;
 695     DrawFrame(isYield);
 696     // TH1 * h = GetHistoRatios(arrPbPb,       AliParticleYield::kCSPbPb, 2760, "V0M0010", "NoLegend", kWhite);
 697     // h->Draw();
 698     // h->GetYaxis()->SetDecimals(1);
 699     // h->GetYaxis()->SetNdivisions(505);
 700     //    h->GetXaxis()->CenterLabels(1);
 701     //    Int_t nratio = h->GetNbinsX();
 702     if(!isYield) {
 703       for(Int_t iratio = 0; iratio < nratio; iratio++){
 704         Double_t exp = TMath::Floor(TMath::Log10(TMath::Abs(scale[iratio])));
 705         Double_t man = scale[iratio] / TMath::Power(10, exp);
 706         if(exp > 2) {
 707           //        TLatex * scaleLabel = new TLatex(iratio+1+0.2,maxy*1.01, Form("#times %0.0f 10^{%0.0f}", man, exp));
 708           TLatex * scaleLabel = new TLatex(iratio+1+0.2,0.005, Form("#times %0.0f 10^{%0.0f}", man, exp));
 709           scaleLabel->SetTextFont(43);
 710           scaleLabel->SetTextSize(20);
 711           scaleLabel->Draw();
 712         } else {
 713           Double_t shiftloc2 = scale[iratio] < 50 ? 0.3 : 0.2;
 714           TLatex * scaleLabel = new TLatex(iratio+1+shiftloc2,  0.005, Form("#times %g", scale[iratio]));
 715           scaleLabel->SetTextFont(43);
 716           scaleLabel->SetTextSize(20);
 717           scaleLabel->Draw();
 718         }
 719
 720       }
 721     }
 722     if(isYield) {
 723       TPaveText *pt = new TPaveText(0.691767, 0.86069, 0.893574, 0.944865,"brNDC");
 724       pt->SetBorderSize(0);
 725       pt->SetFillColor(0);
 726       pt->SetTextAlign(12);
 727       pt->SetTextFont(43);
 728       pt->SetTextSize(23);
 729       pt->AddText("ALICE Preliminary");
 730       pt->Draw();
 731     }
 732     else {
 733       TPaveText *pt = new TPaveText(    0.176, 0.842881, 0.378514, 0.929595,"brNDC");
 734       pt->SetBorderSize(0);
 735       pt->SetFillColor(0);
 736       pt->SetLineWidth(1);
 737       pt->SetTextFont(43);
 738       pt->SetTextSize(23);
 739       pt->AddText("ALICE Preliminary");
 740       pt->Draw();
 741     }
 742
 743     gPad->SetGridx();
 744
 745   }
 746   else if (what == "PbPb_0010") {
 747     array = arrPbPb;
 748     system = 2; energy = 2760.; centrality = "V0M0010";
 749     label = "Pb-Pb #sqrt{s}_{NN} = 2.76 TeV, 0-10%";
 750     color = kRed+1;
 751     marker = kFullCircle;
 752     if(!shift)    shiftloc =  0;
 753     correlatedUnc = correlatedUncLocalPbPbOnlyKStarPhi;
 754
 755   }
 756
 757   else if (what == "PbPb_6080") {
 758     array = arrPbPb;
 759     system = 2; energy = 2760.; centrality = "V0M6080";
 760     label = "Pb-Pb #sqrt{s}_{NN} = 2.76 TeV, 60-80%";
 761     color  = kBlue+1;
 762     marker = kFullCircle;
 763     if(!shift)    shiftloc =  0.0;
 764     correlatedUnc = correlatedUncLocalPbPbOnlyKStarPhi;
 765   }
 766   else if (what == "PbPb_2040") {
 767     array = arrPbPb;
 768     system = 2; energy = 2760.; centrality = "V0M2040";
 769     label = "Pb-Pb #sqrt{s}_{NN} = 2.76 TeV, 20-40%";
 770     color  = kBlue+1;
 771     marker = kFullCircle;
 772     if(!shift)    shiftloc =  0.0;
 773     correlatedUnc = correlatedUncLocalPbPbOnlyKStarPhi;
 774   }
 775   else if (what == "pp7") {
 776     array = arrpp7;
 777     system = 0; energy = 7000.; centrality = "";
 778     label = "pp #sqrt{s} = 7 TeV";
 779     color  = kMagenta+1;
 780     marker = kFullCircle;
 781     if(!shift)    shiftloc =  0.2;
 782     //    correlatedUnc = correlatedUncLocalPP;
 783   }
 784   else if (what == "pp900") {
 785     array = arrpp900;
 786     system = 0; energy = 900.; centrality = "";
 787     label = "pp #sqrt{s} = 0.9 TeV";
 788     color  = kCyan+2;
 789     marker = kFullCircle;
 790     if(!shift)    shiftloc =  -0.2;
 791     //    correlatedUnc = correlatedUncLocalPP;
 792   }
 793   else if (what == "pp276") {
 794     array = arrpp276;
 795     system = 0; energy = 2760.; centrality = "";
 796     label = "pp #sqrt{s} = 2.76 TeV";
 797     color  = kYellow+2;
 798     marker = kFullCircle;
 799     if(!shift)    shiftloc = 0;
 800     //    correlatedUnc = correlatedUncLocalPP;
 801   }
 802   else if (what == "pPb0005") {
 803     array = arrpPb;
 804     system = 1; energy = 5020.; centrality = "V0A0005";
 805     //    label = "p-Pb #sqrt{s}_{NN} = 5.02 TeV, V0A 0-5%";
 806     label = "p-Pb #sqrt{s}_{NN} = 5.02 TeV, V0A Multiplicity (Pb-Side) 0-5%";
 807     color  = kBlack;
 808     marker = kFullCircle;
 809     if(!shift)    shiftloc = -0.2;
 810     //    correlatedUnc = correlatedUncLocalPP;
 811   }
 812   else if (what == "pPb2040") {
 813     array = arrpPb;
 814     system = 1; energy = 5020.; centrality = "V0A2040";
 815     label = "p-Pb #sqrt{s}_{NN} = 5.02 TeV, V0A Multiplicity (Pb-Side) 20-40%";
 816     color  = kBlue+1;
 817     marker = kFullCircle;
 818     if(!shift)    shiftloc = 0.;
 819     //    correlatedUnc = correlatedUncLocalPP;
 820   }
 821   else if (what == "pPb6080") {
 822     array = arrpPb;
 823     system = 1; energy = 5020.; centrality = "V0A6080";
 824     label = "p-Pb #sqrt{s}_{NN} = 5.02 TeV, V0A Multiplicity (Pb-Side) 60-80%";
 825     color  = kGreen+3;
 826     marker = kFullCircle;
 827     if(!shift)    shiftloc = 0.;
 828     //    correlatedUnc = correlatedUncLocalPP;
 829   }
 830   else if (what == "PbPbSTAR") {
 831     array = arrSTARPbPb;
 832     system = 2; energy = 200.; centrality = "0005";
 833     label = "STAR, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-5%";
 834     color  = kBlack;
 835     marker = kOpenStar;
 836     if(!shift)    shiftloc = +0.2;
 837     //    correlatedUnc = correlatedUncZero;
 838   }
 839   else if (what == "PbPbPHENIX") {
 840     array = arrPHENIXPbPb;
 841     system = 2; energy = 200.; centrality = "0005";
 842     label = "PHENIX, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-5%";
 843     color  = kBlack;
 844     marker = kOpenSquare;
 845     if(!shift)    shiftloc = -0.15;
 846     //    correlatedUnc = correlatedUncZero;
 847   }
 848   else if (what == "PbPbBRAHMS") {
 849     array = arrBRAHMSPbPb;
 850     system = 2; energy = 200.; centrality = "0010";
 851     label = "BRAHMS, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-10%";
 852     color  = kBlack;
 853     marker = kOpenCross;
 854     if(!shift)    shiftloc = -0.3;
 855     //    correlatedUnc = correlatedUncZero;
 856   }
 857   else if (what == "ppSTAR") {
 858     array = arrSTARpp;
 859     system = 0; energy = 200.; centrality = "";
 860     label = "STAR, pp #sqrt{s} = 0.2 TeV";
 861     color  = kBlack;
 862     marker = kOpenStar;
 863     if(!shift)    shiftloc = -0.15;
 864     //    correlatedUnc = correlatedUncZero;
 865   }
 866   else if (what == "ppPHENIX") {
 867     array = arrPHENIXpp;
 868     system = 0; energy = 200.; centrality = "";
 869     label = "PHENIX, pp #sqrt{s} = 0.2 TeV";
 870     color  = kBlack;
 871     marker = kOpenSquare;
 872     if(!shift)    shiftloc = -0.2;
 873     //    correlatedUnc = correlatedUncZero;
 874   }
 875   // From here on, it's meta names, to draw several series of ratios
 876   else if (what == "allpp"){
 877     DrawRatio("frame");
 878     DrawRatio("pp900");
 879     DrawRatio("pp276");
 880     DrawRatio("pp7");
 881     array =0;
 882     NewLegendQM(0.62249, 0.635734, 0.910643, 0.94673);
 883
 884     SaveCanvas("Ratios_pponly");
 885   }
 886   else if (what == "allppWithRHIC"){
 887     DrawRatio("frame");
 888     DrawRatio("pp900", 0, 0.05);
 889     DrawRatio("pp276", 0, 0.20);
 890     DrawRatio("pp7", 0, 0.35);
 891     DrawRatio("ppPHENIX", 0, -0.35);
 892     DrawRatio("ppSTAR", 0, -0.15);
 893
 894     array =0;
 895
 896     NewLegendQM(0.588353, 0.636857, 0.910643, 0.948352);
 897
 898     SaveCanvas("Ratios_pponly_withRHIC");
 899   }
 900   else if (what == "PbPbWithPP7TeV"){
 901
 902     DrawRatio("frame");
 903     DrawRatio("PbPb_0010", 0, 0.1);
 904     DrawRatio("pp7", 0, -0.1);
 905     array =0;
 906
 907
 908     NewLegendQM(0.538153, 0.749397, 0.893574, 0.950362);
 909     DrawExtrapolatedSymbolsAndLegendPbPb0010();
 910
 911     SaveCanvas("Ratios_withpp7tev");
 912
 913   } else if(what == "allsyst") {
 914
 915     DrawRatio("frame");
 916     DrawRatio("pp7", 0, -0.2);
 917     DrawRatio("pPb0005", 0, 0.00001);
 918     DrawRatio("PbPb_0010", 0, 0.2);
 919     array =0;
 920
 921
 922     NewLegendQM(0.462851, 0.631722, 0.89257, 0.936697);
 923     DrawExtrapolatedSymbolsAndLegendPbPb0010();
 924     DrawExtrapolatedSymbolsAndLegendpPb0005();
 925
 926     SaveCanvas("Ratios_allsystems");
 927
 928   }else if(what =="PbPb6080andpPb0005") {
 929     DrawRatio("frame");
 930     DrawRatio("pPb0005",0, -0.1);
 931     DrawRatio("PbPb_6080", 0, 0.1);
 932     array=0;
 933
 934     NewLegendQM(    0.46988, 0.730036, 0.910643, 0.948736);
 935     DrawExtrapolatedSymbolsAndLegendpPb0005();
 936     SaveCanvas("Ratios_6080vspPb");
 937
 938
 939   }else if(what =="pp_vsRHIC") {
 940     DrawRatio("frame");
 941     DrawRatio("pp7");
 942     DrawRatio("ppPHENIX", 0, 0.00001);
 943     DrawRatio("ppSTAR");
 944     array=0;
 945
 946     NewLegendQM(    0.554217, 0.677869, 0.910643, 0.948736);
 947
 948     SaveCanvas("Ratios_vsRHIC_pp");
 949
 950
 951   } else if (what =="PbPb_vsRHIC") {
 952     DrawRatio("frame");
 953     DrawRatio("PbPb_0010");
 954     DrawRatio("PbPbSTAR");
 955     DrawRatio("PbPbPHENIX");
 956     DrawRatio("PbPbBRAHMS");
 957     array = 0;
 958
 959
 960
 961     NewLegendQM(    0.434739, 0.591593, 0.939759, 0.936697);
 962     DrawExtrapolatedSymbolsAndLegendPbPb0010();
 963
 964     SaveCanvas("Ratios_vsRHIC_PbPb");
 965   } else if (what == "aliceall") {
 966     DrawRatio("frame");
 967     DrawRatio("PbPb_0010");
 968     DrawRatio("PbPb_6080");
 969     DrawRatio("pPb0005");
 970     DrawRatio("pp7");
 971     DrawRatio("pp276");
 972     DrawRatio("pp900");
 973   }
 974
 975   // FROM HERE: IT's yields
 976   else if( what == "fit_ReferenceFit_PbPb0010") {
 977
 978     DrawRatio("frame",1);
 979     DrawRatio("PbPb_0010",1);
 980     particlesToExcludeFromChi2="[313]"; // do not consider K*
 981     legThermal->SetNColumns(4);
 982     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
 983     // FIXME: sistemare valori rapporti
 984     shiftRatioDataModel=-0.2;
 985     PlotThermusYields("lhc2760_final_0005_single_gc_output_gs1_wevc_nkst.txt" , kBlack     , kSolid      , "THERMUS 2.3 , 155 #pm 2 , 5924 #pm 543 , 23.6/9");
 986     shiftRatioDataModel =0;
 987     PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat"                        , kOrange-1  , kDashed     , "GSI      , 156 #pm 2  , 5330 #pm 505 , 17.4/9");
 988     shiftRatioDataModel = 0.2;
 989     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_with_nuclei.txt"     , kBlue+1    , kDashDotted , "SHARE 3     , 156 #pm 3 , 4476 #pm 696 , 14.1/9");
 990
 991
 992     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
 993     DrawMarkerKStarNoFit() ;
 994     DrawExtrapolatedSymbolsYieldsPbPb0010();
 995
 996     SaveCanvas("Fit_PbPb0010_Reference");
 997
 998     array =0;
 999   }  else if( what == "fit_ReferenceFit_GSIONLY_PbPb0010") {
1000
1001     DrawRatio("frame",1);
1002     DrawRatio("PbPb_0010",1);
1003     particlesToExcludeFromChi2="[313]"; // do not consider K*
1004     legThermal->SetNColumns(4);
1005     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1006     // FIXME: sistemare valori rapporti
1007     shiftRatioDataModel =0;
1008     PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat"                        , kBlack  , kSolid     , "GSI      , 156 #pm 2  , 5330 #pm 505 , 17.4/9");
1009
1010     myPadHisto->cd();
1011     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1012     DrawMarkerKStarNoFit() ;
1013     DrawExtrapolatedSymbolsYieldsPbPb0010();
1014
1015     SaveCanvas("Fit_PbPb0010_Reference_GSI");
1016
1017     array =0;
1018   } else if( what == "fit_ReferenceFit_GSITHERMUS_PbPb0010") {
1019
1020     DrawRatio("frame",1);
1021     DrawRatio("PbPb_0010",1);
1022     particlesToExcludeFromChi2="[313]"; // do not consider K*
1023     legThermal->SetNColumns(4);
1024     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1025     // FIXME: sistemare valori rapporti
1026     shiftRatioDataModel =-0.1;
1027     PlotThermusYields("lhc2760_final_0005_single_gc_output_gs1_wevc_nkst.txt" , kBlack     , kSolid      , "THERMUS 2.3 , 155 #pm 2 , 5924 #pm 543 , 23.6/9");
1028     shiftRatioDataModel =0.1;
1029     PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat"                        , kOrange-1  , kDashed     , "GSI      , 156 #pm 2  , 5330 #pm 505 , 17.4/9");
1030
1031
1032     myPadHisto->cd();
1033     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1034     DrawMarkerKStarNoFit() ;
1035     DrawExtrapolatedSymbolsYieldsPbPb0010();
1036
1037     SaveCanvas("Fit_PbPb0010_Reference_GSITHERMUS");
1038     //    SaveCanvas("");
1039
1040     array =0;
1041   }else if( what == "fitSHARE_NoPionsNoProtons_PbPb0010") {
1042     array =0;
1043
1044     DrawRatio("frame",1);
1045     DrawRatio("PbPb_0010",1);
1046     legThermal->SetNColumns(4);
1047     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1048
1049     particlesToExcludeFromChi2 = "[313]";
1050     shiftRatioDataModel=-0.2;
1051     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_with_nuclei.txt"         ,kBlack    , kSolid      ,"SHARE 3          , 156 #pm 3 , 4476 #pm 696 , 14.1/9 ");
1052     particlesToExcludeFromChi2 = "[313][2212]";
1053     shiftRatioDataModel=0;
1054     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_with_nuclei_protons.txt" ,kOrange-1 , kDashed     ,"SHARE 3 (no p)   , 156 #pm 3 , 4520 #pm 623  , 8.2/8");
1055     particlesToExcludeFromChi2 = "[313][211]";
1056     shiftRatioDataModel=0.2;
1057     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_with_nuclei_pions.txt"   ,kBlue+1   , kDashDotted ,"SHARE 3 (no #pi) , 157 #pm 3 , 4103 #pm 677 , 12.2/8");
1058
1059     myPadHisto->cd();
1060     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1061
1062     // Add markers for additional particles not in fit
1063     DrawMarkerKStarNoFit() ;
1064     DrawExtrapolatedSymbolsYieldsPbPb0010();
1065
1066     SaveCanvas("Fit_PbPb0010_SHARE_NoPiNoP");
1067
1068
1069   } else if( what == "fitGSI_NoPionsNoProtons_PbPb0010") {
1070     array =0;
1071
1072     DrawRatio("frame",1);
1073     DrawRatio("PbPb_0010",1);
1074     legThermal->SetNColumns(4);
1075     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1076
1077     particlesToExcludeFromChi2 = "[313]";
1078     shiftRatioDataModel=-0.2;
1079     PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat"     , kBlack    , kSolid      , "GSI          , 156 #pm 1.5 , 5330 #pm 505  , 18.1/9" );
1080     particlesToExcludeFromChi2 = "[313][2212]";
1081     shiftRatioDataModel=0.;
1082     PlotGSIYields("data+therm_fit2_s2760_0-10qm14NOp.dat"  , kOrange-1 , kDashed     , "GSI (no p)   , 156 #pm 2   , 5590 #pm 330  ,  7.7/8" );
1083     particlesToExcludeFromChi2 = "[313][211]";
1084     shiftRatioDataModel=0.2;
1085     PlotGSIYields("data+therm_fit2_s2760_0-10qm14NOpi.dat" ,  kBlue+1  , kDashDotted , "GSI (no #pi) , 157 #pm 2   ,  4990 #pm 630 , 16.5/8");
1086
1087     myPadHisto->cd();
1088     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1089
1090     // Add markers for additional particles not in fit
1091     DrawMarkerKStarNoFit() ;
1092     DrawExtrapolatedSymbolsYieldsPbPb0010();
1093
1094     SaveCanvas("Fit_PbPb0010_GSI_NoPiNoP");
1095
1096
1097   } else if (what == "fitShareWithWithoutNuclei_PbPb0010") {
1098     DrawRatio("frame",1);
1099     DrawRatio("PbPb_0010",1);
1100     legThermal->SetNColumns(4);
1101     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1102
1103     shiftRatioDataModel = -0.1;
1104     particlesToExcludeFromChi2="[313]"; // do not consider K* and nuclei
1105     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_with_nuclei.txt"  , kBlack  , kSolid      , "SHARE 3             , 156 #pm 3 , 4476 #pm 696 , 14.1/9");
1106     shiftRatioDataModel = 0.1;
1107     particlesToExcludeFromChi2="[313][1000010020][1000020030][1010010030]"; // do not consider K* and nuclei
1108     PlotThermusYields("NEW_fit_gamma_q_fixed_PbPb_0010_without_nuclei.txt" , kCyan+2 , kDashDotted , "SHARE 3  (no nuclei), 156 #pm 4 , 4364 #pm 848 , 9.6/6");
1109
1110     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1111
1112     DrawMarkerKStarNoFit() ;
1113     DrawExtrapolatedSymbolsYieldsPbPb0010();
1114
1115     SaveCanvas("Fit_PbPb0010_SHARE_WithWoNuclei");
1116
1117
1118   } else if( what == "fitThermus_GammaSFree_PbPb0010") {
1119     std::cout << "MISSING DATA" << std::endl;
1120     array =0;
1121   } else if( what == "fitShare_GammaSGammaQFree_PbPb0010") {
1122     std::cout << "MISSING DATA" << std::endl;
1123     array =0;
1124   } else if( what == "fitShare_All_PbPb0010") {
1125     DrawRatio("frame",1);
1126     DrawRatio("PbPb_0010",1);
1127
1128
1129     legThermal->SetNColumns(6);
1130     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #gamma_{s}, #gamma_{q}, #chi^{2}/NDF", "0");
1131
1132     particlesToExcludeFromChi2="[313][1000010020][1000020030][1010010030]"; // do not consider K* and nuclei
1133     shiftRatioDataModel = -0.3;
1134     PlotThermusYields("NEW_fit_gamma_q_s_fixed_PbPb_0010_without_nuclei.txt" , kBlack      , kSolid      , "SHARE 3               , 156 #pm 4  , 4364 #pm 848  , 1 (fix)        , 1 (fix)        , 12.4/6");
1135     shiftRatioDataModel = -0.1;
1136     PlotThermusYields("NEW_fit_gamma_q_fixed_PbPb_0010_without_nuclei.txt"   , kRed+1      , kDashDotted , "SHARE 3               , 155 #pm 3  , 4406 #pm 766  , 1.07 #pm 0.05  , 1 (fix)        , 9.6/5");
1137     shiftRatioDataModel = 0.1;
1138     PlotThermusYields("NEW_fit_gamma_q_s_free_PbPb_0010_without_nuclei.txt"  , kCyan+2     , kDashed     , "SHARE 3               ,  138 #pm 6 , 3064 #pm 1319 ,  1.98 #pm 0.68 ,  1.63 #pm 0.38 , 3.1/4");
1139     shiftRatioDataModel = 0.3;
1140     particlesToExcludeFromChi2="[313]"; // do not consider K*
1141     PlotThermusYields("NEW_fit_gamma_q_s_free_PbPb_0010_with_nuclei.txt"     , kOrange - 1 , kDotted     , "SHARE 3 (with nuclei) , 152 #pm 8  , 4445 #pm 743  , 1.16 #pm 0.20  ,  1.06 #pm 0.12 ,  9.0/7");
1142
1143
1144     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1145     DrawMarkerKStarNoFit() ;
1146     DrawMarkerNucleiNoFit();
1147     DrawExtrapolatedSymbolsYieldsPbPb0010();
1148
1149
1150     legThermal->SetX1(0.12249  );
1151     legThermal->SetY1(0.0454769);
1152     legThermal->SetX2(0.821285 );
1153     legThermal->SetY2(0.383481 );
1154
1155     SaveCanvas("Fit_PbPb0010_SHARE_All");
1156     array =0;
1157   } else if( what == "fitGSI_PbPb6080") {
1158     DrawRatio("frame",1);
1159     DrawRatio("PbPb_6080",1);
1160     legThermal->SetNColumns(4);
1161     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1162     particlesToExcludeFromChi2 = "[313][1000020030][1010010030]";
1163     PlotGSIYields    ("data+therm_fit2_s2760_60-80qm14.dat", kBlack, kSolid, "GSI, 157 #pm 2, 210 #pm 20, 8.2/7");
1164
1165     legThermal->SetX1(0.143574  );
1166     legThermal->SetY1(0.0731318 );
1167     legThermal->SetX2(0.659639  );
1168     legThermal->SetY2(0.245206  );
1169
1170     myPadHisto->cd();
1171     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1172
1173     DrawMarkerKStarNoFit(1) ;
1174
1175     //    DrawExtrapolatedSymbolsYieldsPbPb0010(0.143574, 0.251352, 0.351406, 0.343535,0);
1176     //DrawExtrapolatedSymbolsYieldsPbPb0010();
1177
1178     SaveCanvas("Fit_PbPb6080_GSI");
1179     array =0;
1180   }else if( what == "fitGSI_PbPb2040") {
1181     DrawRatio("frame",1);
1182     DrawRatio("PbPb_2040",1);
1183     legThermal->SetNColumns(4);
1184     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1185     particlesToExcludeFromChi2 = "[313][1000020030][1010010030]";
1186     PlotGSIYields    ("data+therm_fit2_s2760_20-40qm14.dat", kBlack, kSolid, "GSI, 161 #pm 3, 1725 #pm 220, 24.6/7");
1187
1188     legThermal->SetX1(0.143574  );
1189     legThermal->SetY1(0.0731318 );
1190     legThermal->SetX2(0.659639  );
1191     legThermal->SetY2(0.245206  );
1192
1193     myPadHisto->cd();
1194     NewLegendQM(0.651606, 0.765993, 0.909639, 0.865951, 1);
1195
1196     DrawMarkerKStarNoFit(1) ;
1197     //    DrawMarkerKStarNoFit() ;
1198
1199     //    DrawExtrapolatedSymbolsYieldsPbPb0010(0.143574, 0.251352, 0.351406, 0.343535,0);
1200     //DrawExtrapolatedSymbolsYieldsPbPb0010();
1201
1202     SaveCanvas("Fit_PbPb2040_GSI");
1203     array =0;
1204   } else if( what == "fitGSI_pPb0005") {
1205     std::cout << "MISSING DATA" << std::endl;
1206     array =0;
1207   } else if( what == "fitGSI_pPb2040") {
1208     std::cout << "MISSING DATA" << std::endl;
1209     array =0;
1210   } else if( what == "fitThermus_GammaSFree_pPb0005") {
1211
1212     DrawRatio("frame",1);
1213     DrawRatio("pPb0005",1, 0.00001);
1214     particlesToExcludeFromChi2="[313][1000020030][1010010030] "; // do not consider K*
1215     legThermal->SetNColumns(6);
1216     AddLineToThermalLegend(legThermal, "Model,T (MeV), #gamma_{s}, r_{C} (fm), r (fm), #chi^{2}/NDF", "0");
1217
1218     shiftRatioDataModel=-0.2;
1219     PlotThermusYields("lhc5020_final_0005_single_gc_output_gs1_wevc_nkst.txt" , kBlack    , kSolid      , "THERMUS 2.3 GC , 158 #pm 2 , 1 (fix)       ,  N/A          , 3.40 #pm 0.11 , 32.0/7");
1220     shiftRatioDataModel =0;
1221     PlotThermusYields("lhc5020_final_0005_single_gc_output_gsf_wevc_nkst.txt" , kOrange-1 , kDashed     , "THERMUS 2.3 GC , 159 #pm 2 , 0.98 #pm 0.03 , N/A           , 3.40 #pm 0.11 , 31.3/6");
1222     shiftRatioDataModel = 0.2;
1223     PlotThermusYields("lhc5020_final_0005_single_sc_output_gs1_woevc_nkst.txt", kCyan+2   , kDashDotted , "THERMUS 2.3 SC , 158 #pm 3 , 1 (fix)       , 4.61 #pm 3.77 , 3.07 #pm 0.13  , 29.6/6");
1224    //lhc5020_final_0005_single_gc_output_gsf.txt
1225     //lhc5020_final_0005_single_sc_output_gs1.txt
1226
1227     NewLegendQM(0.650602, 0.694971, 0.909639, 0.865951, 1);
1228     DrawExtrapNotInFitpPb0005() ;
1229
1230     legThermal->SetX1(0.121486 );
1231     legThermal->SetY1(0.0741793);
1232     legThermal->SetX2(0.759036 );
1233     legThermal->SetY2(0.384575 );
1234
1235      SaveCanvas("Fit_pPb0005_THERMUS");
1236
1237     //DrawExtrapolatedSymbolsAndLegendPbPb0010();
1238
1239     array =0;
1240   } else if( what == "fitThermus_GammaSFree_pPb2040") {
1241     DrawRatio("frame",1);
1242     DrawRatio("pPb2040",1, 0.00001);
1243     particlesToExcludeFromChi2="[313][1000020030][1010010030] "; // do not consider K*
1244     legThermal->SetNColumns(6);
1245     AddLineToThermalLegend(legThermal, "Model,T (MeV), #gamma_{s}, r_{C} (fm), r (fm), #chi^{2}/NDF", "0");
1246
1247     shiftRatioDataModel=-0.2;
1248     PlotThermusYields("lhc5020_final_2040_single_gc_output_gs1_wevc_nkst.txt" , kBlack    , kSolid      , "THERMUS 2.3 GC , 155 #pm 2 , 1 (fix)       ,  N/A          , 2.83 #pm 0.08 , 40.1/7");
1249     shiftRatioDataModel =0;
1250     PlotThermusYields("lhc5020_final_2040_single_gc_output_gsf_wevc_nkst.txt" , kOrange-1 , kDashed     , "THERMUS 2.3 GC , 156 #pm 2 , 0.93 #pm 0.03 , N/A           , 2.83 #pm 0.08 , 34.6/6");
1251     shiftRatioDataModel = 0.2;
1252     PlotThermusYields("lhc5020_final_2040_single_sc_output_gs1_woevc_nkst.txt", kCyan+2   , kDashDotted , "THERMUS 2.3 SC , 156 #pm 2 , 1 (fix)       , 3.45 #pm 0.77 , 2.54 #pm 0.09 , 35.5/6");
1253    //lhc5020_final_0005_single_gc_output_gsf.txt
1254     //lhc5020_final_0005_single_sc_output_gs1.txt
1255
1256     NewLegendQM(0.650602, 0.694971, 0.909639, 0.865951, 1);
1257     DrawExtrapNotInFitpPb0005(0) ;
1258
1259     legThermal->SetX1(0.121486 );
1260     legThermal->SetY1(0.0741793);
1261     legThermal->SetX2(0.759036 );
1262     legThermal->SetY2(0.384575 );
1263
1264
1265
1266     SaveCanvas("Fit_pPb2040_THERMUS");
1267
1268     array =0;
1269   } else if( what == "fitThermus_GammaSFree_pPb6080") {
1270
1271     DrawRatio("frame",1);
1272     DrawRatio("pPb6080",1, 0.00001);
1273     particlesToExcludeFromChi2="[313][1000020030][1010010030][1000010020]"; // do not consider K*
1274     legThermal->SetNColumns(6);
1275     AddLineToThermalLegend(legThermal, "Model,T (MeV), #gamma_{s}, r_{C} (fm), r (fm), #chi^{2}/NDF", "0");
1276
1277     shiftRatioDataModel=-0.2;
1278
1279     PlotThermusYields("lhc5020_final_6080_single_gc_output_gs1_wevc_nkst.txt" , kBlack    , kSolid      , "THERMUS 2.3 GC , 152 #pm 2 , 1 (fix)    , N/A , 2.18 #pm 0.06 , 48.2/7");
1280     shiftRatioDataModel =0;
1281     PlotThermusYields("lhc5020_final_6080_single_gc_output_gsf_wevc_nkst.txt" , kOrange-1 , kDashed     , "THERMUS 2.3 GC , 154 #pm 2 , 0.88 #pm 3 , N/A , 2.21 #pm 0.07 , 28.8/6");
1282     shiftRatioDataModel = 0.2;
1283     PlotThermusYields("lhc5020_final_6080_single_sc_output_gs1_woevc_nkst.txt", kCyan+2   , kDashDotted , "THERMUS 2.3 SC , 154 #pm 2 ,1 (fix) ,3.18 #pm 0.57 ,1.96 #pm 0.07 ,40.5/6");
1284
1285
1286    //lhc5020_final_0005_single_gc_output_gsf.txt
1287     //lhc5020_final_0005_single_sc_output_gs1.txt
1288
1289     NewLegendQM(0.650602, 0.694971, 0.909639, 0.865951, 1);
1290     DrawExtrapNotInFitpPb0005(0) ;
1291
1292     legThermal->SetX1(0.121486  );
1293     legThermal->SetY1(0.0268308 );
1294     legThermal->SetX2(0.758032  );
1295     legThermal->SetY2(0.337226  );
1296
1297
1298     SaveCanvas("Fit_pPb6080_THERMUS");
1299
1300
1301
1302     array =0;
1303   }  else if( what == "fitThermus_RC_pPb0005") {
1304     std::cout << "MISSING DATA" << std::endl;
1305     array =0;
1306   } else if( what == "fitThermus_RC_pPb2040") {
1307     std::cout << "MISSING DATA" << std::endl;
1308     array =0;
1309   } else if( what == "fitThermus_RC_pPb6080") {
1310     std::cout << "MISSING DATA" << std::endl;
1311     array =0;
1312   } else if( what == "fitShare_pPb0005") {
1313
1314
1315     DrawRatio("frame",1);
1316     DrawRatio("pPb0005",1, 0.00001);
1317     legThermal->SetNColumns(6);
1318     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #gamma_{s}, #gamma_{q}, #chi^{2}/NDF", "0");
1319     particlesToExcludeFromChi2="[313][1000010020][1000020030][1010010030]"; // do not consider K* and nuclei
1320
1321     shiftRatioDataModel = -0.2;
1322     PlotThermusYields("NEW_fit_gamma_q_s_fixed_pPb_0005_without_nuclei.txt" , kBlack    , kSolid      , "SHARE 3 , 158 #pm 3 , 121 #pm 18  , 1 (fix)       , 1 (fix) , 24.3/6");
1323     shiftRatioDataModel = 0;
1324     PlotThermusYields("NEW_fit_gamma_q_fixed_pPb_0005_without_nuclei.txt"   , kBlue+1   , kDashDotted , "SHARE 3 , 161 #pm 3 , 115 #pm 17  , 0.93 #pm 0.04 , 1 (fix) , 20.3/5");
1325     shiftRatioDataModel = 0.2;
1326     PlotThermusYields("NEW_fit_gamma_q_s_free_pPb_0005_without_nuclei.txt"  , kOrange+2 , kDashDotted , "SHARE 3 , 144 #pm 1 , 81 #pm 25   ,  1.599 #pm 0 , 1.71 #pm 0.06 , 11.4/4");
1327
1328     NewLegendQM(0.650602, 0.694971, 0.909639, 0.865951, 1);
1329     DrawExtrapNotInFitpPb0005() ;
1330     DrawMarkerNucleiNoFit();
1331     legThermal->SetX1(0.124498 );
1332     legThermal->SetY1(0.0715488);
1333     legThermal->SetX2(0.672691 );
1334     legThermal->SetY2(0.384575 );
1335
1336     SaveCanvas("Fit_pPb0005_SHARE");
1337
1338
1339     array =0;
1340
1341   } else if( what == "fitShare_pPb0005_NoOmega") {
1342
1343
1344     DrawRatio("frame",1);
1345     DrawRatio("pPb0005",1, 0.00001);
1346     legThermal->SetNColumns(6);
1347     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #gamma_{s}, #gamma_{q}, #chi^{2}/NDF", "0");
1348     particlesToExcludeFromChi2="[313][1000010020][1000020030][1010010030]"; // do not consider K* and nuclei
1349
1350     shiftRatioDataModel = -0.3;
1351     PlotThermusYields("NEW_fit_gamma_q_s_fixed_pPb_0005_without_nuclei.txt"                 , kBlack    , kSolid      , "SHARE 3                 , 158 #pm 3  , 121 #pm 18  , 1 (fix)       , 1 (fix), 24.3/6");
1352     shiftRatioDataModel = -.1;
1353     PlotThermusYields("NEW_fit_gamma_q_s_fixed_pPb_0005_without_nuclei_excluding_Omega.txt" , kBlue+1   , kDashDotted , "SHARE 3 (No #Omega)     , 163 #pm 4  ,  99 #pm 17  , 1 (fix)       , 1 (fix), 15.6/5");
1354     shiftRatioDataModel = 0.1;
1355     PlotThermusYields("NEW_fit_gamma_q_fixed_pPb_0005_without_nuclei_excluding_Omega.txt"   , kOrange+2 , kDashed     , "SHARE 3 (No #Omega)     , 163 #pm 3  , 100 #pm 16  , 0.96 #pm 0.04 , 1 (fix), 14.7/4");
1356     shiftRatioDataModel = 0.3;
1357     PlotThermusYields("NEW_fit_gamma_q_s_fixed_pPb_0005_with_nuclei_excluding_Omega.txt"    , kCyan+2   , kDotted     , "SHARE 3 (No #Omega + d) , 160 #pm 3  , 114 #pm 17  , 1 (fix)       , 1 (fix), 18.6/6");
1358
1359     NewLegendQM(0.650602, 0.694971, 0.909639, 0.865951, 1);
1360     DrawExtrapNotInFitpPb0005() ;
1361     DrawMarkerNucleiNoFit();
1362     legThermal->SetX1(0.124498 );
1363     legThermal->SetY1(0.0715488);
1364     legThermal->SetX2(0.672691 );
1365     legThermal->SetY2(0.384575 );
1366
1367     SaveCanvas("Fit_pPb0005_SHARE_NoOmega");
1368
1369
1370     array =0;
1371
1372   }  else if( what == "fitShare_pPb2040") {
1373     std::cout << "MISSING DATA" << std::endl;
1374     array =0;
1375   } else if( what == "fitShare_pPb6080") {
1376     std::cout << "MISSING DATA" << std::endl;
1377     array =0;
1378   } else if( what == "fitGSI_pp") {
1379     DrawRatio("frame",1);
1380     DrawRatio("pp7",1, 0.00001);
1381     legThermal->SetNColumns(5);
1382     AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #gamma_{s}, #chi^{2}/NDF", "0");
1383     particlesToExcludeFromChi2 = "[313][1000020030][1010010030]";
1384     shiftRatioDataModel=-0.1;
1385     PlotGSIYields    ("data+therm_fit2_s7000.dat", kBlack, kSolid, "GSI, 146 #pm 2, 25 #pm 2, 1 (fix),  78.2/7", 0);
1386     shiftRatioDataModel=0.1;
1387     PlotGSIYields    ("data+therm_fit2_s7000gs.dat", kRed+1, kDashed, "GSI, 150 #pm 2, 23 #pm 2, 0.88, 45.6/7", 0);
1388     legThermal->SetX1(0.143574 );
1389     legThermal->SetY1(0.0731318);
1390     legThermal->SetX2(0.659639 );
1391     legThermal->SetY2(0.374263 );
1392     myPadHisto->cd();
1393
1394     NewLegendQM(0.725904, 0.761431, 0.874498, 0.841323, 1);
1395
1396     DrawMarkerKStarNoFit();
1397     TLegend * leg  = new TLegend(0.143574, 0.380408, 0.351406, 0.4603, NULL, "brNDC");
1398     leg->SetBorderSize(0);
1399     leg->SetTextFont(43);
1400     leg->SetTextSize(14);
1401     leg->SetLineColor(1);
1402     leg->SetLineStyle(1);
1403     leg->SetLineWidth(1);
1404     leg->SetFillColor(0);
1405     leg->SetFillStyle(1001);
1406     TLegendEntry * entry=leg->AddEntry("TMarker","Not in fit","p");
1407     entry->SetLineColor(1);
1408     entry->SetLineStyle(1);
1409     entry->SetLineWidth(1);
1410     entry->SetMarkerColor(1);
1411     entry->SetMarkerStyle(markerNoFit);
1412     entry->SetMarkerSize(1.2);
1413     entry->SetTextFont(43);
1414     leg->Draw();
1415
1416     SaveCanvas("Fit_pp7000_GSI");
1417     //DrawExtrapNotInFitpPb0005(0);
1418     array =0;
1419   }else if( what == "fitFlorence_pp") {
1420     DrawRatio("frame",1);
1421     DrawRatio("pp7",1, 0.00001);
1422     legThermal->SetNColumns(1);
1423     AddLineToThermalLegend(legThermal, "Model", "0");
1424     shiftRatioDataModel = -0.1;
1425     particlesToExcludeFromChi2 = "[1000020030][1010010030][1000010020]";
1426     PlotFlorenceYields("pp7000_Florence.txt", kBlack, kSolid, "Becattini et al.");
1427     shiftRatioDataModel = 0.1;
1428     particlesToExcludeFromChi2 = "[313][1000020030][1010010030]";
1429     PlotGSIYields    ("data+therm_fit2_s7000gs.dat", kRed+1, kDashed, "GSI GC",0);
1430     myPadHisto->cd();
1431     NewLegendQM(0.725904, 0.761431, 0.874498, 0.841323, 1);
1432
1433     SaveCanvas("Fit_pp7000_Florence");
1434     //DrawExtrapNotInFitpPb0005(0);
1435     array =0;
1436   } else if( what == "fitGSI_fullCanonical") {
1437     std::cout << "MISSING DATA" << std::endl;
1438     array =0;
1439   } else if( what == "fitThermus_RC_pp") {
1440     std::cout << "MISSING DATA" << std::endl;
1441     array =0;
1442   } else if( what == "fitShare_pp") {
1443     std::cout << "MISSING DATA" << std::endl;
1444     array =0;
1445   }
1446
1447
1448
1449   else {
1450     std::cout << "Unknown ratio " << what.Data() << std::endl;
1451     return;
1452   }
1453
1454   if(!correlatedUnc) {
1455     std::cout << "correlatedUnc not set!" << std::endl;
1456
1457   }
1458   std::cout << "CORR: " << correlatedUnc[1] << std::endl;
1459
1460   if(array) {
1461     if(isYield) {
1462       TGraphErrors * hstat =
1463       AliPWGHistoTools::GetGraphFromHisto(GetHistoYields(array,  system,  energy, centrality, label, color, marker, kStatError, shiftloc)
1464                                           ,0);
1465       hstat->SetMarkerSize(1.2);
1466       hstat->Draw("PZ");
1467       AliPWGHistoTools::GetGraphFromHisto(GetHistoYields(array,  system,  energy, centrality, label+"NoLegend", color, marker, kSystError, shiftloc)
1468                                           ,0)->Draw("[]");
1469     } else {
1470       AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array,  system,  energy, centrality, label, color, marker, kStatError, shiftloc)
1471                                           ,0)->Draw("PZ");
1472       AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array,  system,  energy, centrality, label+"NoLegend", color, marker, kSystError, shiftloc)
1473                                           ,0)->Draw("[]");
1474     }
1475   }
1476
1477
1478 }
1479
1480 void DrawFrame(Bool_t isYield) {
1481
1482   myCan = new TCanvas("myCan","ThermalFits",50,10,1000,isYield ? 950 : 650);
1483   myCan->Draw();
1484   myCan->cd();
1485   // Set the Pads
1486   Double_t boundaryLabels = isYield ? 0.88 : 0.85;//0.92;
1487   myPadLabel = new TPad("myPadLabel","myPadLabel",0.0, boundaryLabels,1.0,1.0);
1488   myPadSetUp(myPadLabel);
1489   myPadLabel->Draw();
1490
1491   myPadHisto = new TPad("myPadHisto","myPadHisto",0.0,isYield ? 0.4 : 0.05 ,1.0, boundaryLabels,0);
1492   myPadSetUp(myPadHisto);
1493   myPadHisto->Draw();
1494
1495   myPadStdDev = new TPad("myPadStdDev","myPadStdDev",0.0,0.035,1.0,0.215,0);
1496   myPadSetUp(myPadStdDev);
1497   if(isYield)  myPadStdDev->Draw();
1498   myPadStdDev->SetGridx();
1499   myPadStdDev->SetGridy();
1500
1501   // This pad is for the ratios data/model
1502   myPadRatio = new TPad("myPadRatio","myPadRatio",0.0,0.22,1.0,0.4,0);
1503   myPadSetUp(myPadRatio);
1504   if(isYield)  myPadRatio->Draw();
1505   myPadRatio->SetGridx();
1506   myPadRatio->SetGridy();
1507
1508
1509
1510   myPadLabel->cd();
1511
1512   double xLabelPosition[nratio] = {0.124498, 0.211847, 0.31, 0.38, 0.465, 0.575, 0.644, 0.72, 0.82, 0.905 };
1513   double xLabelYields[npart]    = {0.115, 0.185, 0.270, 0.330, 0.422, 0.485, 0.570, 0.625, 0.695, 0.785, 0.835, 0.915};
1514   double yLabelPosition     = 0.40;
1515
1516   // labels
1517   if(isYield) {
1518     for(Int_t ipart = 0; ipart < npart; ipart++){
1519       TLatex *myPart = new TLatex(xLabelYields[ipart],yLabelPosition, yieldsLabels[ipart]);
1520       myLatexDraw(myPart,20);
1521     }
1522   }  else {
1523     for(Int_t iratio = 0; iratio < nratio; iratio++){
1524       TLatex *myRatio = new TLatex(xLabelPosition[iratio],yLabelPosition, ratiosLabels[iratio]);
1525       myLatexDraw(myRatio,20);
1526     }
1527   }
1528   // Xi's and Omega's bar (there was no way to convince root to draw it properly)
1529   if(isYield) {
1530     TLine *line = new TLine(0.653,0.660,0.663,0.660);
1531     line->SetLineWidth(2);
1532     line->Draw();
1533     line = new TLine(0.725,0.660,0.738,0.660);
1534     line->SetLineWidth(2);
1535     line->Draw();
1536     line = new TLine(0.362,0.660,0.374,0.660);
1537     line->SetLineWidth(2);
1538     line->Draw();
1539
1540
1541   }
1542   else {
1543     TLine *line = new TLine(0.408,0.68,0.418,0.68);
1544     line->SetLineWidth(2);
1545     line->Draw();
1546     line = new TLine(0.498,0.68,0.513,0.68);
1547     line->SetLineWidth(2);
1548     line->Draw();
1549
1550
1551   }
1552
1553
1554   if(isYield) {
1555     myPadStdDev->cd();
1556     myPadStdDev->SetTopMargin(0.0);
1557     myPadStdDev->SetTicks(1,1);
1558
1559     Float_t devMax = 4.7;
1560
1561     TH2F *myBlankStdDev = new TH2F("myBlankStdDev","myBlankStdDev",npart,1,npart+1,10,-devMax,+devMax);
1562     myBlankStdDev->GetXaxis()->SetLabelFont(43); // precision 3: size will be in pixels
1563     myBlankStdDev->GetYaxis()->SetLabelFont(43);
1564     myBlankStdDev->GetYaxis()->SetTitleFont(43);
1565     myBlankStdDev->SetLabelSize(23,"xy");
1566     myBlankStdDev->SetTitleSize(20,"y");
1567     myBlankStdDev->SetNdivisions(20,"x");
1568     myBlankStdDev->SetNdivisions(605,"y");
1569     myBlankStdDev->SetLabelOffset(0.012,"xy");
1570     myBlankStdDev->SetYTitle("(mod.-data)/#sigma_{data}");
1571     myBlankStdDev->SetTitleOffset(1.65,"y");
1572     myBlankStdDev->Draw();
1573
1574     TH2F *myBlankRatio = new TH2F("myBlankRatio","myBlankRatio",npart,1,npart+1,10,-0.8,0.8);
1575     myBlankRatio->GetXaxis()->SetLabelFont(43); // precision 3: size will be in pixels
1576     myBlankRatio->GetYaxis()->SetLabelFont(43);
1577     myBlankRatio->GetYaxis()->SetTitleFont(43);
1578     myBlankRatio->SetLabelSize(23,"xy");
1579     myBlankRatio->SetTitleSize(20,"y");
1580     myBlankRatio->SetNdivisions(20,"x");
1581     myBlankRatio->SetNdivisions(605,"y");
1582     myBlankRatio->SetLabelOffset(0.012,"xy");
1583     myBlankRatio->SetYTitle("(mod.-data)/mod.");
1584     myBlankRatio->SetTitleOffset(1.65,"y");
1585     myPadRatio->cd();
1586     //    myPadRatio->SetLogy();
1587     myBlankRatio->Draw();
1588
1589   }
1590
1591   myPadHisto->cd();
1592   myPadHisto->SetBottomMargin(isYield ? 0.01 : 0.04);
1593   //  myPadHisto->SetLogy();
1594   myPadHisto->SetTicks(1,1);
1595   TH1 *myBlankHisto = 0;
1596   if(isYield) {
1597     myBlankHisto =  new TH2F("NoLegend","NoLegend",npart,1,npart+1,10,  0.00002, maxy+0.01  );
1598     gPad->SetLogy();
1599
1600   } else {
1601    myBlankHisto =  new TH2F("NoLegend","NoLegend",nratio,1,nratio+1,10,  0, maxy+0.01  );
1602   }
1603   myBlankHisto->GetXaxis()->SetLabelFont(43); // precision 3: size will be in pixels
1604   myBlankHisto->GetYaxis()->SetLabelFont(43);
1605   myBlankHisto->GetYaxis()->SetTitleFont(43);
1606   myBlankHisto->SetLabelSize(23,"y");
1607   myBlankHisto->SetLabelSize(0,"x");
1608   myBlankHisto->SetTitleSize(26,"y");
1609   myBlankHisto->SetMaximum(10);
1610   myBlankHisto->SetMinimum(0);
1611   myBlankHisto->SetNdivisions(isYield? 20 :10,"x");
1612   myBlankHisto->SetNdivisions(505,"y");
1613   if(isYield) myBlankHisto->SetYTitle("d#it{N}/d#it{y}");
1614   myBlankHisto->SetLabelOffset(0.012,"xy");
1615   myBlankHisto->SetTitleOffset(isYield ? 1.3 : 1,"y");
1616   myBlankHisto->Draw();
1617
1618   if(isYield) {
1619     legThermal = new TLegend(0.144578, 0.0702247, 0.659639, 0.383226);
1620     legThermal->SetBorderSize(1);
1621     legThermal->SetTextFont(43);
1622     legThermal->SetTextSize(18);
1623     //    legThermal->SetNColumns(6);
1624     legThermal->SetFillColor(0);
1625     legThermal->SetLineWidth(1);
1626     legThermal->Draw();
1627   }
1628 }
1629
1630 void myLatexDraw(TLatex *currentLatex, Float_t currentSize, Int_t currentColor){
1631   currentLatex->SetTextFont(43);
1632   currentLatex->SetTextSizePixels(Int_t(currentSize));
1633   //  currentLatex->SetTextAngle(0);
1634   currentLatex->Draw();
1635   return;
1636 }
1637
1638 void myPaveSetup(float rRatio, float rRange3, float rRange5,
1639                  int rFillColor){
1640
1641   float cHiRange = 0, cLoRange = 0;
1642   if(rRange3<rRange5) {cHiRange=rRange5;cLoRange=rRange3;}
1643   else {cHiRange=rRange3;cLoRange=rRange5;}
1644
1645   TPave *cPave= new TPave(rRatio-0.25,cLoRange,rRatio+0.25,cHiRange,0,"br");
1646   cPave->SetFillColor(rFillColor);
1647   cPave->SetLineColor(1);
1648   cPave->Draw();
1649 }
1650
1651 void myPadSetUp(TPad *currentPad){
1652   currentPad->SetLeftMargin(0.10);
1653   currentPad->SetRightMargin(0.02);
1654   currentPad->SetTopMargin(0.02);
1655   currentPad->SetBottomMargin(0.02);
1656   return;
1657 }
1658 TGraphErrors*  PlotThermusYields(const char * filename, Int_t color, Int_t lineStyle,
1659                                  const char * tag) {
1660
1661   Int_t lw = lineStyle == kSolid ? 2 : 3; // Set line width
1662
1663
1664   std::map<Int_t,Double_t> mapYields;
1665   std::map<Int_t,Double_t> mapStdDev;
1666
1667   Int_t pdg;
1668   Double_t yield, stddev;
1669   ifstream thermusFile(filename);
1670   TString line;
1671   Double_t chi2 = 0;
1672   std::cout << "---"<<tag<<"---" << std::endl;
1673
1674   //  std::istream is(thermusFile);
1675   // Read the std dev and the ratio in 2 maps, then plot them in a graph.
1676   while(line.ReadLine(thermusFile, kTRUE)) {
1677     if(line.BeginsWith("#")) continue;
1678     TObjArray * tokens = line.Tokenize(" \t");
1679     if(tokens->GetEntries() != 3) continue;// not a line with data
1680     //    thermusFile >> pdg >> yield >> stddev;
1681     pdg    = ((TObjString*)tokens->At(0))->String().Atof();
1682     yield  = ((TObjString*)tokens->At(1))->String().Atof();
1683     stddev = ((TObjString*)tokens->At(2))->String().Atof();
1684
1685     if( thermusFile.eof() ) break;
1686     std::cout << "PDG " << pdg << " " << yield << " " << stddev << std::endl;
1687
1688     mapYields[TMath::Abs(pdg)] += yield;
1689     mapStdDev[TMath::Abs(pdg)] += stddev;
1690     if(pdg < 0) {
1691       // found the antiparticle: now compute the mean
1692       mapYields[TMath::Abs(pdg)] /=2;
1693       mapStdDev[TMath::Abs(pdg)] /=2;
1694     }
1695     delete tokens;
1696   }
1697
1698   // Now plot
1699   TGraphErrors * gThermus = new TGraphErrors;
1700   TGraphErrors * gThermusStdDev = new TGraphErrors;
1701   for(Int_t ipart = 0; ipart < npart; ipart++){
1702     gThermus->SetPoint(ipart, ipart+1.5, mapYields[particleYields[ipart]]);
1703     gThermus->SetPointError(ipart, 0.3, 0);
1704     gThermusStdDev->SetPoint(ipart, ipart+1.5, mapStdDev[particleYields[ipart]]);
1705     gThermusStdDev->SetPointError(ipart, 0.3, 0);
1706   }
1707
1708   myPadHisto->cd();
1709   gThermus->Draw("PZ");
1710   gThermus->SetLineWidth(lw);
1711
1712   gThermus->SetLineColor(color);
1713   gThermus->SetLineStyle(lineStyle);
1714   gThermus->SetTitle("NoLegend");
1715
1716
1717   //  myPadStdDev->cd();
1718   // gThermusStdDev->Draw("PZ");
1719   // gThermusStdDev->SetLineWidth(lw);
1720   // gThermusStdDev->SetLineColor(color);
1721   // gThermusStdDev->SetLineStyle(lineStyle);
1722   TGraphErrors* gStdDev2 = 0;
1723   TGraphErrors* gRatio   = 0;
1724   std::cout << "CHI2: "
1725             <<  GetGraphRatioAndStdDev(gThermus, gRatio, gStdDev2)
1726             << std::endl;
1727   myPadRatio->cd();
1728   gRatio->Draw("PZ");
1729   myPadStdDev->cd();
1730   gStdDev2->Draw("PZ");
1731   myPadHisto->cd();
1732   AddLineToThermalLegend(gThermus, tag, "l");
1733   return gThermus;
1734 }
1735
1736 TGraphErrors*  PlotFlorenceYields(const char * filename, Int_t color, Int_t lineStyle,
1737                                  const char * tag) {
1738
1739   Int_t lw = lineStyle == kSolid ? 2 : 3; // Set line width
1740
1741
1742   std::map<Int_t,Double_t> mapYields;
1743   std::map<Int_t,Double_t> mapStdDev;
1744
1745   Int_t pdg;
1746   Double_t yield, yieldbar, stddev;
1747   ifstream thermusFile(filename);
1748   TString line;
1749   std::cout << "---"<<tag<<"---" << std::endl;
1750
1751   //  std::istream is(thermusFile);
1752   // Read the std dev and the ratio in 2 maps, then plot them in a graph.
1753   while(line.ReadLine(thermusFile, kTRUE)) {
1754     if(line.BeginsWith("#")) continue;
1755     TObjArray * tokens = line.Tokenize(" \t");
1756     if(tokens->GetEntries() != 4) continue;// not a line with data
1757     //    thermusFile >> pdg >> yield >> stddev;
1758     pdg      = ((TObjString*)tokens->At(0))->String().Atof();
1759     yield    = ((TObjString*)tokens->At(2))->String().Atof();
1760     yieldbar = ((TObjString*)tokens->At(3))->String().Atof(); // Antiparticle yield
1761     if(pdg == 0) {
1762       // not a line with data
1763       delete tokens;
1764       continue ;
1765     }
1766     if( thermusFile.eof() ) break;
1767     std::cout << "PDG " << pdg << " " << yield << " " << yieldbar << std::endl;
1768
1769     mapYields[TMath::Abs(pdg)] += yieldbar ? ((yield+yieldbar)/2) : yield; // If the antiparticle exists, use the average
1770
1771     delete tokens;
1772   }
1773
1774   // Now plot
1775   TGraphErrors * gFlorence = new TGraphErrors;
1776   for(Int_t ipart = 0; ipart < npart; ipart++){
1777     gFlorence->SetPoint(ipart, ipart+1.5, mapYields[particleYields[ipart]]);
1778     gFlorence->SetPointError(ipart, 0.3, 0);
1779   }
1780
1781   myPadHisto->cd();
1782   gFlorence->Draw("PZ");
1783   gFlorence->SetLineWidth(lw);
1784
1785   gFlorence->SetLineColor(color);
1786   gFlorence->SetLineStyle(lineStyle);
1787   gFlorence->SetTitle("NoLegend");
1788
1789
1790
1791   TGraphErrors* gStdDev2 = 0;
1792   TGraphErrors* gRatio   = 0;
1793   std::cout << "CHI2: "
1794             <<  GetGraphRatioAndStdDev(gFlorence, gRatio, gStdDev2)
1795             << std::endl;
1796   myPadRatio->cd();
1797   gRatio->Draw("PZ");
1798   myPadStdDev->cd();
1799   gStdDev2->Draw("PZ");
1800   myPadHisto->cd();
1801   AddLineToThermalLegend(gFlorence, tag, "l");
1802   return gFlorence;
1803 }
1804
1805
1806 TGraphErrors*  PlotGSIYields(const char * filename, Int_t color, Int_t lineStyle,
1807                              const char * tag, Bool_t isPbPb) {
1808
1809   // tag is a comma separated list of elements to be added to the legend as diferent columns
1810
1811   Int_t lw = lineStyle == kSolid ? 2 : 3; // Set line width
1812
1813
1814   const Int_t pdgPbPb0010[] = {211, -211, 321, -321, 310, 313, 333, 2212, -2212, 3122, 3312, -3312, 3334, -3334, 1000010020, 1000020030, 1010010030, -1010010030};
1815   const Int_t pdgpp[] = {211, 321, 310, 313, 333, 2212, 3122, 3312, 3334, 1000010020};
1816
1817   const Int_t *pdgOrder = isPbPb ? pdgPbPb0010 : pdgpp;
1818
1819   std::map<Int_t,Double_t> mapYields;
1820   std::map<Int_t,Double_t> mapStdDev;
1821   std::map<Int_t,Double_t> mapUncert;
1822   std::map<Int_t,Double_t> mapData;
1823
1824   Double_t data, uncert, model;
1825
1826   ifstream gsiFile(filename);
1827   //  std::istream is(thermusFile);
1828   std::cout << "----- "<<tag<<" -----" << std::endl;
1829
1830   // Read the std dev and the ratio in 2 maps, then plot them in a graph.
1831   Int_t ipart = 0;
1832   while(gsiFile) {
1833     gsiFile >> data >> uncert >> model;
1834     if( gsiFile.eof() ) break;
1835     Int_t pdg = pdgOrder[ipart];
1836     std::cout << "PDG " << pdg << " " << data << std::endl;
1837     mapYields[TMath::Abs(pdg)] += model;
1838     mapUncert[TMath::Abs(pdg)] += uncert;
1839     mapData[TMath::Abs(pdg)]   += data;
1840     if(pdg < 0) {
1841       // found the antiparticle: now compute the mean
1842       mapYields[TMath::Abs(pdg)] /=2;
1843       mapData[TMath::Abs(pdg)] /=2;
1844       mapUncert[TMath::Abs(pdg)] /=2;
1845
1846     }
1847     ipart++;
1848   }
1849
1850   // Now plot
1851   TGraphErrors * gGsi = new TGraphErrors;
1852   TGraphErrors * gGsiStdDev = new TGraphErrors;
1853   std::cout << "PDG    \tmodel\tdata\tuncert\tstddev" << std::endl;  // header
1854   for(Int_t ipart = 0; ipart < npart; ipart++){ // Here we use npart, as this is what we wnat to plot!
1855     Int_t pdg = particleYields[ipart];
1856     mapStdDev[TMath::Abs(pdg)]  = ( mapYields[TMath::Abs(pdg)] - mapData[TMath::Abs(pdg)]) /  mapUncert[TMath::Abs(pdg)] ;
1857     std::cout << "PDG " << pdg <<"\t"
1858               << mapYields[TMath::Abs(pdg)] << "\t" << mapData[TMath::Abs(pdg)] <<"\t"
1859               << mapUncert[TMath::Abs(pdg)] << "\t" << mapStdDev[TMath::Abs(pdg)]
1860               << std::endl;
1861
1862     if(!mapYields[particleYields[ipart]]) mapYields[particleYields[ipart]] = -10;
1863     if(!mapStdDev[particleYields[ipart]]) mapStdDev[particleYields[ipart]] = -10;
1864
1865
1866     gGsi->SetPoint(ipart, ipart+1.5, mapYields[particleYields[ipart]]);
1867     gGsi->SetPointError(ipart, 0.3, 0);
1868     gGsiStdDev->SetPoint(ipart, ipart+1.5, mapStdDev[particleYields[ipart]]);
1869     gGsiStdDev->SetPointError(ipart, 0.3, 0);
1870   }
1871
1872
1873
1874   myPadHisto->cd();
1875   gGsi->Draw("PZ");
1876   gGsi->SetLineWidth(lw);
1877
1878   gGsi->SetLineColor(color);
1879   gGsi->SetLineStyle(lineStyle);
1880   gGsi->SetTitle("NoLegend");
1881
1882
1883   // myPadStdDev->cd();
1884   // gGsiStdDev->Draw("PZ");
1885   // gGsiStdDev->SetLineWidth(lw);
1886   // gGsiStdDev->SetLineColor(color);
1887   // gGsiStdDev->SetLineStyle(lineStyle);
1888   // myPadHisto->cd();
1889
1890   TGraphErrors* gStdDev2 = 0;
1891   TGraphErrors* gRatio   = 0;
1892   std::cout << "CHI2: "
1893             <<  GetGraphRatioAndStdDev(gGsi, gRatio, gStdDev2)
1894             << std::endl;
1895   myPadRatio->cd();
1896   gRatio->Draw("PZ");
1897   myPadStdDev->cd();
1898   gStdDev2->Draw("PZ");
1899
1900   AddLineToThermalLegend(gGsi, tag, "L");
1901
1902   return gGsi;
1903 }
1904
1905
1906 void DrawExtrapolatedSymbolsAndLegendPbPb0010() {
1907     myPadHisto->cd();
1908
1909
1910     TLegend *leg = new TLegend(    0.149598, 0.782203, 0.415663, 0.858447,NULL,"pNDC");
1911     leg->SetBorderSize(0);
1912     leg->SetTextFont(43);
1913     leg->SetTextSize(18);
1914     leg->SetLineColor(1);
1915     leg->SetLineStyle(1);
1916     leg->SetLineWidth(2);
1917     leg->SetFillColor(0);
1918     leg->SetFillStyle(1001);
1919     leg->SetMargin(0.1);
1920
1921     TLegendEntry *entry=leg->AddEntry("TMarker","Extrapolated (Pb-Pb 0-10%)","p");
1922     entry->SetLineColor(1);
1923     entry->SetLineStyle(1);
1924     entry->SetLineWidth(1);
1925     entry->SetMarkerColor(1);
1926     entry->SetMarkerStyle(27);
1927     entry->SetMarkerSize(1.2);
1928     entry->SetTextFont(43);
1929     leg->Draw();
1930     myPadLabel->cd();
1931     // Markers for extrapolated points
1932
1933     TMarker *marker = new TMarker(0.666667,0.111825,27);
1934     marker->SetMarkerStyle(27);
1935     marker->SetMarkerSize(1.2);
1936     marker->Draw();
1937     marker = new TMarker(0.920683,0.111825,27);
1938     marker->SetMarkerStyle(27);
1939     marker->SetMarkerSize(1.2);
1940     marker->Draw();
1941
1942     // BR for 3He
1943
1944     myPadLabel->cd();
1945     TLatex *   tex = new TLatex(0.73,0.05,"BR = 25%");
1946     tex->SetNDC();
1947     tex->SetTextFont(43);
1948     tex->SetTextSize(14);
1949     tex->SetLineWidth(2);
1950     tex->Draw();
1951
1952
1953
1954 }
1955 void DrawExtrapolatedSymbolsAndLegendpPb0005() {
1956     myPadHisto->cd();
1957
1958
1959     TLegend *leg = new TLegend(    0.149598, 0.709972, 0.415663, 0.786216,NULL,"pNDC");
1960     leg->SetBorderSize(0);
1961     leg->SetTextFont(43);
1962     leg->SetTextSize(18);
1963     leg->SetLineColor(1);
1964     leg->SetLineStyle(1);
1965     leg->SetLineWidth(2);
1966     leg->SetFillColor(0);
1967     leg->SetFillStyle(1001);
1968     leg->SetMargin(0.1);
1969     TLegendEntry *entry=leg->AddEntry("TMarker","Extrapolated (p-Pb 0-5%)","p");
1970     entry->SetLineColor(1);
1971     entry->SetLineStyle(1);
1972     entry->SetLineWidth(1);
1973     entry->SetMarkerColor(1);
1974     entry->SetMarkerStyle(28);
1975     entry->SetMarkerSize(1.2);
1976     entry->SetTextFont(43);
1977     leg->Draw();
1978     myPadLabel->cd();
1979
1980    TMarker *marker = new TMarker(0.590361,0.111825,28);
1981    marker->SetMarkerStyle(28);
1982    marker->SetMarkerSize(1.2);
1983    marker->Draw();
1984    marker = new TMarker(0.938755,0.111825,28);
1985    marker->SetMarkerStyle(28);
1986    marker->SetMarkerSize(1.2);
1987    marker->Draw();
1988 }
1989
1990 void AddLineToThermalLegend(TObject * obj, TString line, const char * optFirst) {
1991
1992   // This should be a comma-separated list of text to be added to the
1993   // columns.  If the number of entries does not match the numer of
1994   // columns, it returns an error
1995   TObjArray * tokens = line.Tokenize(",");
1996   if(tokens->GetEntries() != legThermal->GetNColumns()) {
1997     std::cout << "Wrong number of columns (" << tokens->GetEntries() << ","<<legThermal->GetNColumns()<<") for the thermal legend, not adding " << line.Data() << std::endl;
1998     return;
1999   }
2000
2001   TIter iter(tokens);
2002   TObjString * col;
2003   Bool_t first = 1;
2004   while((col = (TObjString*)iter.Next())) {
2005     // Add entry removing whitespaces
2006     legThermal->AddEntry(obj, col->String().Strip(TString::kBoth, ' ').Data(), first ? optFirst : "0");
2007     if (first) first = 0;
2008   }
2009 }
2010
2011 void DrawExtrapolatedSymbolsYieldsPbPb0010(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Bool_t plotExtraploatedLegend){
2012     // Markers for extrapolated points
2013   myPadLabel->cd();
2014   TMarker * marker = new TMarker(0.36245,0.111825,markerExtrap);
2015   marker->SetMarkerStyle(markerExtrap);
2016   marker->SetMarkerSize(1.2);
2017   marker->Draw();
2018   marker = new TMarker(0.945783,0.111825,markerExtrap);
2019   marker->SetMarkerStyle(markerExtrap);
2020   marker->SetMarkerSize(1.2);
2021   marker->Draw();
2022   // BR for 3He
2023   myPadHisto->cd();
2024   TLatex *   tex = new TLatex(11.15, 2.5e-5,"BR = 25%");
2025   tex->SetTextFont(43);
2026   tex->SetTextSize(14);
2027   tex->SetLineWidth(2);
2028   tex->Draw();
2029
2030   myPadHisto->cd();
2031
2032
2033
2034     TLegend * leg = new TLegend(x1,y1,x2,y2,NULL,"brNDC");
2035   leg->SetBorderSize(0);
2036   leg->SetTextFont(43);
2037   leg->SetTextSize(14);
2038   leg->SetLineColor(1);
2039   leg->SetLineStyle(1);
2040   leg->SetLineWidth(1);
2041   leg->SetFillColor(0);
2042   leg->SetFillStyle(1001);
2043   TLegendEntry * entry=leg->AddEntry("TMarker","Not in fit","p");
2044   entry->SetLineColor(1);
2045   entry->SetLineStyle(1);
2046   entry->SetLineWidth(1);
2047   entry->SetMarkerColor(1);
2048   entry->SetMarkerStyle(markerNoFit);
2049   entry->SetMarkerSize(1.2);
2050   entry->SetTextFont(43);
2051   if(plotExtraploatedLegend) {
2052     entry=leg->AddEntry("TMarker","Extrapolated","p");
2053     entry->SetLineColor(1);
2054     entry->SetLineStyle(1);
2055     entry->SetLineWidth(1);
2056     entry->SetMarkerColor(1);
2057     entry->SetMarkerStyle(markerExtrap);
2058     entry->SetMarkerSize(1.2);
2059     entry->SetTextFont(43);
2060   }
2061   leg->Draw();
2062   leg->SetMargin(0.1);
2063
2064
2065 }
2066
2067 void DrawMarkerKStarNoFit(Bool_t plotLegend) {
2068   myPadLabel->cd();
2069   TMarker *marker = new TMarker(0.344378,0.111825,markerNoFit);
2070   marker->SetMarkerStyle(markerNoFit);
2071   marker->SetMarkerSize(1.2);
2072   marker->Draw();
2073   myPadHisto->cd();
2074
2075   if(plotLegend) {
2076
2077
2078     TLegend * leg = new TLegend(0.126506, 0.253051, 0.335341, 0.345118,NULL,"brNDC");
2079     leg->SetBorderSize(0);
2080     leg->SetTextFont(43);
2081     leg->SetTextSize(14);
2082     leg->SetLineColor(1);
2083     leg->SetLineStyle(1);
2084     leg->SetLineWidth(1);
2085     leg->SetFillColor(0);
2086     leg->SetFillStyle(1001);
2087     TLegendEntry * entry=leg->AddEntry("TMarker","Not in fit","p");
2088     entry->SetLineColor(1);
2089     entry->SetLineStyle(1);
2090     entry->SetLineWidth(1);
2091     entry->SetMarkerColor(1);
2092     entry->SetMarkerStyle(markerNoFit);
2093     entry->SetMarkerSize(1.2);
2094     entry->SetTextFont(43);
2095     leg->Draw();
2096   }
2097 }
2098
2099 void DrawMarkerNucleiNoFit() {
2100   myPadLabel->cd();
2101   TMarker *marker = new TMarker(0.928715,0.111825,markerNoFit);
2102   marker->SetMarkerStyle(markerNoFit);
2103   marker->SetMarkerSize(1.2);
2104   marker->Draw();
2105
2106   marker = new TMarker(0.791751,0.111825,markerNoFit);
2107   marker->SetMarkerStyle(markerNoFit);
2108   marker->SetMarkerSize(1.2);
2109   marker->Draw();
2110
2111   marker = new TMarker(0.866466,0.111825,markerNoFit);
2112   marker->SetMarkerStyle(markerNoFit);
2113   marker->SetMarkerSize(1.2);
2114   marker->Draw();
2115   myPadHisto->cd();
2116
2117 }
2118
2119 Double_t GetGraphRatioAndStdDev(TGraphErrors * gModel, TGraphErrors * &gRatio, TGraphErrors *&gStdDev) {
2120
2121   // Plots the ratios data/model
2122
2123   // I recomputed the stddev and the chi2 here, because some values
2124   // changed slightly due to typos in some of the input files with
2125   // respect to the ones used for those fits.
2126
2127
2128   if(!gModel) {
2129     std::cout << "EMPTY MODEL" << std::endl;
2130     return 0;
2131   }
2132   // 0. set up the graphs which we will need
2133   TGraphErrors * gTemp = 0;
2134   TGraphErrors * gStat = 0;
2135   TGraphErrors * gSyst = 0;
2136   // the cloning below takes care of the style
2137   gRatio = (TGraphErrors*)gModel->Clone();
2138   gRatio->Clear();
2139   gStdDev = (TGraphErrors*)gModel->Clone();
2140   gStdDev->Clear();
2141
2142   // 1. Find the data graphs. We need both the stat and the syst, since we want the ratio with the total uncertainty (or not?)
2143   TVirtualPad * currentPad = gPad;
2144   myPadHisto->cd();
2145   TList * list = myPadHisto->GetListOfPrimitives();
2146   TIterator * iter = list->MakeIterator();
2147
2148   while ((gTemp = (TGraphErrors*) iter->Next())){
2149     if(gTemp->InheritsFrom("TGraphErrors")) {
2150       // Found a graph, it is the correct one?
2151       TString name = gTemp->GetName();
2152       TString title = gTemp->GetTitle();
2153       std::cout << "name " << name.Data() << std::endl;
2154
2155       if (name.Contains("hPart")) {
2156         // ok, it's the data
2157         if (title.Contains("NoLegend")) gSyst = gTemp; // it's the syst error FIXME: it would be better to add the error type to the data
2158         else gStat = gTemp;
2159       }
2160       if(gStat && gSyst) {
2161         std::cout << "Cannot find gStat or gSyst (" << gStat << "," << gSyst <<")" << std::endl;
2162         break; // found both stat and syst uncertainties
2163       }
2164     }
2165   }
2166
2167
2168
2169   // Compute the ratio, the stddev and the chi2
2170   Int_t npoint = gModel->GetN(); // We are sure that data and model have the same number of points in the same order, because they were created using the particleYields array (FIXME: is this also true for GSI?)
2171   Double_t chi2 = 0;
2172   for(Int_t ipoint = 0; ipoint < npoint; ipoint++){
2173
2174     Double_t yield = gStat->GetY()[ipoint];
2175     Double_t stat  = gStat->GetEY()[ipoint];
2176     Double_t syst  = gSyst->GetEY()[ipoint];
2177     Double_t error = TMath::Sqrt(stat*stat+syst*syst);
2178     Double_t model = gModel->GetY()[ipoint];
2179     Double_t width = gModel->GetEX()[ipoint];
2180     Double_t ratio  = (model-yield)/model;
2181     Double_t stddev = (model-yield)/error;
2182     if(!particlesToExcludeFromChi2.Contains(Form("[%d]", particleYields[ipoint]))) chi2 += TMath::Power(stddev,2);
2183     else std::cout << "Ecluding PDG "<< particleYields[ipoint] <<" from chi2 calculation" << std::endl;
2184     Double_t errorRatio = error/model;
2185     gRatio->SetPoint(ipoint, gModel->GetX()[ipoint]+shiftRatioDataModel, ratio);
2186     gRatio->SetPointError(ipoint, 0, errorRatio);
2187     if(model) {
2188       gStdDev->SetPoint(ipoint, gModel->GetX()[ipoint], stddev);
2189       gStdDev->SetPointError(ipoint, width, 0);
2190     }
2191     // The commented block down here was used to compare to the estimate provided directly by Boris and Benjamin
2192     // gStdDev->SetPoint(ipoint, gModel->GetX()[ipoint]-0.4, stddev);
2193     // gStdDev->SetPointError(ipoint, 0.2, 0);
2194     // std::cout << "PDG " << particleYields[ipoint] <<"\t"
2195     //        << model << "\t" << yield <<"\t"
2196     //        << error << "\t" << stddev
2197     //        << std::endl;
2198   }
2199
2200   gRatio->SetLineStyle(kSolid);
2201   gRatio->SetMarkerStyle(kOpenSquare);
2202   gRatio->SetMarkerColor(gRatio->GetLineColor());
2203   gRatio->SetLineWidth(2);
2204
2205   currentPad->cd();
2206   return chi2;
2207 }
2208
2209
2210 void SaveCanvas(const char * name) {
2211
2212   if(!saveCanvas) return;
2213
2214   std::cout << "Saving " << name << ".{eps,pdf,root,C}" << std::endl;
2215
2216   myCan->Update();
2217   gSystem->ProcessEvents();
2218   myCan->Print(Form("%s.eps",name));
2219   myCan->Print(Form("%s.root",name));
2220   myCan->Print(Form("%s.C",name));
2221   gSystem->Exec(Form("epstopdf %s.eps", name));
2222   gSystem->Exec(Form("if [ \"$USER\" = \"mfloris\" ]; then cp %s.{eps,pdf,root,C} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ",name));
2223
2224 }
2225
2226 void DrawExtrapNotInFitpPb0005(Bool_t drawExtrap) {
2227   myPadLabel->cd();
2228   TMarker *marker = new TMarker(0.344378,0.111825,28);
2229   marker->SetMarkerStyle(28);
2230   marker->SetMarkerSize(1.2);
2231   marker->Draw();
2232   // marker = new TMarker(0.7851406,0.111825,28);
2233   // marker->SetMarkerStyle(28);
2234   // marker->SetMarkerSize(1.2);
2235   // marker->Draw();
2236   if(drawExtrap) {
2237     marker = new TMarker(0.364458,0.111825,27);
2238     marker->SetMarkerStyle(27);
2239     marker->SetMarkerSize(1.2);
2240     marker->Draw();
2241     marker = new TMarker(0.792,0.111825,27);
2242     marker->SetMarkerStyle(27);
2243     marker->SetMarkerSize(1.2);
2244     marker->Draw();
2245   }
2246   myPadHisto->cd();
2247
2248
2249
2250   TLegend * leg = 0;
2251
2252   if (drawExtrap) leg = new TLegend(0.123494, 0.400358, 0.331325, 0.534512,NULL,"brNDC");
2253   else leg = new TLegend(0.123494, 0.395097, 0.332329, 0.474011, NULL, "brNDC");
2254   leg->SetBorderSize(0);
2255   leg->SetTextFont(43);
2256   leg->SetTextSize(14);
2257   leg->SetLineColor(1);
2258   leg->SetLineStyle(1);
2259   leg->SetLineWidth(1);
2260   leg->SetFillColor(0);
2261   leg->SetFillStyle(1001);
2262   TLegendEntry * entry=leg->AddEntry("TMarker","Not in fit","p");
2263   entry->SetLineColor(1);
2264   entry->SetLineStyle(1);
2265   entry->SetLineWidth(1);
2266   entry->SetMarkerColor(1);
2267   entry->SetMarkerStyle(markerNoFit);
2268   entry->SetMarkerSize(1.2);
2269   entry->SetTextFont(43);
2270
2271   if(drawExtrap) {
2272     entry=leg->AddEntry("TMarker","Extrapolated","p");
2273     entry->SetLineColor(1);
2274     entry->SetLineStyle(1);
2275     entry->SetLineWidth(1);
2276     entry->SetMarkerColor(1);
2277     entry->SetMarkerStyle(markerExtrap);
2278     entry->SetMarkerSize(1.2);
2279     entry->SetTextFont(43);
2280   }
2281   leg->Draw();
2282   leg->SetMargin(0.1);
2283 }