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