1 #if !defined (__CINT__) || (defined(__MAKECINT__))
3 #include "TClonesArray.h"
4 #include "AliParticleYield.h"
10 enum MyParticles { kPDGPi = 211, kPDGK = 321, kPDGProton = 2212, kPDGKS0 = 310, kPDGLambda=3122, kPDGXi=3312,kPDGOmega=3334,kPDGPhi=333,kPDGKStar=313,kPDGDeuteron=1000010020,kPDGHE3 = 1000020030, kPDGHyperTriton = 1010010030, kPDGSigmaStarPlus=3224,kPDGSigmaStarMinus=3114,kPDGXiStar=3324};
12 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) ;
13 void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges=0) ;
15 // Plots ratios for QM and saves input files for thermal models
17 TClonesArray * PlotRatiosForQM14() {
19 TClonesArray * arr = AliParticleYield::ReadFromASCIIFile("PbPb_2760_Cascades.txt");
20 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_DeuHelium3.txt"));
21 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Hypertriton.txt"));
22 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_Kstar892.txt"));
23 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_LambdaK0.txt"));
24 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_PiKaPr.txt"));
25 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_phi1020.txt"));
26 arr->AbsorbObjects( AliParticleYield::ReadFromASCIIFile("PbPb_2760_AveragedNumbers.txt"));
28 TClonesArray * arrpp7 = AliParticleYield::ReadFromASCIIFile("pp_7000.txt");
30 // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/0);
31 // PrepareThermalModelsInputFiles(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010", /*separateCharges*/1);
33 TCanvas * c1 = new TCanvas("Ratios", "Ratios", 1400, 600);
35 // GetHistoRatios(arr, AliParticleYield::kCSPbPb, 2760, "V0M0010")->Draw();
36 GetHistoRatios(arr, AliParticleYield::kCSpp, 7000, ".*")->Draw("");
40 TH1F * GetHistoRatios(TClonesArray * arr, Int_t system, Float_t energy, TString centrality) {
42 const Int_t nratio = 10;
43 Int_t num [nratio] = {kPDGK , kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron , kPDGHE3 , kPDGHyperTriton , kPDGPhi , kPDGKStar};
44 Int_t denum[nratio] = {kPDGPi , kPDGPi , kPDGKS0 , kPDGPi , kPDGPi , kPDGProton , kPDGDeuteron , kPDGPi , kPDGK , -kPDGK};
45 Int_t isSum[nratio] = {1 ,1 ,0 ,1 ,1 ,0 ,0 ,1 ,0 ,1};
46 TH1F * h = new TH1F("hRatio", "hRatio", nratio, 1, nratio+1);
48 // Double_t isSum = -1; // if this is -1, then the sum criterion is ignored
49 for(Int_t iratio = 1; iratio <= nratio; iratio++){
50 AliParticleYield * ratio = AliParticleYield::FindRatio(arr, num[iratio-1], denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
51 std::cout << num[iratio-1] << " " << denum[iratio-1]<< " " ;
52 if(ratio)ratio->Print("short");
56 // If the ratio is not found, try to build it!
57 AliParticleYield * part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality, isSum[iratio-1]);
58 if(!part1) {// Try with the !sum, if part 1 is not found
59 part1 = AliParticleYield::FindParticle(arr, num[iratio-1], system, energy, centrality,!isSum[iratio-1]);
61 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,isSum[iratio-1]);
62 if(!part2) {// Try with the !sum, if part 2 is not found
63 part2 = AliParticleYield::FindParticle(arr, denum[iratio-1], system, energy, centrality,!isSum[iratio-1]);
65 ratio = AliParticleYield::Divide(part1, part2);
67 std::cout << "WARNING: building ratio " << num[iratio-1] <<"/"<<denum[iratio-1]<<": Check uncertainties!!" << std::endl;
68 ratio->Print("short");
72 h->SetBinContent(iratio, ratio->GetYield());
73 h->SetBinError (iratio, ratio->GetTotalError(0/* 0 = no normalization error */));
74 h->GetXaxis()->SetBinLabel(iratio, Form("#splitline{%s}{%s}",ratio->GetCentr().Data() , ratio->GetLatexName()));
77 h->GetXaxis()->SetBinLabel(iratio, Form("#frac{%d}{%d}",num[iratio-1], denum[iratio-1]));
88 void PrepareThermalModelsInputFiles(TClonesArray * arr, Int_t system, Float_t energy, TString centrality, Bool_t separateCharges) {
89 // If "Separate charges" is true, tries to dump both charges are dumped
90 TClonesArray * arrOut = new TClonesArray("AliParticleYield");
91 const Int_t npart = 12;
92 Int_t particles [npart] = {kPDGPi ,kPDGK ,kPDGKS0, kPDGKStar, kPDGPhi, kPDGProton , kPDGLambda , kPDGXi , kPDGOmega , kPDGDeuteron, kPDGHyperTriton, kPDGHE3 };
93 Int_t isSum[npart] = {1 ,1 ,0 , 1 , 0 , 1 ,0 ,1 ,1 ,0 , 1 , 0 };
95 Int_t ipartOut = 0; // Index for the array
96 for(Int_t ipart = 0; ipart < npart; ipart++){
97 if(!separateCharges) {
98 AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, isSum[ipart]);
99 if(!part && isSum[ipart]) {
100 //Could not find the particle, but the sum was requested: build the sum!
101 part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
102 AliParticleYield * part2 = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
103 if(part2 && part) part = AliParticleYield::Add(part, part2);
106 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
109 // ignore isSum and try to find both particles
111 AliParticleYield * part = AliParticleYield::FindParticle(arr, particles[ipart], system, energy, centrality, 0);
112 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
114 // Try to find antiparticle (-pdg code)
115 part = AliParticleYield::FindParticle(arr, -particles[ipart], system, energy, centrality, 0);
116 if(part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
118 // If neither charge was found, check if we at least have the sum
119 part = AliParticleYield::FindParticle(arr, abs(particles[ipart]), system, energy, centrality, 1);
120 if (part) new((*arrOut)[ipartOut++]) AliParticleYield(*part);
125 std::cout << "Particles for thermal model fits:" << std::endl;
126 arrOut->Print("short");
127 std::cout << "" << std::endl;
128 // Write GSI input file
130 AliParticleYield * part = 0;
131 ofstream fout(Form("gsi_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));
132 while ((part = (AliParticleYield*) it.Next())){
133 fout << part->GetYield() << " " << part->GetTotalError() << std::endl;
136 // Write thermus file
137 AliParticleYield::WriteThermusFile(arrOut, Form("thermus_System_%d_Energy_%0.0f_Centr_%s_BothCharges_%d", system, energy, centrality.Data(), separateCharges));