B2 analysis code
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / B2 / macros / B2Mult.C
CommitLineData
2403d402 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16// B2 as a function of multiplicity
17// author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
18
19#include <Riostream.h>
20#include <TSystem.h>
21#include <TROOT.h>
22#include <TFileMerger.h>
23#include <TString.h>
24#include <TFile.h>
25#include <TGraphErrors.h>
26
27#include "AliLnB2.h"
28#include "B2.h"
29#include "Config.h"
30
31Double_t GetCd(Double_t z)
32{
33//
34// parameterization of <Cd> as a function of multiplicity
35// from ALICE Rlong and Rside measurements
36//
37 return 0.046133 + 0.0484458*z;
38}
39
40Int_t B2Mult(const TString& pSpectra = "~/alice/output/Proton-lhc10d-Mult-Spectra.root",
41 const TString& ptag = "lhc10d",
42 const TString& dSpectra = "~/alice/output/Deuteron-lhc10d-Mult-Spectra.root",
43 const TString& dtag = "lhc10d",
44 const TString& outputMultPt = "~/alice/output/B2-Mult-Pt.root",
45 const TString& outputPtMult = "~/alice/output/B2-Pt-Mult.root",
46 const TString& otag = "lhc10d")
47{
48//
49// B2 as a function of multiplicity
50//
51 using namespace B2mult;
52
53 const Int_t kNpart = 2;
54
55 const Int_t kNMinPt = 0;
56 const Int_t kNMaxPt = 6;
57
58 const TString kPrefix[] = { "", "Anti"};
59 const TString kSuffix[] = { "", "bar" };
60 Int_t kCharge[] = {1, -1};
61
62 // B2 as a function of pt for each multiplicity class
63
64 for(Int_t i=0; i<kNmult; ++i)
65 {
66 TFileMerger m;
67
68 for(Int_t j=0; j<kNpart; ++j)
69 {
70 TString b2file = kPrefix[j] + "B2.root";
71
72 AliLnB2 b2(pSpectra, ptag + "-" + kMultClass[i], dSpectra, dtag + "-" + kMultClass[i], b2file, otag + "-" + kMultClass[i], 2, kCharge[j]);
73
74 b2.SetCd(GetCd(kKNOmult[i]));
75
76 b2.Run();
77
78 m.AddFile(b2file.Data(),0);
79 }
80
81 // merge B2 and B2bar
82
83 TString outputfile = otag + "-" + kMultClass[i] + "-B2.root";
84
85 m.OutputFile(outputfile.Data());
86 m.Merge();
87
88 gSystem->Exec("rm -f B2.root AntiB2.root");
89 }
90
91 // merge multiplicity classes
92
93 TFileMerger m;
94
95 for(Int_t i=0; i<kNmult; ++i)
96 {
97 TString b2 = otag + "-" + kMultClass[i] + "-B2.root";
98 m.AddFile(b2.Data(),0);
99 }
100
101 m.OutputFile(outputMultPt.Data());
102 m.Merge();
103
104 // delete tmp files
105
106 for(Int_t i=0; i<kNmult; ++i)
107 {
108 gSystem->Exec(Form("rm -f %s-%s-B2.root",otag.Data(),kMultClass[i].Data()));
109 }
110
111 // B2 as a function of multiplicity for each pt
112
113 TFile* finput = new TFile(outputMultPt.Data());
114 if (finput->IsZombie()) exit(1);
115
116 TGraphErrors* grB2pt[kNpart][kNmult];
117 TGraphErrors* grR3pt[kNpart][kNmult];
118
119 for(Int_t i=0; i<kNpart; ++i)
120 {
121 for(Int_t j=0; j<kNmult; ++j)
122 {
123 grB2pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("B2%s_Pt", kSuffix[i].Data()));
124 grR3pt[i][j] = (TGraphErrors*)FindObj(finput, otag + "-" + kMultClass[j], Form("R3%s_Pt", kSuffix[i].Data()));
125 }
126 }
127
128 TFile* foutput = new TFile(outputPtMult.Data(),"recreate");
129
130 Double_t* pt = grB2pt[0][0]->GetX();
131 TString ptLabel[kNMaxPt];
132
133 if(kNMaxPt > grB2pt[0][0]->GetN())
134 {
135 std::cerr << "max pt too big" << std::endl;
136 exit(1);
137 }
138
139 for(Int_t i=kNMinPt; i<kNMaxPt; ++i)
140 {
141 ptLabel[i] = Form("pt%.02fA",pt[i]);
142 foutput->mkdir(ptLabel[i].Data());
143 }
144
145 for(Int_t i=0; i<kNpart; ++i)
146 {
147 for(Int_t j=kNMinPt; j<kNMaxPt; ++j)
148 {
149 Double_t B2[kNmult];
150 Double_t B2StatErr[kNmult];
151
152 Double_t R3[kNmult];
153 Double_t R3StatErr[kNmult];
154
155 for(Int_t k=0; k<kNmult; ++k)
156 {
157 Double_t x, y, ey;
158 grB2pt[i][k]->GetPoint(j,x,y);
159 ey = grB2pt[i][k]->GetErrorY(j);
160
161 B2[k] = y;
162 B2StatErr[k] = ey;
163
164 grR3pt[i][k]->GetPoint(j,x,y);
165 ey = grR3pt[i][k]->GetErrorY(j);
166
167 R3[k] = y;
168 R3StatErr[k] = ey;
169 }
170
171 TGraphErrors* grB2Mult = new TGraphErrors(kNmult, kKNOmult, B2, kKNOmultErr, B2StatErr);
172 grB2Mult->SetName(Form("B2%s_Zmult", kSuffix[i].Data()));
173
174 TGraphErrors* grR3Mult = new TGraphErrors(kNmult, kKNOmult, R3, kKNOmultErr, R3StatErr);
175 grR3Mult->SetName(Form("R3%s_Zmult", kSuffix[i].Data()));
176
177 foutput->cd(ptLabel[j].Data());
178
179 grB2Mult->Write();
180 grR3Mult->Write();
181
182 delete grB2Mult;
183 delete grR3Mult;
184 }
185 }
186
187 delete foutput;
188 delete finput;
189
190 // particle ratios
191
192 gROOT->ProcessLine(Form(".x RatioMult.C+g(\"%s\",\"%s\",\"%s\",\"%s\")",pSpectra.Data(), ptag.Data(),dSpectra.Data(), dtag.Data()));
193
194 // draw B2 as a function of pt
195
196 for(Int_t i=0; i<kNpart; ++i)
197 {
198 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Pt\",\"\",0,2, 1.e-3, 7.e-2,\"p_{T}/A (GeV/c)\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2pt\",\"B2%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
199
200 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Pt\",\"\",0,2, 0, 1.7,\"p_{T}/A (GeV/c)\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3pt\",\"R3%spt\")", outputMultPt.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
201 }
202
203 // draw B2 as a function of z
204
205 for(Int_t i=0; i<kNpart; ++i)
206 {
207 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"B2%s_Zmult\",\"\",0,7, 1.e-3, 7.e-2,\"z\",\"B_{2} (GeV^{2}/c^{3})\", 0,\"c%d.B2z\",\"B2%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
208
209 gROOT->ProcessLine(Form(".x DrawDir.C+g(\"%s\",\"R3%s_Zmult\",\"\",0,8, 0, 4,\"z\",\"R_{side}^{2} R_{long} (fm^{3})\", 0,\"c%d.R3z\",\"R3%sZ\")", outputPtMult.Data(), kSuffix[i].Data(), i, kSuffix[i].Data()));
210 }
211
212 return 0;
213}