]> 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
Macros improvements
[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"
7#include <fstream>
8#endif
9
10enum 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};
11
ad431d97 12TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ;
13TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) ;
9d84731e 14void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
15
16// Plots ratios for QM and saves input files for thermal models
17
18TClonesArray * PlotRatiosForQM14() {
ad431d97 19#if !(!defined (__CINT__) || (defined(__MAKECINT__)))
20 LoadLibs();
21#endif
9d84731e 22 TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
23 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
24 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
25 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
26 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
27 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
28 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
29 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
30
ad431d97 31 TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
32
33 TClonesArray * arrpp276 = AliParticleYield::ReadFromASCIIFile("pp_2760.txt");
34 TClonesArray * arrpp900 = AliParticleYield::ReadFromASCIIFile("pp_900.txt");
35
36 TClonesArray * arrpPb = AliParticleYield::ReadFromASCIIFile("pPb_5020_MultiStrange.txt");
37 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_PiKaPrLamndaK0.txt"));
38 arrpPb->AbsorbObjects(AliParticleYield::ReadFromASCIIFile("pPb_5020_deuteron.txt"));
39
40 TClonesArray * arrThermus = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Thermus_Boris_20140407.txt");
9d84731e 41
ad431d97 42 // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges0);
43 // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", separateCharges1);
44 // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/0);
45 // PrepareThermalModelsInputFiles(arrpp7, AliParticleYield::kCSpp, 7000, "", /*separateCharges*/1);
9d84731e 46
47 TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
ad431d97 48 c1->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368);
49
50 // c1->SetLogy();
51 // CENTRAL
52 TH1 * h = GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%");
53 h->GetYaxis()->SetRangeUser(0, 0.4);
54 h->Draw();
55 // //GetHistoRatios(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
56 // GetHistoRatios(arrpp7, AliParticleYield::kCSpp, 7000, "" , "pp #sqrt{s} = 7 TeV" )->Draw("same");
57 // GetHistoRatios(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same");
58 // // GetHistoRatios(arrpp276, AliParticleYield::kCSpp, 2760, "" , "pp #sqrt{s} = 2.76 TeV" )->Draw("same");
59 // // GetHistoRatios(arrpp900, AliParticleYield::kCSpp, 900, "" , "pp #sqrt{s} = 0.9 TeV" )->Draw("same");
60 // NewLegend("", "lp",0,1,1);
61
62 // Peripheral
63 GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M6080", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 60-80%")->Draw("same");
64 //GetHistoRatios(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
65 // GetHistoRatios(arrpp7, AliParticleYield::kCSpp, 7000, "" , "pp #sqrt{s} = 7 TeV" )->Draw("same");
66 GetHistoRatios(arrpPb, AliParticleYield::kCSpPb, 5020, "V0A0005", "p-Pb, #sqrt{s_{NN}} = 5.02 TeV, V0A 0-5%")->Draw("same");
67 // GetHistoRatios(arrpp276, AliParticleYield::kCSpp, 2760, "" , "pp #sqrt{s} = 2.76 TeV" )->Draw("same");
68 // GetHistoRatios(arrpp900, AliParticleYield::kCSpp, 900, "" , "pp #sqrt{s} = 0.9 TeV" )->Draw("same");
69 NewLegend("", "lp",0,1,1);
70
71 //return;
72 TCanvas * c2 = new TCanvas("Yields", "Yields", 1400, 600);
73 c2->SetMargin( 0.0744986, 0.0329513, 0.225131, 0.0593368);
74
75 GetHistoYields(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Pb-Pb, #sqrt{s_{NN}} = 2.76 TeV, 0-10%")->Draw();
76 GetHistoYields(arrThermus, AliParticleYield::kCSPbPb, 2760, "V0M0010", "Thermus")->Draw("same");
77 NewLegend("", "lp",0,1,1);
9d84731e 78 return arr;
79}
80
ad431d97 81TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) {
82 // FIXME: THIS SHOULD BE REVIEWED TO MAKE SURE THE PLOTS ARE LABELLED CORRECTLY
9d84731e 83
84 const Int_t nratio = 10;
ad431d97 85 Int_t num [nratio] = {kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar};
86 // Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGPi , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK};
87 Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGProton , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK};
88 Int_t isSum[nratio] = {1 ,1 ,0 ,1 ,1 ,0 ,0 ,1 ,0 ,1 };
89 Double_t scale[nratio] = {1 ,3 ,1 ,30 ,250 ,50 ,100 ,4e5 ,2 ,2 };
90 TH1F * h = new TH1F(Form("hRatio_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, nratio, 1, nratio+1);
9d84731e 91
92 // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
93 for(Int_t iratio = 1; iratio <= nratio; iratio++){
94 AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
95 std::cout << num[iratio-1] << " " << denum[iratio-1]<< " " ;
96 if(ratio)ratio->Print("short");
97
98
99 if(!ratio) {
100 // If the ratio is not found, try to build it!
101 AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality, isSum[iratio-1]);
ad431d97 102 // Try with the !sum, if part 1 is not found
103 if(!part1) {
9d84731e 104 part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
ad431d97 105 // If the sum was requested, try to recover it!
106 if(isSum[iratio-1]) {
107 AliParticleYield * part1_bar = AliParticleYield::FindParticle(arr, -num[iratio-1], system, energy, centrality,0);
108 if(part1 && part1_bar) {
109 std::cout << "Adding " << part1_bar->GetPartName() << " " << part1->GetPartName() << std::endl;
110
111 part1 = AliParticleYield::Add(part1, part1_bar);
112
113 }
114 } else if(part1) {
115 // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
116 part1->Scale(0.5);
117 }
118
9d84731e 119 }
120 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
121 if(!part2) {// Try with the !sum, if part 2 is not found
122 part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
ad431d97 123 if(isSum[iratio-1]) {
124 AliParticleYield * part2_bar = AliParticleYield::FindParticle(arr, -denum[iratio-1], system, energy, centrality,0);
125 if(part2 && part2_bar) part2 = AliParticleYield::Add(part2, part2_bar);
126 } else if(part2){
127 // if the not sum was requested, but the sum is found, divide by 2 so that it is comparable
128 part2->Scale(0.5);
129 }
130
9d84731e 131 }
ad431d97 132 ratio = AliParticleYield::Divide(part1, part2, 0, "YQ");
9d84731e 133 if(ratio) {
ad431d97 134 std::cout << "" << std::endl;
9d84731e 135 std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
ad431d97 136 ratio->Print("");
137 std::cout << "" << std::endl;
9d84731e 138 }
139 }
140 if(ratio){
ad431d97 141 ratio->Scale(scale[iratio-1]);
9d84731e 142 h->SetBinContent(iratio, ratio->GetYield());
143 h->SetBinError (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
ad431d97 144 h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",Form("#times%2.2f", scale[iratio-1]), ratio->GetLatexName()));
9d84731e 145 }
146 else {
147 h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
148
149 }
150 }
151
152 return h;
153
154
155
156}
157
158void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges) {
159 // If "Separate charges" is true, tries to dump both charges are dumped
160 TClonesArray * arrOut = new TClonesArray("AliParticleYield");
161 const Int_t npart = 12;
162 Int_t particles [npart] = {kPDGPi ,kPDGK ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3 };
163 Int_t isSum[npart] = {1 ,1 ,0 , 1 , 0 , 1 ,0 ,1 ,1 ,0 , 1 , 0 };
164
165 Int_t ipartOut = 0; // Index for the array
166 for(Int_t ipart = 0; ipart < npart; ipart++){
167 if(!separateCharges) {
168 AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, isSum[ipart]);
169 if(!part && isSum[ipart]) {
170 //Could not find the particle, but the sum was requested: build the sum!
171 part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
172 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
173 if(part2 && part) part = AliParticleYield::Add(part, part2);
174 else part = 0;
175 }
176 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
177 }
178 else {
179 // ignore isSum and try to find both particles
180 Bool_t notFound = 0;
181 AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
182 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
183 else notFound=1;
184 // Try to find antiparticle (-pdg code)
185 part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
186 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
187 else if (notFound) {
188 // If neither charge was found, check if we at least have the sum
189 part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality, 1);
190 if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
191 }
192
193 }
194 }
195 std::cout << "Particles for thermal model fits:" << std::endl;
196 arrOut->Print("short");
197 std::cout << "" << std::endl;
198 // Write GSI input file
199 TIter it(arrOut);
200 AliParticleYield * part = 0;
201 ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
202 while ((part = (AliParticleYield*) it.Next())){
203 fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
204 }
205 fout.close();
206 // Write thermus file
207 AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
208}
ad431d97 209
210
211TH1F * GetHistoYields(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, const char * histotitle) {
212
213 const Int_t npart = 11;
214 Int_t pdg [npart] = {kPDGPi, kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar};
215 // Int_t isSum[npart] = {1 ,1 ,1 ,0 ,1 ,1 ,0 ,0 ,1 ,0 ,1 };
216 Int_t isSum[npart] = {0,0,0,0,0,0,0,0,0,0,1};
217 // Double_t scale[npart] = {1 ,1 ,3 ,1 ,30 ,250 ,50 ,10 ,4e5 ,2 ,2 };
218 Double_t scale[npart] = {1,5,30,30,200,1000,4000,2e6,2e6,20,20,};
219 TH1F * h = new TH1F(Form("hPart_%d_%0.0f_%s",system,energy,centrality.Data()), histotitle, npart, 1, npart+1);
220
221 // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
222 for(Int_t ipart = 1; ipart <= npart; ipart++){
223 AliParticleYield * part = AliParticleYield::FindParticle(arr, pdg[ipart-1], system, energy, centrality,isSum[ipart-1]);
224 if(!part) continue;
225 part->Scale(scale[ipart-1]);
226 h->SetBinContent(ipart, part->GetYield());
227 h->SetBinError (ipart, part->GetTotalError(0/* 0 = no normalization error */));
228 h->GetXaxis()->SetBinLabel(ipart, Form("#splitline{%s}{%s}",Form("#times%2.0g", scale[ipart-1]), part->GetLatexName()));
229
230 }
231 return h;
232}