]> 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 interpolation for the 20-40% bin
[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
31
9d84731e 32enum 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};
33
046fa53e 34typedef enum {kStatError, kSystError, kTotalError} myerror_t;
35
36TH1F * 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 37TH1F * 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 38void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
046fa53e 39void SetStyle(Bool_t graypalette=0) ;
b926af7d 40TLegend * NewLegendQM(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Bool_t isYield = 0) ;
41void DrawRatio(TString what, Bool_t isYield = kFALSE, Double_t shiftloc=0.);
15fccaa0 42
69f3eeda 43void DrawFrame(Bool_t yields = 0) ;
b926af7d 44void DrawExtrapolatedSymbolsAndLegendPbPb0010() ;
45void DrawExtrapolatedSymbolsAndLegendpPb0005() ;
4910214c 46void DrawMarkerKStarNoFit() ;
47void DrawExtrapolatedSymbolsYieldsPbPb0010();
48
b926af7d 49
046fa53e 50void LoadArrays() ;
69f3eeda 51//void AddLabel(Float_t x, Float_t y, TString text);
52void myLatexDraw(TLatex *currentLatex, Float_t currentSize=0.5, Int_t currentColor=1);
53void myPaveSetup(float rRatio=0, float rRange3=0, float rRange5=0,
54 int rFillColor=0);
55void myPadSetUp(TPad *currentPad);
b926af7d 56TGraphErrors* PlotThermusYields(const char * filename, Int_t color, Int_t lineStyle,
4910214c 57 const char * tag);
b926af7d 58TGraphErrors * PlotGSIYields(const char * fileName, Int_t color=kBlack, Int_t lineStyle = kSolid,
4910214c 59 const char * tag ="");
60
61void AddLineToThermalLegend(TObject * obj, TString line, const char * optFirst = "L");
69f3eeda 62
63// Ratios to be draw. Remember to change the labels in DrawFrame if you change this
64const Int_t nratio = 10;
65// Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGPi , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK};
66
67Int_t num [nratio] = {kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar};
68Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGProton , kPDGDeuteron , kPDGPi , kPDGK , kPDGK};
15fccaa0 69Int_t isSum[nratio] = {1 , 1 , 1 , 1 , 1 , 1 , 0 , 1 , 1 , 1 };
15fccaa0 70const char * ratiosLabels[] = {"#frac{K^{+}+K^{-}}{#pi^{+}+#pi^{-}}",
71 "#frac{p+#bar{p}}{#pi^{+}+#pi^{-}}",
b926af7d 72 "#frac{2#Lambda}{K_{S}^{0}}",
15fccaa0 73 "#frac{#Xi^{-}+#Xi^{+}}{#pi^{+}+#pi^{-}}",
74 "#frac{#Omega^{-}+#Omega^{+}}{#pi^{+}+#pi^{-}}",
75 "#frac{d}{p+#bar{p}}",
76 "#frac{{}^{3}He }{d}",
77 "#frac{{}^{3}_{#Lambda}H+{}^{3}_{#Lambda}#bar{H} }{#pi^{+}+#pi^{-}}",
78 "#frac{#phi}{K^{+}+K^{-}}",
79 "#frac{K*+#bar{K}*}{K^{+}+K^{-}}",};
69f3eeda 80static const Double_t scale[] = {1 , 3 , 0.5 , 30 , 250 , 50 , 100 , 4e5 , 2 , 1 };
81//static const Double_t scale[] = {1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,};
b926af7d 82const Int_t npart = 12;
83Int_t particleYields [npart] = {kPDGPi ,kPDGK ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3 };
84Int_t isSumYields[npart] = {1 ,1 ,0 , 1 , 0 , 1 ,1 ,1 ,1 ,0 , 1 , 0 };
85//Int_t isSumInputFiles[npart] = {1 ,1 ,0 , 1 , 0 , 1 ,1 ,1 ,1 ,0 , 1 , 0 };
86const char * yieldsLabels[] = {"#frac{#pi^{+}+#pi^{-}}{2}",
87 "#frac{K^{+}+K^{-}}{2}",
88 "K_{S}^{0}",
89 "#frac{K*+#bar{K}*}{2}",
90 "#phi",
91 "#frac{p+#bar{p}}{2}",
92 "#Lambda",
93 "#frac{#Xi^{-}+#Xi^{+}}{2}",
94 "#frac{#Omega^{-}+#Omega^{+}}{2}",
95 "d",
96 "#frac{{}^{3}_{#Lambda}H+{}^{3}_{#Lambda}#bar{H}}{2}",
97 "{}^{3}He",
98};
9d84731e 99
046fa53e 100// Preferred colors and markers
046fa53e 101const Int_t markers[] = {kFullCircle, kFullSquare,kOpenCircle,kOpenSquare,kOpenDiamond,kOpenCross,kFullCross,kFullDiamond,kFullStar,kOpenStar,0};
69f3eeda 102//
b926af7d 103Int_t markerNoFit = 28;
104Int_t markerExtrap = 27;
15fccaa0 105Double_t maxy = 0.5;
9d84731e 106
046fa53e 107// Data arrays;
108TClonesArray *arrPbPb=0, *arrpp7=0, *arrpPb=0, * arrpp276=0, * arrpp900=0, * arrThermus=0;
109TClonesArray *arrSTARPbPb=0, *arrPHENIXPbPb=0, *arrBRAHMSPbPb=0;
110TClonesArray *arrSTARpp =0, *arrPHENIXpp=0;
111
69f3eeda 112//const Double_t *scaleRatios = 0;
4c7bd9d1 113Double_t *correlatedUnc = 0;
114Double_t correlatedUncLocalPbPb[14] = {0.0424 , 0.0424 , 0.041 , 0 , 0 , 0.0424 , 0.0424 , 0 , 0.05 , 0.05 };
115Double_t correlatedUncLocalPP [14] = {0.0424 , 0.0424 , 0.041 , 0 , 0 , 0.0424 , 0.0424 , 0 , 0.0424 , 0.0424 };
116Double_t correlatedUncZero[14] = {0,0,0,0,0,0,0,0,0,0,0,0,0,0};
117
69f3eeda 118TCanvas *myCan = 0;
b926af7d 119TPad *myPadStdDev =0;
120TPad *myPadHisto =0;
121TPad *myPadLabel =0;
122TLegend * legThermal = 0;
4c7bd9d1 123
9d84731e 124TClonesArray * PlotRatiosForQM14() {
ad431d97 125#if !(!defined (__CINT__) || (defined(__MAKECINT__)))
046fa53e 126 LoadLibs();
ad431d97 127#endif
9d84731e 128
046fa53e 129 //
130 LoadArrays();
ad431d97 131
046fa53e 132 // Uncomment stuff in this section to save the inputs for thermal models
4910214c 133#define SAVE_INPUT_THERMAL_MODEL
046fa53e 134#ifdef SAVE_INPUT_THERMAL_MODEL
b926af7d 135 PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
4910214c 136 // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/1);
137 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", /*separateCharges*/1);
138 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A2040", /*separateCharges*/1);
139 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A6080", /*separateCharges*/1);
140 // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M6080", /* separateCharges*/1);
141 // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M2040", /*separateCharges*/1);
69f3eeda 142
b926af7d 143 PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
4910214c 144 // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0);
145 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", /*separateCharges*/0);
146 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A2040", /*separateCharges*/0);
147 // PrepareThermalModelsInputFiles(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A6080", /*separateCharges*/0);
148 // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M6080", /*separateCharges*/0);
149 // PrepareThermalModelsInputFiles(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M2040", /*separateCharges*/0);
9d84731e 150
046fa53e 151 return 0;
152#endif
9d84731e 153
046fa53e 154 SetStyle();
ad431d97 155
b926af7d 156 //DrawRatio("allpp");
15fccaa0 157 // DrawRatio("PbPbWithPP7TeV");
4910214c 158 // DrawRatio("allsyst");
15fccaa0 159 // DrawRatio("PbPb6080andpPb0005");
b926af7d 160 //DrawRatio("pp_vsRHIC");
4910214c 161 // DrawRatio("PbPb_vsRHIC");
15fccaa0 162 // DrawRatio("aliceall");
ad431d97 163
b926af7d 164
165 // Yields and FITS
4910214c 166 maxy=20000;
167 //DrawRatio("fit_ReferenceFit_PbPb0010", 1);
168 // DrawRatio("fitSHARE_NoPionsNoProtons_PbPb0010",1);
169 // DrawRatio("fitGSI_NoPionsNoProtons_PbPb0010", 1);
b926af7d 170 // DrawRatio("fitShare_pPb0005", 1);
171 //DrawRatio("fitShare_All_PbPb0010", 1);
4910214c 172 // DrawRatio("fitShareWithWithoutNuclei_PbPb0010", 1);
173 // DrawRatio("fitGSI_PbPb6080", 1);
15fccaa0 174 // NewLegendQM();
046fa53e 175 return arrPbPb;
9d84731e 176}
177
046fa53e 178TH1F * 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 179 // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY
9d84731e 180
4c7bd9d1 181
69f3eeda 182 // scaleRatios = scale;
046fa53e 183 TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s_%d",system,energy,centrality.Data(),errorType), histotitle, nratio, 1+shift, nratio+1+shift);
184
185 TClonesArray * arrRatios = new TClonesArray ("AliParticleYield");// We save to disk the precomputed ratios
186 Int_t iratioArr = 0;// Index used only to add particles to the array above
9d84731e 187
188 // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
189 for(Int_t iratio = 1; iratio <= nratio; iratio++){
79936afc 190 std::cout << "******** " << num[iratio-1] << "/" << denum[iratio-1]<< " ("<<isSum[iratio-1]<<")*******" << std::endl ;
9d84731e 191
79936afc 192 AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
193 if(ratio)
194 {
195 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)
196 std::cout << " " ;
197 ratio->Print("short");
198 }
199
9d84731e 200 if(!ratio) {
201 // If the ratio is not found, try to build it!
79936afc 202 std::cout << " Looking for " << num[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
9d84731e 203 AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality, isSum[iratio-1]);
69f3eeda 204 if(part1) {
205 part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
206 if(isSum[iratio-1] && part1->IsTypeAverage()) {
207 std::cout << "Sum requested, found average, scaling x2" << std::endl;
208 part1->Scale(2.);
209 }
210 }
ad431d97 211 // Try with the !sum, if part 1 is not found
212 if(!part1) {
79936afc 213 std::cout << " Looking for " << num[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
9d84731e 214 part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
79936afc 215 if(part1)
216 {
217 part1 = new AliParticleYield(*part1); // We need to clone it to avoid a mess if we need to use this particle again later
218 // If the sum was requested, try to recover it!
4c7bd9d1 219 if(isSum[iratio-1]) {
79936afc 220 std::cout << " Looking for " << -num[iratio-1] <<std::endl;
15fccaa0 221 // If it's lambda and ALICE, use 2L instead of L + Lbar // FIXME: do the same for deuterons?
222 if((num[iratio-1] == kPDGLambda || num[iratio-1] == kPDGDeuteron) && energy > 300) {
223 std::cout << " It's Lambda or deuteron ALICE: Scaling x2 " << std::endl;
224 part1->Print();
225 part1->Scale(2.);
226 } else {
227 AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0);
228 if(part1 && part1_bar) {
229 std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl;
230 part1 = AliParticleYield::Add(part1, part1_bar);
231
232 } else if (TDatabasePDG::Instance()->GetParticle(-num[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
233 std::cout << "WARNING: Sum requested but not found, scaling x2 " << part1->GetName() << std::endl;
234 part1->Scale(2);
235 }
79936afc 236 }
237 } else if(part1) { // if we are here, it means the sum was *not* requested (isSum=0), but we found something with (!isSum) = 1
238 // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
239 std::cout << "WARNING: Using sum/2 for " << part1->GetName() << std::endl;
046fa53e 240
79936afc 241 part1->Scale(0.5);
242 }
ad431d97 243 }
244
9d84731e 245 }
15fccaa0 246
247
79936afc 248 std::cout << " Looking for " << denum[iratio-1] << " ("<<isSum[iratio-1]<<")"<<std::endl;
9d84731e 249 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
69f3eeda 250 if(part2) {
251 part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
252 if(isSum[iratio-1] && part2->IsTypeAverage()) {
253 std::cout << "Sum requested, found average, scaling x2" << std::endl;
254 part2->Scale(2.);
255 }
256 }
9d84731e 257 if(!part2) {// Try with the !sum, if part 2 is not found
79936afc 258 std::cout << " Looking for " << denum[iratio-1] << " ("<<!isSum[iratio-1]<<")"<<std::endl;
9d84731e 259 part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
79936afc 260 if(part2)
261 {
262 part2 = new AliParticleYield(*part2); // We need to clone it to avoid a mess if we need to use this particle again later
4c7bd9d1 263 if(isSum[iratio-1]) {
79936afc 264 std::cout << " Looking for " << -denum[iratio-1] << std::endl;
265 AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0);
266 if(part2 && part2_bar){
267 std::cout << "Adding " << part2_bar->GetPartName() << " " << part2->GetPartName() << std::endl;
268 part2 = AliParticleYield::Add(part2, part2_bar);
4c7bd9d1 269 } else if (TDatabasePDG::Instance()->GetParticle(-denum[iratio-1])){ // Before scaling x2 check it it makes sense (antiparticle exists)
79936afc 270 std::cout << "WARNING: Sum requested but not found, scaling x2 " << part2->GetName() << std::endl;
271 part2->Scale(2);
272 }
273 } else if(part2){
274 // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
275 std::cout << "WARNING: Using sum/2 for " << part2->GetName() << std::endl;
276 part2->Scale(0.5);
277 }
278 }
ad431d97 279
9d84731e 280 }
4c7bd9d1 281 ratio = AliParticleYield::Divide(part1, part2, correlatedUnc[iratio-1], "YQ"); // Assume by that the systematics of part1 and part2 are uncorrelated.
69f3eeda 282
9d84731e 283 if(ratio) {
ad431d97 284 std::cout << "" << std::endl;
9d84731e 285 std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
ad431d97 286 ratio->Print("");
287 std::cout << "" << std::endl;
9d84731e 288 }
289 }
290 if(ratio){
ad431d97 291 ratio->Scale(scale[iratio-1]);
9d84731e 292 h->SetBinContent(iratio, ratio->GetYield());
046fa53e 293 if(errorType == kTotalError) {
294 h->SetBinError (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
295 } else if (errorType == kStatError) {
296 h->SetBinError (iratio, ratio->GetStatError());
297 } else if (errorType == kSystError) {
298 h->SetBinError (iratio, ratio->GetSystError());
299 } else {
300 std::cout << "ERROR: Unknown Error Type " << errorType << std::endl;
301 }
302
303 // h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f", scale[iratio-1]), ratio->GetLatexName()));
304 h->GetXaxis()->SetBinLabel(iratio, ratio->GetLatexName());
305 new ((*arrRatios)[iratioArr++]) AliParticleYield(*ratio);
9d84731e 306 }
307 else {
308 h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
309
310 }
79936afc 311 std::cout << "*** END OF " << num[iratio-1] << "/" << denum[iratio-1]<< " *******" << std::endl ;
312
9d84731e 313 }
314
046fa53e 315 h->GetYaxis()->SetRangeUser(0, maxy);
316 h->SetLineColor(icolor);
317 h->SetMarkerColor(icolor);
318 h->SetMarkerStyle(imarker);
69f3eeda 319 // the "if" avoids saving twice the same ratios
320 if(errorType == kSystError) AliParticleYield::SaveAsASCIIFile(arrRatios, TString("ratios_")+h->GetName());
9d84731e 321 return h;
322
323
324
325}
326
327void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges) {
328 // If "Separate charges" is true, tries to dump both charges are dumped
329 TClonesArray * arrOut = new TClonesArray("AliParticleYield");
69f3eeda 330 TClonesArray * arrOutGSI = new TClonesArray("AliParticleYield"); // We add dummy lines to the GSI output file if needed!
f5b6253e 331;
9d84731e 332
333 Int_t ipartOut = 0; // Index for the array
69f3eeda 334 Int_t ipartOutGSI = 0; // Index for the array
335
9d84731e 336 for(Int_t ipart = 0; ipart < npart; ipart++){
337 if(!separateCharges) {
f5b6253e 338 AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality, isSumYields[ipart]);
339 if(!part && isSumYields[ipart]) {
9d84731e 340 //Could not find the particle, but the sum was requested: build the sum!
f5b6253e 341 part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality, 0);
342 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particleYields[ipart], system, energy, centrality, 0);
69f3eeda 343 if(part2 && part) part = AliParticleYield::Add(part, part2);
344 else if(part) part->Scale(2.); // If we only found a particle, we can scale it by a factor 2.
9d84731e 345 else part = 0;
346 }
69f3eeda 347 // We want to save the average of particle and antiparticle in this case
348 if(part) {
f5b6253e 349 if(isSumYields[ipart] && !part->IsTypeAverage()) part->Scale(0.5); // If it's not already an average, but just a sum, divide by 2
69f3eeda 350 new((*arrOut )[ipartOut++]) AliParticleYield(*part);
351 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
352 } else { // Add dummy particle to the GSI list
f5b6253e 353 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
69f3eeda 354 }
9d84731e 355 }
356 else {
f5b6253e 357 // ignore isSumYields and try to find both particleYields
9d84731e 358 Bool_t notFound = 0;
f5b6253e 359 AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart], system, energy, centrality, 0);
360 if(part) {
361 new((*arrOut)[ipartOut++]) AliParticleYield(*part);
362 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
363 }
364 else {
365 // std::cout << "ADDING DUMMY part " << particleYields[ipart] << std::endl;
366 // new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
367 notFound=1;
368 }
9d84731e 369 // Try to find antiparticle (-pdg code)
f5b6253e 370 part = AliParticleYield::FindParticle(arr, -particleYields[ipart], system, energy, centrality, 0);
69f3eeda 371 if(part) {
372 new((*arrOut)[ipartOut++]) AliParticleYield(*part);
373 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
374 }
9d84731e 375 else if (notFound) {
376 // If neither charge was found, check if we at least have the sum
f5b6253e 377 part = AliParticleYield::FindParticle(arr, abs(particleYields[ipart]), system, energy, centrality, 1);
378 if (part) {
379 if(!part->IsTypeAverage()) part->Scale(0.5); // If it's a sum (not an average) divide by 2
69f3eeda 380 new((*arrOut)[ipartOut++]) AliParticleYield(*part);
381 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(*part);
f5b6253e 382 if(TDatabasePDG::Instance()->GetParticle(-particleYields[ipart]) &&
383 (particleYields[ipart] != kPDGLambda) &&
384 (particleYields[ipart] != kPDGKStar)
385 ){// if only the sum was found, add a dummy entry to the
386 // GSI file, so that anton always has the same # of
387 // lines. However, we don't do this for the Lambda and
388 // KStar (always one of the charges or the average)
389 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(-particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
390 }
69f3eeda 391 }
392 else {
f5b6253e 393 std::cout << "ADDING DUMMY sum " << particleYields[ipart] << std::endl;
394 new((*arrOutGSI)[ipartOutGSI++]) AliParticleYield(particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
395 if (particleYields[ipart]==kPDGHyperTriton) {
396 // If is the 3LH, add another one for the antiparticle
397 AliParticleYield(-particleYields[ipart], system, energy, -10, -10, -10, -10, -10, -10, 5, 256, "DUMMY", 1, "ALICE");
398 }
69f3eeda 399 }
9d84731e 400 }
401
402 }
403 }
404 std::cout << "Particles for thermal model fits:" << std::endl;
69f3eeda 405 arrOut->Print("short");
406 // arrOut->Print("");
9d84731e 407 std::cout << "" << std::endl;
408 // Write GSI input file
69f3eeda 409 TIter it(arrOutGSI);
9d84731e 410 AliParticleYield * part = 0;
411 ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
412 while ((part = (AliParticleYield*) it.Next())){
413 fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
414 }
415 fout.close();
416 // Write thermus file
417 AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
418}
ad431d97 419
420
b926af7d 421TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle,
422 Int_t icolor, Int_t imarker, Int_t errorsType, Float_t shift) {
ad431d97 423
b926af7d 424 TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1+shift, npart+1+shift);
ad431d97 425
ad431d97 426 for(Int_t ipart = 1; ipart <= npart; ipart++){
4910214c 427 std::cout << "----- Searching " << particleYields[ipart-1] << " -------" << std::endl;
428
b926af7d 429 AliParticleYield * part = AliParticleYield::FindParticle(arr, particleYields[ipart-1], system, energy, centrality,isSumYields[ipart-1]);
4910214c 430 if(part) {
431 std::cout << "found" << std::endl;
432 part->Print();
433 }
b926af7d 434 if(!part && isSumYields[ipart-1]) {
435 //Could not find the particle, but the sum was requested: build the sum!
436 part = AliParticleYield::FindParticle(arr, particleYields[ipart-1], system, energy, centrality, 0);
437 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particleYields[ipart-1], system, energy, centrality, 0);
4910214c 438 if(part2 && part) {
439 std::cout << " Building sum" << std::endl;
440 part->Print();
441 part2->Print();
442 part = AliParticleYield::Add(part, part2);
443 }
444 else if(part) {
445 std::cout << "Scaling part" << std::endl;
446 part->Print();
447 part = new AliParticleYield(*part); // Always clone before scaling
448 part->Scale(2.); // If we only found a particle, we can scale it by a factor 2.
449 }
450 else part = 0;
b926af7d 451 }
452 if(!part){
453 std::cout << "Cannot find " << particleYields[ipart-1] << std::endl;
454 continue;
455 }
4910214c 456 if(isSumYields[ipart-1] && !part->IsTypeAverage()) {
457 std::cout << " scaling /2" << std::endl;
458 part = new AliParticleYield(*part); // Always clone before scaling
459 part->Scale(0.5); // take average
460 }
b926af7d 461 std::cout << " Plotting " << particleYields[ipart-1] << std::endl;
4910214c 462 part->Print();
b926af7d 463 // part->Scale(scale[ipart-1]);
ad431d97 464 h->SetBinContent(ipart, part->GetYield());
b926af7d 465 if(errorsType == kTotalError) {
466 h->SetBinError (ipart, part->GetTotalError(0/* 0 = no normalization error */));
467 } else if (errorsType == kSystError) {
468 h->SetBinError (ipart, part->GetSystError());
469 } else if (errorsType == kStatError) {
470 h->SetBinError (ipart, part->GetStatError());
471 }
472 h->GetXaxis()->SetBinLabel(ipart, part->GetLatexName());
ad431d97 473
474 }
b926af7d 475 h->SetMarkerStyle(imarker);
476 h->SetMarkerColor(icolor);
4910214c 477 h->SetMarkerSize(1.4);
b926af7d 478
ad431d97 479 return h;
480}
046fa53e 481
482void LoadArrays() {
483 arrPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
484 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
485 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
486 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
487 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
488 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
489 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
490 arrPbPb->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
491
492 arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
493
494 arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt");
495 arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt");
496
497 arrpPb = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt");
498 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
499 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"));
500 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_AveragedNumbers.txt"));
69f3eeda 501 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_phi.txt"));
502 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_Kstar.txt"));
046fa53e 503 arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt");
504
505 // RHIC data
506 arrSTARPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_STAR-AntonQM12.txt");
507 arrPHENIXPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_PHENIX-AntonQM12.txt");
508 arrBRAHMSPbPb = AliParticleYield::ReadFromASCIIFile("PbPb_200_BRAHMS-AntonQM12.txt");
509 arrSTARpp = AliParticleYield::ReadFromASCIIFile("pp_200_STAR.txt");
510 arrPHENIXpp = AliParticleYield::ReadFromASCIIFile("pp_200_PHENIX.txt");
511
512}
513
514void SetStyle(Bool_t graypalette) {
515 std::cout << "Setting style!" << std::endl;
516
517 gStyle->Reset("Plain");
518 gStyle->SetOptTitle(0);
519 gStyle->SetOptStat(0);
520 if(graypalette) gStyle->SetPalette(8,0);
521 else gStyle->SetPalette(1);
69f3eeda 522 gStyle->SetDrawBorder(0);
046fa53e 523 gStyle->SetCanvasColor(10);
524 gStyle->SetCanvasBorderMode(0);
69f3eeda 525 gStyle->SetPadBorderMode(0);
046fa53e 526 gStyle->SetFrameLineWidth(1);
527 gStyle->SetFrameFillColor(kWhite);
528 gStyle->SetPadColor(10);
529 gStyle->SetPadTickX(1);
530 gStyle->SetPadTickY(1);
531 gStyle->SetPadBottomMargin(0.15);
532 gStyle->SetPadLeftMargin(0.15);
533 gStyle->SetHistLineWidth(1);
534 gStyle->SetHistLineColor(kRed);
535 gStyle->SetFuncWidth(2);
536 gStyle->SetFuncColor(kGreen);
537 gStyle->SetLineWidth(2);
538 gStyle->SetLabelSize(0.045,"yz");
539 gStyle->SetLabelSize(0.06,"x");
540 gStyle->SetLabelOffset(0.01,"y");
541 gStyle->SetLabelOffset(0.01,"x");
542 gStyle->SetLabelColor(kBlack,"xyz");
543 gStyle->SetTitleSize(0.05,"xyz");
544 gStyle->SetTitleOffset(1.25,"y");
545 gStyle->SetTitleOffset(1.2,"x");
546 gStyle->SetTitleFillColor(kWhite);
547 gStyle->SetTextSizePixels(26);
548 gStyle->SetTextFont(42);
549 // gStyle->SetTickLength(0.04,"X"); gStyle->SetTickLength(0.04,"Y");
550
551 gStyle->SetLegendBorderSize(0);
552 gStyle->SetLegendFillColor(kWhite);
553 // gStyle->SetFillColor(kWhite);
554 gStyle->SetLegendFont(42);
555
556 gStyle->SetErrorX(0);
557 gStyle->SetEndErrorSize(5);
558}
559
b926af7d 560TLegend * NewLegendQM(Double_t x1, Double_t y1, Double_t x2, Double_t y2, Bool_t isYield) {
046fa53e 561
b926af7d 562 const char * style = "p";
046fa53e 563 const char ** labels=0;
79936afc 564 // Bool_t beautify=kFALSE;
046fa53e 565 Bool_t useTitle=kTRUE;
046fa53e 566
15fccaa0 567 TLegend * l = new TLegend(x1, y1, x2, y2);
568 l->SetFillColor(kWhite);
569 l->SetTextFont(43);
570 l->SetTextSize(25);
571 l->SetBorderSize(1);
572 l->SetLineWidth(1);
b926af7d 573 l->SetMargin(0.1);
046fa53e 574 // const Int_t markers[] = {20,24,21,25,23,28,33,20,24,21,25,23,28,33};
575 // const Int_t colors[] = {1,2,3,4,6,7,8,9,10,11,1,2,3,4,6,7,8,9,10};
576
577 TList * list = gPad->GetListOfPrimitives();
578 TIterator * iter = list->MakeIterator();
579 TObject * obj = 0;
580 Int_t ilabel = -1;
581 while ((obj = (TObject*) iter->Next())){
582 if (obj->InheritsFrom("TH1") || obj->InheritsFrom("TGraph") || obj->InheritsFrom("TF1")) {
583 if( (TString(obj->GetName()) == "hframe" ) ) continue;
584 ilabel++;
585 if (labels != 0)
586 l->AddEntry(obj, labels[ilabel], style);
587 else{
588 if (useTitle) {
589 if(TString(obj->GetTitle()).Contains("NoLegend")) continue;
590 l->AddEntry(obj, obj->GetTitle(), style);
591 }
592 else
593 l->AddEntry(obj, obj->GetName(), style);
594 }
79936afc 595 // if(beautify) {
046fa53e 596
79936afc 597 // if(!obj->InheritsFrom("TF1")){
598 // ((TH1*)obj)->SetLineColor(colors[ilabel]);
599 // ((TH1*)obj)->SetMarkerColor(colors[ilabel]);
600 // ((TH1*)obj)->SetMarkerStyle(markers[ilabel]);
601 // } else {
602 // ((TF1*)obj)->SetLineColor(colors[ilabel]);
603 // }
604
046fa53e 605 }
606 }
b926af7d 607 if(isYield) {
608 // Add some details on excluded stuff
609 l->SetTextSize(16);
b926af7d 610 }
611 l->Draw();
612 return l;
046fa53e 613}
614
615
15fccaa0 616
b926af7d 617void DrawRatio(TString what, Bool_t isYield, Double_t shift) {
046fa53e 618 // This is used to simplify the code above
619 // In order to draw syst error bars, we need to convert to graphs the syst errors histos
b926af7d 620 // if isYield == true plots yields rather than ratios
046fa53e 621
4c7bd9d1 622 // Sample colors
623 // const Int_t colors[] = {kBlack, kRed+1 , kBlue+1, kGreen+3, kMagenta+1, kOrange-1,kCyan+2,kYellow+2 , kWhite};
624
046fa53e 625 TClonesArray * array = 0;
626 Int_t system, color, marker;
b926af7d 627 Float_t energy = 0, shiftloc = shift;
046fa53e 628 TString centrality, label;
4c7bd9d1 629 // FIXME: move this in the different sections below
b926af7d 630 correlatedUnc = correlatedUncZero;
4c7bd9d1 631 std::cout << "Plotting " << what.Data() << std::endl;
632
046fa53e 633
634 if (what == "frame" ) {
4c7bd9d1 635 correlatedUnc = correlatedUncZero;
b926af7d 636 DrawFrame(isYield);
69f3eeda 637 // TH1 * h = GetHistoRatios(arrPbPb, AliParticleYield::kCSPbPb, 2760, "V0M0010", "NoLegend", kWhite);
638 // h->Draw();
639 // h->GetYaxis()->SetDecimals(1);
640 // h->GetYaxis()->SetNdivisions(505);
641 // h->GetXaxis()->CenterLabels(1);
642 // Int_t nratio = h->GetNbinsX();
b926af7d 643 if(!isYield) {
644 for(Int_t iratio = 0; iratio < nratio; iratio++){
645 Double_t exp = TMath::Floor(TMath::Log10(TMath::Abs(scale[iratio])));
646 Double_t man = scale[iratio] / TMath::Power(10, exp);
647 if(exp > 2) {
648 // TLatex * scaleLabel = new TLatex(iratio+1+0.2,maxy*1.01, Form("#times %0.0f 10^{%0.0f}", man, exp));
649 TLatex * scaleLabel = new TLatex(iratio+1+0.2,0.005, Form("#times %0.0f 10^{%0.0f}", man, exp));
650 scaleLabel->SetTextFont(43);
651 scaleLabel->SetTextSize(20);
652 scaleLabel->Draw();
653 } else {
654 Double_t shiftloc2 = scale[iratio] < 50 ? 0.3 : 0.2;
655 TLatex * scaleLabel = new TLatex(iratio+1+shiftloc2, 0.005, Form("#times %g", scale[iratio]));
656 scaleLabel->SetTextFont(43);
657 scaleLabel->SetTextSize(20);
658 scaleLabel->Draw();
659 }
660
661 }
046fa53e 662 }
4910214c 663 TPaveText *pt = new TPaveText( 0.176, 0.842881, 0.378514, 0.929595,"brNDC");
15fccaa0 664 pt->SetBorderSize(0);
665 pt->SetFillColor(0);
666 pt->SetLineWidth(1);
667 pt->SetTextFont(43);
668 pt->SetTextSize(23);
669 pt->AddText("ALICE Preliminary");
670 pt->Draw();
671 // pt = new TPaveText( 0.167671, 0.815206, 0.375502, 0.865021,"brNDC");
672 // pt->SetBorderSize(0);
673 // pt->SetFillColor(0);
674 // pt->SetLineWidth(1);
675 // pt->SetTextFont(43);
676 // pt->SetTextSize(20);
677 // pt->AddText("particle+antiparticle");
678 // pt->Draw();
046fa53e 679
4c7bd9d1 680 gPad->SetGridx();
046fa53e 681
682 }
683 else if (what == "PbPb_0010") {
684 array = arrPbPb;
685 system = 2; energy = 2760.; centrality = "V0M0010";
4910214c 686 label = "Pb-Pb #sqrt{s}_{NN} = 2.76 TeV, 0-10%";
79936afc 687 color = kRed+1;
046fa53e 688 marker = kFullCircle;
b926af7d 689 if(!shift) shiftloc = 0;
4c7bd9d1 690 correlatedUnc = correlatedUncLocalPbPb;
691
046fa53e 692 }
693
694 else if (what == "PbPb_6080") {
79936afc 695 array = arrPbPb;
696 system = 2; energy = 2760.; centrality = "V0M6080";
4910214c 697 label = "Pb-Pb #sqrt{s}_{NN} = 2.76 TeV, 60-80%";
79936afc 698 color = kBlue+1;
699 marker = kFullCircle;
b926af7d 700 if(!shift) shiftloc = 0.0;
4c7bd9d1 701 correlatedUnc = correlatedUncLocalPbPb;
046fa53e 702 }
703 else if (what == "pp7") {
79936afc 704 array = arrpp7;
705 system = 0; energy = 7000.; centrality = "";
706 label = "pp #sqrt{s} = 7 TeV";
4c7bd9d1 707 color = kMagenta+1;
79936afc 708 marker = kFullCircle;
b926af7d 709 if(!shift) shiftloc = 0.2;
4c7bd9d1 710 correlatedUnc = correlatedUncLocalPP;
046fa53e 711 }
712 else if (what == "pp900") {
79936afc 713 array = arrpp900;
714 system = 0; energy = 900.; centrality = "";
715 label = "pp #sqrt{s} = 0.9 TeV";
4c7bd9d1 716 color = kCyan+2;
79936afc 717 marker = kFullCircle;
b926af7d 718 if(!shift) shiftloc = -0.2;
4c7bd9d1 719 correlatedUnc = correlatedUncLocalPP;
046fa53e 720 }
721 else if (what == "pp276") {
79936afc 722 array = arrpp276;
723 system = 0; energy = 2760.; centrality = "";
724 label = "pp #sqrt{s} = 2.76 TeV";
4c7bd9d1 725 color = kYellow+2;
79936afc 726 marker = kFullCircle;
b926af7d 727 if(!shift) shiftloc = 0;
4c7bd9d1 728 correlatedUnc = correlatedUncLocalPP;
046fa53e 729 }
730 else if (what == "pPb0005") {
79936afc 731 array = arrpPb;
4c7bd9d1 732 system = 1; energy = 5020.; centrality = "V0A0005";
4910214c 733 label = "p-Pb #sqrt{s}_{NN} = 5.02 TeV, V0A 0-5%";
79936afc 734 color = kBlack;
735 marker = kFullCircle;
b926af7d 736 if(!shift) shiftloc = -0.2;
4c7bd9d1 737 correlatedUnc = correlatedUncLocalPP;
79936afc 738 }
739 else if (what == "PbPbSTAR") {
046fa53e 740 array = arrSTARPbPb;
4c7bd9d1 741 system = 2; energy = 200.; centrality = "0005";
4910214c 742 label = "STAR, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-5%";
046fa53e 743 color = kBlack;
744 marker = kOpenStar;
b926af7d 745 if(!shift) shiftloc = +0.2;
4c7bd9d1 746 correlatedUnc = correlatedUncZero;
79936afc 747 }
748 else if (what == "PbPbPHENIX") {
749 array = arrPHENIXPbPb;
4c7bd9d1 750 system = 2; energy = 200.; centrality = "0005";
4910214c 751 label = "PHENIX, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-5%";
79936afc 752 color = kBlack;
753 marker = kOpenSquare;
b926af7d 754 if(!shift) shiftloc = -0.15;
4c7bd9d1 755 correlatedUnc = correlatedUncZero;
79936afc 756 }
757 else if (what == "PbPbBRAHMS") {
758 array = arrBRAHMSPbPb;
4c7bd9d1 759 system = 2; energy = 200.; centrality = "0010";
4910214c 760 label = "BRAHMS, Au-Au #sqrt{s}_{NN} = 0.2 TeV, 0-10%";
79936afc 761 color = kBlack;
762 marker = kOpenCross;
b926af7d 763 if(!shift) shiftloc = -0.3;
4c7bd9d1 764 correlatedUnc = correlatedUncZero;
765 }
766 else if (what == "ppSTAR") {
767 array = arrSTARpp;
768 system = 0; energy = 200.; centrality = "";
4910214c 769 label = "STAR, pp #sqrt{s} = 0.2 TeV";
4c7bd9d1 770 color = kBlack;
771 marker = kOpenStar;
b926af7d 772 if(!shift) shiftloc = 0.;
4c7bd9d1 773 correlatedUnc = correlatedUncZero;
774 }
775 else if (what == "ppPHENIX") {
776 array = arrPHENIXpp;
777 system = 0; energy = 200.; centrality = "";
4910214c 778 label = "PHENIX, pp #sqrt{s} = 0.2 TeV";
4c7bd9d1 779 color = kBlack;
780 marker = kOpenSquare;
b926af7d 781 if(!shift) shiftloc = -0.2;
4c7bd9d1 782 correlatedUnc = correlatedUncZero;
69f3eeda 783 }
784 // From here on, it's meta names, to draw several series of ratios
785 else if (what == "allpp"){
4910214c 786 DrawRatio("frame");
15fccaa0 787 DrawRatio("pp900");
4910214c 788 DrawRatio("pp276");
789 DrawRatio("pp7");
15fccaa0 790 array =0;
791 NewLegendQM(0.62249, 0.635734, 0.910643, 0.94673);
792
793 myCan->Update();
794 gSystem->ProcessEvents();
795 myCan->Print("Ratios_pponly.eps");
796 gSystem->Exec("epstopdf Ratios_pponly.eps");
797 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_pponly.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
798 }
799 else if (what == "PbPbWithPP7TeV"){
800
801 DrawRatio("frame");
802 DrawRatio("PbPb_0010");
803 DrawRatio("pp7");
804 array =0;
805
806 NewLegendQM(0.413655, 0.748094, 0.910643, 0.948736);
807
808 myCan->Update();
809 gSystem->ProcessEvents();
810 myCan->Print("Ratios_withpp7tev.eps");
811 gSystem->Exec("epstopdf Ratios_withpp7tev.eps");
812 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_withpp7tev.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
813 } else if(what == "allsyst") {
814
815 DrawRatio("frame");
b926af7d 816 DrawRatio("pp7", 0, -0.2);
817 DrawRatio("pPb0005", 0, 0.00001);
4910214c 818 DrawRatio("PbPb_0010", 0, 0.2);
69f3eeda 819 array =0;
15fccaa0 820
15fccaa0 821
b926af7d 822 NewLegendQM(0.462851, 0.631722, 0.89257, 0.936697);
823 DrawExtrapolatedSymbolsAndLegendPbPb0010();
824 DrawExtrapolatedSymbolsAndLegendpPb0005();
15fccaa0 825 myCan->Update();
826 gSystem->ProcessEvents();
827 myCan->Print("Ratios_allsystems.eps");
828 gSystem->Exec("epstopdf Ratios_allsystems.eps");
829 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_allsystems.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
830
831 }else if(what =="PbPb6080andpPb0005") {
832 DrawRatio("frame");
15fccaa0 833 DrawRatio("pPb0005");
4910214c 834 DrawRatio("PbPb_6080");
15fccaa0 835 array=0;
836
4910214c 837 NewLegendQM( 0.46988, 0.730036, 0.910643, 0.948736);
b926af7d 838 DrawExtrapolatedSymbolsAndLegendpPb0005();
15fccaa0 839 myCan->Update();
840 gSystem->ProcessEvents();
841 myCan->Print("Ratios_6080vspPb.eps");
842 gSystem->Exec("epstopdf Ratios_6080vspPb.eps");
843 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_6080vspPb.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
844
845
846 }else if(what =="pp_vsRHIC") {
847 DrawRatio("frame");
848 DrawRatio("pp7");
849 DrawRatio("ppSTAR");
850 DrawRatio("ppPHENIX");
851 array=0;
852
b926af7d 853 NewLegendQM( 0.554217, 0.677869, 0.910643, 0.948736);
15fccaa0 854
855 myCan->Update();
856 gSystem->ProcessEvents();
857 myCan->Print("Ratios_vsRHIC_pp.eps");
858 gSystem->Exec("epstopdf Ratios_vsRHIC_pp.eps");
859 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_vsRHIC_pp.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
860
861 } else if (what =="PbPb_vsRHIC") {
862 DrawRatio("frame");
863 DrawRatio("PbPb_0010");
864 DrawRatio("PbPbSTAR");
865 DrawRatio("PbPbPHENIX");
866 DrawRatio("PbPbBRAHMS");
867 array = 0;
868
869
b926af7d 870
871 NewLegendQM( 0.434739, 0.591593, 0.939759, 0.936697);
872 DrawExtrapolatedSymbolsAndLegendPbPb0010();
15fccaa0 873
874 myCan->Update();
875 gSystem->ProcessEvents();
876 myCan->Print("Ratios_vsRHIC_PbPb.eps");
877 gSystem->Exec("epstopdf Ratios_vsRHIC_PbPb.eps");
878 gSystem->Exec("if [ \"$USER\" = \"mfloris\" ]; then cp Ratios_vsRHIC_PbPb.{eps,pdf} /Users/mfloris/Documents/PapersNTalks/ALICE/ThermalFits/img/; fi ");
879 } else if (what == "aliceall") {
880 DrawRatio("frame");
881 DrawRatio("PbPb_0010");
882 DrawRatio("PbPb_6080");
883 DrawRatio("pPb0005");
884 DrawRatio("pp7");
885 DrawRatio("pp276");
886 DrawRatio("pp900");
046fa53e 887 }
888
b926af7d 889 // FROM HERE: IT's yields
890 else if( what == "fit_ReferenceFit_PbPb0010") {
891
892 DrawRatio("frame",1);
893 DrawRatio("PbPb_0010",1);
4910214c 894 legThermal->SetNColumns(4);
895 AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
b926af7d 896 // FIXME: sistemare valori rapporti
4910214c 897 PlotThermusYields("lhc2760_final_0005_single_gc_output_gs1_wevc_nkst.txt" , kBlack , kSolid , "Thermus 2.3 , 155 #pm 2 , 5882 #pm 532 , 24.6/9");
898 PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat" , kRed+1 , kDashed , "GSI , 156 #pm 1.5 , 5330 #pm 505 , 17.4/9" );
899 PlotThermusYields("fit_gamma_q_s_fixed_with_nuclei_PbPb_0010.txt" , kBlue+1 , kDashDotted , "SHARE 3 , 156 #pm 2 , 4479 #pm 639 , 12.5/8");
b926af7d 900
4910214c 901 TLegend * l = NewLegendQM(0.682731, 0.701578, 0.939759, 0.920947, 1);
902 TMarker * m = new TMarker(0.1, 0.1,markerNoFit);
903 m->SetMarkerSize(1.2);
904 l->AddEntry( m, "Not in fit", "p");
905 m = new TMarker(0.1,0.1,markerExtrap);
906 m->SetMarkerSize(1.2);
907 l->AddEntry( m, "Extrapolated (Pb-Pb 0-10%)", "p");
908 DrawMarkerKStarNoFit() ;
909 DrawExtrapolatedSymbolsYieldsPbPb0010();
b926af7d 910
911
912 array =0;
4910214c 913 } else if( what == "fitSHARE_NoPionsNoProtons_PbPb0010") {
b926af7d 914 array =0;
4910214c 915
916 DrawRatio("frame",1);
917 DrawRatio("PbPb_0010",1);
918 legThermal->SetNColumns(4);
919 AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
920
921 PlotThermusYields("fit_gamma_q_s_fixed_PbPb_0010.txt" ,kBlack , kSolid ,"SHARE 3 , 155 #pm 4 , 4583 #pm 853 , 10.1/5 ");
922 PlotThermusYields("fit_gamma_q_s_fixed_wo_pions_PbPb_0010.txt" ,kRed+1 , kDashed ,"SHARE 3 (no #pi) , 158 #pm 4 , 3801 #pm 738 , 7.2/4");
923 PlotThermusYields("fit_gamma_q_s_fixed_wo_protons_PbPb_0010.txt" ,kBlue+1 , kDashDotted ,"SHARE 3 (no p) , 159 #pm 3 ,3898 #pm 676 , 2.7/4");
924
925
926 TLegend * l = NewLegendQM(0.682731, 0.701578, 0.939759, 0.920947, 1);
927 TMarker * m = new TMarker(0.1, 0.1,markerNoFit);
928 m->SetMarkerSize(1.2);
929 l->AddEntry( m, "Not in fit", "p");
930 m = new TMarker(0.1,0.1,markerExtrap);
931 m->SetMarkerSize(1.2);
932 l->AddEntry( m, "Extrapolated (Pb-Pb 0-10%)", "p");
933
934 // Add markers for additional particles not in fit
935 DrawMarkerKStarNoFit() ;
936 DrawExtrapolatedSymbolsYieldsPbPb0010();
b926af7d 937
4910214c 938 myPadLabel->cd();
939 Double_t ymarker = 0.111825;
940 TMarker * marker = new TMarker(0.4267068,ymarker,28);
941 marker->SetMarkerStyle(28);
942 marker->SetMarkerSize(1.2);
943 marker->Draw();
944 marker = new TMarker(0.7881526,ymarker,28);
945 marker->SetMarkerStyle(28);
946 marker->SetMarkerSize(1.2);
947 marker->Draw();
948 marker = new TMarker(0.8644578,ymarker,28);
949 marker->SetMarkerStyle(28);
950 marker->SetMarkerSize(1.2);
951 marker->Draw();
952 marker = new TMarker(0.9558233,ymarker,28);
953 marker->SetMarkerStyle(28);
954 marker->SetMarkerSize(1.2);
955 marker->Draw();
956 myPadHisto->cd();
957
958
959 } else if( what == "fitGSI_NoPionsNoProtons_PbPb0010") {
b926af7d 960 array =0;
4910214c 961
962 DrawRatio("frame",1);
963 DrawRatio("PbPb_0010",1);
964 legThermal->SetNColumns(4);
965 AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
966
967 PlotGSIYields("data+therm_fit2_s2760_0-10qm14.dat" , kBlack , kSolid , "GSI , 156 #pm 1.5 , 5330 #pm 505 , 17.4/9" );
968 PlotGSIYields("data+therm_fit2_s2760_0-10qm14NOp.dat" , kRed+1 , kDashed , "GSI (no p) , 156 #pm 2 , 5590 #pm 330, 7.5/8" );
969 PlotGSIYields("data+therm_fit2_s2760_0-10qm14NOpi.dat", kBlue+1, kDashDotted , "GSI (no #pi), 157 #pm 2, 4990 #pm 630, 15.9/8");
970
971
972 TLegend * l = NewLegendQM(0.682731, 0.701578, 0.939759, 0.920947, 1);
973 TMarker * m = new TMarker(0.1, 0.1,markerNoFit);
974 m->SetMarkerSize(1.2);
975 l->AddEntry( m, "Not in fit", "p");
976 m = new TMarker(0.1,0.1,markerExtrap);
977 m->SetMarkerSize(1.2);
978 l->AddEntry( m, "Extrapolated (Pb-Pb 0-10%)", "p");
979
980 // Add markers for additional particles not in fit
981 DrawMarkerKStarNoFit() ;
982 DrawExtrapolatedSymbolsYieldsPbPb0010();
983
984
985 } else if (what == "fitShareWithWithoutNuclei_PbPb0010") {
986 DrawRatio("frame",1);
987 DrawRatio("PbPb_0010",1);
988 legThermal->SetNColumns(4);
989 AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
990
991 PlotThermusYields("fit_gamma_q_s_fixed_PbPb_0010.txt" ,kBlack , kSolid ,"SHARE 3 (no nuclei) , 155 #pm 4 , 4583 #pm 853 , 10.1/5 ");
992 PlotThermusYields("fit_gamma_q_s_fixed_with_nuclei_PbPb_0010.txt" , kBlue+1 , kDashDotted , "SHARE 3 , 156 #pm 2 , 4479 #pm 639 , 12.5/8");
993
994 TLegend * l = NewLegendQM(0.682731, 0.701578, 0.939759, 0.920947, 1);
995 TMarker * m = new TMarker(0.1, 0.1,markerNoFit);
996 m->SetMarkerSize(1.2);
997 l->AddEntry( m, "Not in fit", "p");
998 m = new TMarker(0.1,0.1,markerExtrap);
999 m->SetMarkerSize(1.2);
1000 l->AddEntry( m, "Extrapolated (Pb-Pb 0-10%)", "p");
1001
1002 DrawMarkerKStarNoFit() ;
1003 DrawExtrapolatedSymbolsYieldsPbPb0010();
b926af7d 1004
1005 } else if( what == "fitThermus_GammaSFree_PbPb0010") {
1006 std::cout << "MISSING DATA" << std::endl;
1007 array =0;
1008 } else if( what == "fitShare_GammaSGammaQFree_PbPb0010") {
1009 std::cout << "MISSING DATA" << std::endl;
1010 array =0;
1011 } else if( what == "fitShare_All_PbPb0010") {
1012 DrawRatio("frame",1);
1013 DrawRatio("PbPb_0010",1);
4910214c 1014 // PlotThermusYields("fit_gamma_q_s_fixed_PbPb_0010.txt" , kBlack , kSolid , "SHARE 3" , 100 , 100 , 100 , 100 , -1 , -1 , -2 , -1 , -1 , -1);
1015 // PlotThermusYields("fit_gamma_q_fixed_PbPb_0010.txt" , kBlue , kSolid , "SHARE 3" , 100 , 100 , 100 , 100 , -1 , -1 , -2 , -1 , -1 , -1);
1016 // PlotThermusYields("fit_gamma_q_s_free_PbPb_0010.txt" , kRed , kSolid , "SHARE 3" , 100 , 100 , 100 , 100 , -1 , -1 , -2 , -1 , -1 , -1);
b926af7d 1017 array =0;
1018 } else if( what == "fitGSI_PbPb6080") {
4910214c 1019 DrawRatio("frame",1);
1020 DrawRatio("PbPb_6080",1);
1021 legThermal->SetNColumns(4);
1022 AddLineToThermalLegend(legThermal, "Model,T (MeV), V (fm^{3}), #chi^{2}/NDF", "0");
1023 PlotGSIYields("data+therm_fit2_s2760_60-80qm14.dat", kBlack, kSolid, "GSI, 157 #pm 2, 210 #pm 20, 9.1/7");
1024
1025 legThermal->SetX1NDC(0.143574);
1026 legThermal->SetY1NDC(0.0702247);
1027 legThermal->SetX2NDC(0.659639);
1028 legThermal->SetY2NDC(0.238764);
1029 legThermal->Draw();
1030
1031 TLegend * l = NewLegendQM( 0.682731, 0.744382, 0.939759, 0.920947, 1);
1032 TMarker * m = new TMarker(0.1, 0.1,markerNoFit);
1033 m->SetMarkerSize(1.2);
1034 l->AddEntry( m, "Not in fit", "p");
1035 DrawMarkerKStarNoFit() ;
1036 //DrawExtrapolatedSymbolsYieldsPbPb0010();
1037
b926af7d 1038 array =0;
1039 } else if( what == "fitGSI_pPb0005") {
1040 std::cout << "MISSING DATA" << std::endl;
1041 array =0;
1042 } else if( what == "fitGSI_pPb2040") {
1043 std::cout << "MISSING DATA" << std::endl;
1044 array =0;
1045 } else if( what == "fitThermus_GammaSFree_pPb0005") {
1046 std::cout << "MISSING DATA" << std::endl;
1047 array =0;
1048 } else if( what == "fitThermus_GammaSFree_pPb2040") {
1049 std::cout << "MISSING DATA" << std::endl;
1050 array =0;
1051 } else if( what == "fitThermus_GammaSFree_pPb6080") {
1052 std::cout << "MISSING DATA" << std::endl;
1053 array =0;
1054 } else if( what == "fitThermus_RC_pPb0005") {
1055 std::cout << "MISSING DATA" << std::endl;
1056 array =0;
1057 } else if( what == "fitThermus_RC_pPb2040") {
1058 std::cout << "MISSING DATA" << std::endl;
1059 array =0;
1060 } else if( what == "fitThermus_RC_pPb6080") {
1061 std::cout << "MISSING DATA" << std::endl;
1062 array =0;
1063 } else if( what == "fitShare_pPb0005") {
1064 DrawRatio("frame",1);
1065 DrawRatio("pPb0005",1);
4910214c 1066 // PlotThermusYields("fit_gamma_q_fixed_pPb_0005.txt" , kBlue , kDashDotted , "SHARE 3" , 161 , 3 , 114 , 17 , 0.93 ,0.04 ,-2 ,0 ,26.0 ,6);
1067 // PlotThermusYields("fit_gamma_q_s_free_pPb_0005.txt" , kRed , kDashDotted , "SHARE 3" , 162 , 3 , 111 , 16 , 0.93 ,0.03 ,-2 ,0 ,61.4 ,6);
1068 // PlotThermusYields("fit_gamma_q_s_fixed_pPb_0005.txt" , kGreen , kDashDotted , "SHARE 3" , 158 , 2 , 125 , 16 , -1 ,-1 ,-2 ,0 ,30.8 ,7);
1069 // PlotThermusYields("fit_gamma_q_s_free_with_d_pPb_0005.txt", kBlack, kSolid, "SHARE 3", 163, 5, 124, 15, 0.85, 0.08, 0.92, 0.06, 29.3, 7);
b926af7d 1070 NewLegendQM(0.613454, 0.701578, 0.940763, 0.918272, 1);
1071 array =0;
1072
1073 } else if( what == "fitShare_pPb2040") {
1074 std::cout << "MISSING DATA" << std::endl;
1075 array =0;
1076 } else if( what == "fitShare_pPb6080") {
1077 std::cout << "MISSING DATA" << std::endl;
1078 array =0;
1079 } else if( what == "fitGSI_pp") {
1080 std::cout << "MISSING DATA" << std::endl;
1081 array =0;
1082 } else if( what == "fitGSI_fullCanonical") {
1083 std::cout << "MISSING DATA" << std::endl;
1084 array =0;
1085 } else if( what == "fitThermus_RC_pp") {
1086 std::cout << "MISSING DATA" << std::endl;
1087 array =0;
1088 } else if( what == "fitShare_pp") {
1089 std::cout << "MISSING DATA" << std::endl;
1090 array =0;
1091 }
1092
046fa53e 1093
15fccaa0 1094
046fa53e 1095 else {
1096 std::cout << "Unknown ratio " << what.Data() << std::endl;
15fccaa0 1097 return;
046fa53e 1098 }
1099
4c7bd9d1 1100 if(!correlatedUnc) {
1101 std::cout << "correlatedUnc not set!" << std::endl;
1102
1103 }
1104 std::cout << "CORR: " << correlatedUnc[1] << std::endl;
1105
046fa53e 1106 if(array) {
b926af7d 1107 if(isYield) {
4910214c 1108 TGraphErrors * hstat =
b926af7d 1109 AliPWGHistoTools::GetGraphFromHisto(GetHistoYields(array, system, energy, centrality, label, color, marker, kStatError, shiftloc)
4910214c 1110 ,0);
1111 hstat->SetMarkerSize(1.2);
1112 hstat->Draw("PZ");
b926af7d 1113 AliPWGHistoTools::GetGraphFromHisto(GetHistoYields(array, system, energy, centrality, label+"NoLegend", color, marker, kSystError, shiftloc)
1114 ,0)->Draw("[]");
1115 } else {
1116 AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array, system, energy, centrality, label, color, marker, kStatError, shiftloc)
1117 ,0)->Draw("PZ");
1118 AliPWGHistoTools::GetGraphFromHisto(GetHistoRatios(array, system, energy, centrality, label+"NoLegend", color, marker, kSystError, shiftloc)
1119 ,0)->Draw("[]");
1120 }
046fa53e 1121 }
1122
1123
1124}
b926af7d 1125
1126void DrawFrame(Bool_t isYield) {
69f3eeda 1127
1128 myCan = new TCanvas("myCan","BD april.2014",50,10,1000,650);
1129 myCan->Draw();
1130 myCan->cd();
1131 // Set the Pads
15fccaa0 1132 Double_t boundaryLabels = 0.85;//0.92;
b926af7d 1133 myPadLabel = new TPad("myPadLabel","myPadLabel",0.0, boundaryLabels,1.0,1.0);
69f3eeda 1134 myPadSetUp(myPadLabel);
1135 myPadLabel->Draw();
1136
b926af7d 1137 myPadHisto = new TPad("myPadHisto","myPadHisto",0.0,isYield ? 0.25 : 0.05 ,1.0, boundaryLabels,0);
69f3eeda 1138 myPadSetUp(myPadHisto);
1139 myPadHisto->Draw();
1140
b926af7d 1141 myPadStdDev = new TPad("myPadStdDev","myPadStdDev",0.0,0.0,1.0,0.25,0);
69f3eeda 1142 myPadSetUp(myPadStdDev);
b926af7d 1143 if(isYield) myPadStdDev->Draw();
1144 myPadStdDev->SetGridx();
1145 myPadStdDev->SetGridy();
69f3eeda 1146 myPadLabel->cd();
1147
15fccaa0 1148 double xLabelPosition[nratio] = {0.124498, 0.211847, 0.31, 0.38, 0.465, 0.575, 0.644, 0.72, 0.82, 0.905 };
b926af7d 1149 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 1150 double yLabelPosition = 0.40;
69f3eeda 1151
1152 // labels
b926af7d 1153 if(isYield) {
1154 for(Int_t ipart = 0; ipart < npart; ipart++){
1155 TLatex *myPart = new TLatex(xLabelYields[ipart],yLabelPosition, yieldsLabels[ipart]);
1156 myLatexDraw(myPart,20);
1157 }
1158 } else {
1159 for(Int_t iratio = 0; iratio < nratio; iratio++){
1160 TLatex *myRatio = new TLatex(xLabelPosition[iratio],yLabelPosition, ratiosLabels[iratio]);
1161 myLatexDraw(myRatio,20);
1162 }
69f3eeda 1163 }
15fccaa0 1164 // Xi's and Omega's bar (there was no way to convince root to draw it properly)
b926af7d 1165 if(isYield) {
1166 TLine *line = new TLine(0.653,0.75,0.663,0.75);
1167 line->SetLineWidth(2);
1168 line->Draw();
1169 line = new TLine(0.728,0.75,0.741,0.75);
1170 line->SetLineWidth(2);
1171 line->Draw();
1172
b926af7d 1173
1174 }
1175 else {
1176 TLine *line = new TLine(0.408,0.75,0.418,0.75);
1177 line->SetLineWidth(2);
1178 line->Draw();
1179 line = new TLine(0.498,0.75,0.513,0.75);
1180 line->SetLineWidth(2);
1181 line->Draw();
1182
1183
1184 }
15fccaa0 1185
69f3eeda 1186
b926af7d 1187 if(isYield) {
69f3eeda 1188 myPadStdDev->cd();
1189 myPadStdDev->SetTopMargin(0.0);
1190 myPadStdDev->SetTicks(1,1);
1191
1192 Float_t devMax = 3.5;
1193
b926af7d 1194 TH2F *myBlankStdDev = new TH2F("myBlankStdDev","myBlankStdDev",npart,1,npart+1,10,-devMax,+devMax);
69f3eeda 1195 myBlankStdDev->GetXaxis()->SetLabelFont(43); // precision 3: size will be in pixels
1196 myBlankStdDev->GetYaxis()->SetLabelFont(43);
1197 myBlankStdDev->GetYaxis()->SetTitleFont(43);
1198 myBlankStdDev->SetLabelSize(23,"xy");
b926af7d 1199 myBlankStdDev->SetTitleSize(20,"y");
1200 myBlankStdDev->SetNdivisions(20,"x");
69f3eeda 1201 myBlankStdDev->SetNdivisions(605,"y");
1202 myBlankStdDev->SetLabelOffset(0.012,"xy");
b926af7d 1203 myBlankStdDev->SetYTitle("(model-data)/#sigma_{data}");
1204 myBlankStdDev->SetTitleOffset(1.3,"y");
69f3eeda 1205 myBlankStdDev->Draw();
1206 }
1207
1208 myPadHisto->cd();
1209 myPadHisto->SetBottomMargin(0.01);
1210 // myPadHisto->SetLogy();
1211 myPadHisto->SetTicks(1,1);
b926af7d 1212 TH1 *myBlankHisto = 0;
1213 if(isYield) {
1214 myBlankHisto = new TH2F("NoLegend","NoLegend",npart,1,npart+1,10, 0.00002, maxy+0.01 );
1215 gPad->SetLogy();
1216
1217 } else {
1218 myBlankHisto = new TH2F("NoLegend","NoLegend",nratio,1,nratio+1,10, 0, maxy+0.01 );
1219 }
69f3eeda 1220 myBlankHisto->GetXaxis()->SetLabelFont(43); // precision 3: size will be in pixels
1221 myBlankHisto->GetYaxis()->SetLabelFont(43);
1222 myBlankHisto->GetYaxis()->SetTitleFont(43);
1223 myBlankHisto->SetLabelSize(23,"xy");
1224 myBlankHisto->SetTitleSize(26,"y");
1225 myBlankHisto->SetMaximum(10);
1226 myBlankHisto->SetMinimum(0);
b926af7d 1227 myBlankHisto->SetNdivisions(isYield? 20 :10,"x");
69f3eeda 1228 myBlankHisto->SetNdivisions(505,"y");
b926af7d 1229 if(isYield) myBlankHisto->SetYTitle("d#it{N}/d#it{y}");
69f3eeda 1230 myBlankHisto->SetLabelOffset(0.012,"xy");
1231 myBlankHisto->SetTitleOffset(1,"y");
1232 myBlankHisto->Draw();
b926af7d 1233
1234 if(isYield) {
1235 legThermal= new TLegend(0.144578, 0.0702247, 0.659639, 0.383226);
1236 legThermal->SetBorderSize(1);
1237 legThermal->SetTextFont(43);
1238 legThermal->SetTextSize(14);
4910214c 1239 // legThermal->SetNColumns(6);
b926af7d 1240 legThermal->SetFillColor(0);
1241 legThermal->SetLineWidth(1);
1242 legThermal->Draw();
b926af7d 1243 }
69f3eeda 1244}
1245
1246void myLatexDraw(TLatex *currentLatex, Float_t currentSize, Int_t currentColor){
15fccaa0 1247 currentLatex->SetTextFont(43);
1248 currentLatex->SetTextSizePixels(Int_t(currentSize));
1249 // currentLatex->SetTextAngle(0);
69f3eeda 1250 currentLatex->Draw();
1251 return;
1252}
1253
1254void myPaveSetup(float rRatio, float rRange3, float rRange5,
1255 int rFillColor){
1256
1257 float cHiRange = 0, cLoRange = 0;
1258 if(rRange3<rRange5) {cHiRange=rRange5;cLoRange=rRange3;}
1259 else {cHiRange=rRange3;cLoRange=rRange5;}
1260
1261 TPave *cPave= new TPave(rRatio-0.25,cLoRange,rRatio+0.25,cHiRange,0,"br");
1262 cPave->SetFillColor(rFillColor);
1263 cPave->SetLineColor(1);
1264 cPave->Draw();
1265}
1266
1267void myPadSetUp(TPad *currentPad){
1268 currentPad->SetLeftMargin(0.10);
1269 currentPad->SetRightMargin(0.02);
1270 currentPad->SetTopMargin(0.02);
1271 currentPad->SetBottomMargin(0.02);
1272 return;
1273}
b926af7d 1274TGraphErrors* PlotThermusYields(const char * filename, Int_t color, Int_t lineStyle,
4910214c 1275 const char * tag) {
b926af7d 1276
1277 Int_t lw = lineStyle == kSolid ? 2 : 3; // Set line width
1278
1279
1280 std::map<Int_t,Double_t> mapYields;
1281 std::map<Int_t,Double_t> mapStdDev;
1282
1283 Int_t pdg;
1284 Double_t yield, stddev;
1285 ifstream thermusFile(filename);
1286 TString line;
1287 std::cout << "---"<<tag<<"---" << std::endl;
1288
1289 // std::istream is(thermusFile);
1290 // Read the std dev and the ratio in 2 maps, then plot them in a graph.
1291 while(line.ReadLine(thermusFile, kTRUE)) {
1292 if(line.BeginsWith("#")) continue;
1293 TObjArray * tokens = line.Tokenize(" \t");
1294 if(tokens->GetEntries() != 3) continue;// not a line with data
1295 // thermusFile >> pdg >> yield >> stddev;
1296 pdg = ((TObjString*)tokens->At(0))->String().Atof();
1297 yield = ((TObjString*)tokens->At(1))->String().Atof();
1298 stddev = ((TObjString*)tokens->At(2))->String().Atof();
1299
1300 if( thermusFile.eof() ) break;
1301 std::cout << "PDG " << pdg << " " << yield << " " << stddev << std::endl;
1302 mapYields[TMath::Abs(pdg)] += yield;
1303 mapStdDev[TMath::Abs(pdg)] += stddev;
1304 if(pdg < 0) {
1305 // found the antiparticle: now compute the mean
1306 mapYields[TMath::Abs(pdg)] /=2;
1307 mapStdDev[TMath::Abs(pdg)] /=2;
1308 }
1309 delete tokens;
1310 }
1311 // Now plot
1312 TGraphErrors * gThermus = new TGraphErrors;
1313 TGraphErrors * gThermusStdDev = new TGraphErrors;
1314 for(Int_t ipart = 0; ipart < npart; ipart++){
1315 gThermus->SetPoint(ipart, ipart+1.5, mapYields[particleYields[ipart]]);
1316 gThermus->SetPointError(ipart, 0.3, 0);
1317 gThermusStdDev->SetPoint(ipart, ipart+1.5, mapStdDev[particleYields[ipart]]);
1318 gThermusStdDev->SetPointError(ipart, 0.3, 0);
1319 }
1320
1321 myPadHisto->cd();
1322 gThermus->Draw("PZ");
1323 gThermus->SetLineWidth(lw);
1324
1325 gThermus->SetLineColor(color);
1326 gThermus->SetLineStyle(lineStyle);
1327 gThermus->SetTitle("NoLegend");
1328
1329
1330 myPadStdDev->cd();
1331 gThermusStdDev->Draw("PZ");
1332 gThermusStdDev->SetLineWidth(lw);
1333 gThermusStdDev->SetLineColor(color);
1334 gThermusStdDev->SetLineStyle(lineStyle);
1335 myPadHisto->cd();
1336
4910214c 1337
1338 AddLineToThermalLegend(gThermus, tag, "l");
b926af7d 1339 return gThermus;
1340}
1341
1342
1343
1344TGraphErrors* PlotGSIYields(const char * filename, Int_t color, Int_t lineStyle,
4910214c 1345 const char * tag) {
1346
1347 // tag is a comma separated list of elements to be added to the legend as diferent columns
b926af7d 1348
1349 Int_t lw = lineStyle == kSolid ? 2 : 3; // Set line width
1350
1351 const Int_t pdgPbPb0010[] = {211, -211, 321, -321, 310, 313, 333, 2212, -2212, 3122, 3312, -3312, 3334, -3334, 1000010020, 1000020030, 1010010030, -1010010030};
1352 const Int_t pdgPbPb6080[] = {211 , -211 , 321 , -321 , 310 , 313 , 333 , 2212 , -2212 , 3122 , 3312 , -3312 , 3334 , -3334 ,1000010020};
1353 const Int_t pdgpPb0005[] = {211, 321, 310, 313, 333, 2212, 3122, 3312, 3334, 1000010020};
1354
1355 std::map<Int_t,Double_t> mapYields;
1356 std::map<Int_t,Double_t> mapStdDev;
1357 std::map<Int_t,Double_t> mapUncert;
1358 std::map<Int_t,Double_t> mapData;
1359
1360 Double_t data, uncert, model;
1361
1362 ifstream gsiFile(filename);
1363 // std::istream is(thermusFile);
4910214c 1364 std::cout << "----- GSI -----" << std::endl;
b926af7d 1365
1366 // Read the std dev and the ratio in 2 maps, then plot them in a graph.
1367 Int_t ipart = 0;
1368 while(gsiFile) {
1369 gsiFile >> data >> uncert >> model;
1370 if( gsiFile.eof() ) break;
1371 Int_t pdg = pdgPbPb0010[ipart];
1372 std::cout << "PDG " << pdg << " " << data << std::endl;;
1373 mapYields[TMath::Abs(pdg)] += model;
1374 mapUncert[TMath::Abs(pdg)] += uncert;
1375 mapData[TMath::Abs(pdg)] += data;
1376 if(pdg < 0) {
1377 // found the antiparticle: now compute the mean
1378 mapYields[TMath::Abs(pdg)] /=2;
1379 mapData[TMath::Abs(pdg)] /=2;
1380 mapUncert[TMath::Abs(pdg)] /=2;
1381
1382 }
1383 ipart++;
1384 }
1385
1386 // Now plot
1387 TGraphErrors * gGsi = new TGraphErrors;
1388 TGraphErrors * gGsiStdDev = new TGraphErrors;
1389 std::cout << "PDG \tmodel\tdata\tuncert\tstddev" << std::endl; // header
1390 for(Int_t ipart = 0; ipart < npart; ipart++){ // Here we use npart, as this is what we wnat to plot!
1391 Int_t pdg = particleYields[ipart];
1392 mapStdDev[TMath::Abs(pdg)] = ( mapYields[TMath::Abs(pdg)] - mapData[TMath::Abs(pdg)]) / mapUncert[TMath::Abs(pdg)] ;
1393 // chi2 +=
1394 std::cout << "PDG " << pdg <<"\t"
1395 << mapYields[TMath::Abs(pdg)] << "\t" << mapData[TMath::Abs(pdg)] <<"\t"
1396 << mapUncert[TMath::Abs(pdg)] << "\t" << mapStdDev[TMath::Abs(pdg)]
1397 << std::endl;
1398
1399 gGsi->SetPoint(ipart, ipart+1.5, mapYields[particleYields[ipart]]);
1400 gGsi->SetPointError(ipart, 0.3, 0);
1401 gGsiStdDev->SetPoint(ipart, ipart+1.5, mapStdDev[particleYields[ipart]]);
1402 gGsiStdDev->SetPointError(ipart, 0.3, 0);
1403 }
1404
1405 myPadHisto->cd();
1406 gGsi->Draw("PZ");
1407 gGsi->SetLineWidth(lw);
15fccaa0 1408
b926af7d 1409 gGsi->SetLineColor(color);
1410 gGsi->SetLineStyle(lineStyle);
1411 gGsi->SetTitle("NoLegend");
1412
1413
1414 myPadStdDev->cd();
1415 gGsiStdDev->Draw("PZ");
1416 gGsiStdDev->SetLineWidth(lw);
1417 gGsiStdDev->SetLineColor(color);
1418 gGsiStdDev->SetLineStyle(lineStyle);
1419 myPadHisto->cd();
1420
4910214c 1421 AddLineToThermalLegend(gGsi, tag, "L");
b926af7d 1422
1423 return gGsi;
1424}
1425
1426
1427void DrawExtrapolatedSymbolsAndLegendPbPb0010() {
1428 myPadHisto->cd();
1429
1430
1431 TLegend *leg = new TLegend( 0.149598, 0.782203, 0.415663, 0.858447,NULL,"pNDC");
1432 leg->SetBorderSize(0);
1433 leg->SetTextFont(43);
1434 leg->SetTextSize(18);
1435 leg->SetLineColor(1);
1436 leg->SetLineStyle(1);
1437 leg->SetLineWidth(2);
1438 leg->SetFillColor(0);
1439 leg->SetFillStyle(1001);
1440 leg->SetMargin(0.1);
1441
1442 TLegendEntry *entry=leg->AddEntry("TMarker","Extrapolated (Pb-Pb 0-10%)","p");
1443 entry->SetLineColor(1);
1444 entry->SetLineStyle(1);
1445 entry->SetLineWidth(1);
1446 entry->SetMarkerColor(1);
1447 entry->SetMarkerStyle(27);
1448 entry->SetMarkerSize(1.2);
1449 entry->SetTextFont(43);
1450 leg->Draw();
1451 myPadLabel->cd();
1452 // Markers for extrapolated points
1453 TMarker *marker = new TMarker(0.666667,0.111825,27);
1454 marker->SetMarkerStyle(27);
1455 marker->SetMarkerSize(1.2);
1456 marker->Draw();
1457 marker = new TMarker(0.920683,0.0904227,27);
1458 marker->SetMarkerStyle(27);
1459 marker->SetMarkerSize(1.2);
1460 marker->Draw();
1461
1462}
1463void DrawExtrapolatedSymbolsAndLegendpPb0005() {
1464 myPadHisto->cd();
1465
1466
1467 TLegend *leg = new TLegend( 0.149598, 0.709972, 0.415663, 0.786216,NULL,"pNDC");
1468 leg->SetBorderSize(0);
1469 leg->SetTextFont(43);
1470 leg->SetTextSize(18);
1471 leg->SetLineColor(1);
1472 leg->SetLineStyle(1);
1473 leg->SetLineWidth(2);
1474 leg->SetFillColor(0);
1475 leg->SetFillStyle(1001);
1476 leg->SetMargin(0.1);
1477 TLegendEntry *entry=leg->AddEntry("TMarker","Extrapolated (p-Pb 0-5%)","p");
1478 entry->SetLineColor(1);
1479 entry->SetLineStyle(1);
1480 entry->SetLineWidth(1);
1481 entry->SetMarkerColor(1);
1482 entry->SetMarkerStyle(28);
1483 entry->SetMarkerSize(1.2);
1484 entry->SetTextFont(43);
1485 leg->Draw();
1486 myPadLabel->cd();
1487
1488 TMarker *marker = new TMarker(0.590361,0.0797218,28);
1489 marker->SetMarkerStyle(28);
1490 marker->SetMarkerSize(1.2);
1491 marker->Draw();
1492 marker = new TMarker(0.938755,0.0797218,28);
1493 marker->SetMarkerStyle(28);
1494 marker->SetMarkerSize(1.2);
1495 marker->Draw();
1496}
4910214c 1497
1498void AddLineToThermalLegend(TObject * obj, TString line, const char * optFirst) {
1499
1500 // This should be a comma-separated list of text to be added to the
1501 // columns. If the number of entries does not match the numer of
1502 // columns, it returns an error
1503 TObjArray * tokens = line.Tokenize(",");
1504 if(tokens->GetEntries() != legThermal->GetNColumns()) {
1505 std::cout << "Wrong number of columns (" << tokens->GetEntries() << ","<<legThermal->GetNColumns()<<") for the thermal legend, not adding " << line.Data() << std::endl;
1506 return;
1507 }
1508
1509 TIter iter(tokens);
1510 TObjString * col;
1511 Bool_t first = 1;
1512 while((col = (TObjString*)iter.Next())) {
1513 // Add entry removing whitespaces
1514 legThermal->AddEntry(obj, col->String().Strip(TString::kBoth, ' ').Data(), first ? optFirst : "0");
1515 if (first) first = 0;
1516 }
1517}
1518
1519void DrawExtrapolatedSymbolsYieldsPbPb0010(){
1520 // Markers for extrapolated points
1521 myPadLabel->cd();
1522 TMarker * marker = new TMarker(0.369478,0.111825,markerExtrap);
1523 marker->SetMarkerStyle(markerExtrap);
1524 marker->SetMarkerSize(1.2);
1525 marker->Draw();
1526 marker = new TMarker(0.938755,0.111825,markerExtrap);
1527 marker->SetMarkerStyle(markerExtrap);
1528 marker->SetMarkerSize(1.2);
1529 marker->Draw();
1530 myPadHisto->cd();
1531}
1532
1533void DrawMarkerKStarNoFit() {
1534 myPadLabel->cd();
1535 TMarker *marker = new TMarker(0.339357,0.111825,markerNoFit);
1536 marker->SetMarkerStyle(28);
1537 marker->SetMarkerSize(1.2);
1538 marker->Draw();
1539 myPadHisto->cd();
1540}