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