]>
Commit | Line | Data |
---|---|---|
2b893216 | 1 | //____________________________________________________________________ |
2 | // | |
3 | // $Id$ | |
4 | // | |
5 | // Script to get the various cross sections, energy loss, and ranges | |
6 | // of a particular particle type in a particular medium. | |
7 | // | |
8 | // This script should be compiled to speed it up. | |
9 | // | |
10 | // It creates a tree on the current output file, with the relevant | |
11 | // information. | |
12 | // | |
13 | // Note, that VMC _must_ be the TGeant3TGeo VMC. | |
14 | // | |
15 | #include <TArrayF.h> | |
16 | #include <TTree.h> | |
17 | #include <TMath.h> | |
18 | #include <iostream> | |
19 | #include <TVirtualMC.h> | |
20 | #include <TDatabasePDG.h> | |
21 | #include <TString.h> | |
22 | #include <TGeoManager.h> | |
23 | #include <TGeoMedium.h> | |
24 | #include <TGeant3.h> | |
25 | #include <TGraph.h> | |
26 | #include <TGraphErrors.h> | |
27 | #include <TAxis.h> | |
28 | #include <TFile.h> | |
29 | #include <TH2.h> | |
30 | #include <TLegend.h> | |
31 | #include <AliRun.h> | |
32 | #include <TCanvas.h> | |
33 | #include <cmath> | |
34 | #include <string> | |
35 | #include <vector> | |
36 | #include <sstream> | |
37 | #include <iomanip> | |
38 | /** @defgroup xsec_script X-section script | |
39 | @ingroup FMD_script | |
40 | */ | |
41 | //____________________________________________________________________ | |
42 | /** @ingroup xsec_script | |
43 | */ | |
44 | struct Mech | |
45 | { | |
46 | char* fName; | |
47 | char* fTitle; | |
48 | char* fUnit; | |
49 | int fColor; | |
50 | int fStyle; | |
51 | }; | |
52 | ||
53 | Mech fgMechs[] = | |
54 | {{ "HADF","FLUKA total hadronic x-section", "cm^{-1}",1,1}, | |
55 | { "INEF","FLUKA hadronic inelastic x-section", "cm^{-1}",2,1}, | |
56 | { "ELAF","FLUKA hadronic elastic x-section", "cm^{-1}",3,1}, | |
57 | { "HADG","GHEISHA total hadronic x-section", "cm^{-1}",4,1}, | |
58 | { "INEG","GHEISHA hadronic inelastic x-section", "cm^{-1}",6,1}, | |
59 | { "ELAG","GHEISHA hadronic elastic x-section", "cm^{-1}",7,1}, | |
60 | { "FISG","GHEISHA nuclear fission x-section", "cm^{-1}",1,2}, | |
61 | { "CAPG","GHEISHA neutron capture x-section", "cm^{-1}",2,2}, | |
62 | { "LOSS","stopping power", "MeV/cm", 3,2}, | |
63 | { "PHOT","photoelectric x-section", "cm^{-1}",4,2}, | |
64 | { "ANNI","positron annihilation x-section", "cm^{-1}",6,2}, | |
65 | { "COMP","Compton effect x-section", "cm^{-1}",7,2}, | |
66 | { "BREM","bremsstrahlung x-section", "cm^{-1}",1,3}, | |
67 | { "PAIR","photon and muon direct- pair x-section","cm^{-1}",2,3}, | |
68 | { "DRAY","delta-rays x-section", "cm^{-1}",3,3}, | |
69 | { "PFIS","photo-fission x-section", "cm^{-1}",4,3}, | |
70 | { "RAYL","Rayleigh scattering x-section", "cm^{-1}",6,3}, | |
71 | { "MUNU","muon-nuclear interaction x-section", "cm^{-1}",7,3}, | |
72 | { "RANG","range", "cm", 1,4}, | |
73 | { "STEP","maximum step", "cm", 2,4}, | |
74 | { 0, 0, 0, 0, 0}}; | |
75 | ||
76 | //____________________________________________________________________ | |
77 | /** @ingroup xsec_script | |
78 | */ | |
79 | struct MechValue | |
80 | { | |
81 | MechValue(const Mech& m, size_t n) | |
82 | : fMech(m), fValues(n), fStatus(0), fELoss(0) | |
83 | { | |
84 | fGraph.SetName(fMech.fName); | |
85 | fGraph.SetTitle(fMech.fTitle); | |
86 | fGraph.GetXaxis()->SetTitle("#beta#gamma"); | |
87 | fGraph.GetYaxis()->SetTitle(fMech.fUnit); | |
88 | fGraph.SetLineColor(fMech.fColor); | |
89 | fGraph.SetLineStyle(fMech.fStyle); | |
90 | fGraph.SetLineWidth(2); | |
91 | std::string name(fMech.fName); | |
92 | if (name != "LOSS") return; | |
93 | ||
94 | fELoss = new TGraphErrors(n); | |
95 | fELoss->SetName("ELOSS"); | |
96 | fELoss->SetTitle(fMech.fTitle); | |
97 | fELoss->GetXaxis()->SetTitle("#beta#gamma"); | |
98 | fELoss->GetYaxis()->SetTitle(fMech.fUnit); | |
99 | fELoss->SetLineColor(fMech.fColor); | |
100 | fELoss->SetLineStyle(fMech.fStyle); | |
101 | fELoss->SetLineWidth(2); | |
102 | fELoss->SetFillColor(fMech.fColor+1); | |
103 | fELoss->SetFillStyle(3001); | |
104 | } | |
105 | bool Get(TGeant3* mc, std::vector<float>& t, | |
106 | std::vector<float>& c, int medNo, int pidNo, | |
107 | float mass) | |
108 | { | |
109 | int ixst; | |
110 | mc->Gftmat(medNo, pidNo, fMech.fName, t.size(), | |
111 | &(t[0]), &(fValues[0]), &(c[0]), ixst); | |
112 | fStatus = ixst; | |
113 | if (!fStatus) return false; | |
114 | fGraph.Set(t.size()); | |
115 | for (size_t i = 0; i < t.size(); i++) { | |
116 | fGraph.SetPoint(i, t[i]/mass, fValues[i]); | |
117 | if (!fELoss) continue; | |
118 | fELoss->SetPoint(i, t[i]/mass, fValues[i]); | |
119 | // ~ 5 sigma | |
120 | fELoss->SetPointError(i, 0, .5 * fValues[i]); | |
121 | } | |
122 | fGraph.Write(); | |
123 | if (fELoss) fELoss->Write(); | |
124 | return true; | |
125 | } | |
126 | TGraph& Draw() | |
127 | { | |
128 | if (fELoss) fELoss->DrawClone("4 same"); | |
129 | fGraph.DrawClone("l same"); | |
130 | return fGraph; | |
131 | } | |
132 | const Mech& fMech; | |
133 | std::vector<float> fValues; | |
134 | bool fStatus; | |
135 | TGraph fGraph; | |
136 | TGraphErrors* fELoss; | |
137 | }; | |
138 | ||
139 | //____________________________________________________________________ | |
140 | /** @ingroup xsec_script | |
141 | */ | |
142 | struct XSections | |
143 | { | |
144 | XSections(size_t n=91, float emin=1e-5, float emax=1e4) | |
145 | : fTKine(n), | |
146 | fCuts(5) | |
147 | { | |
148 | float dp = 1. / log10(emax/emin); | |
149 | float pmin = log10(emin); | |
150 | fTKine[0] = emin; | |
151 | for (size_t i = 1; i < fTKine.size(); i++) { | |
152 | float el = pmin + i * dp; | |
153 | fTKine[i] = pow(10, el); | |
154 | } | |
155 | for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i) | |
156 | *i = 1e-4; | |
157 | Mech* mech = &(fgMechs[0]); | |
158 | while (mech->fName) { | |
159 | fMechs.push_back(new MechValue(*mech, n)); | |
160 | mech++; | |
161 | } | |
162 | } | |
163 | void Run(const std::string& medName, const std::string& pdgName) | |
164 | { | |
165 | TGeant3* mc = static_cast<TGeant3*>(gMC); | |
166 | if (!mc) { | |
167 | std::cerr << "Couldn't get VMC" << std::endl; | |
168 | return; | |
169 | } | |
170 | TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str()); | |
171 | if (!medium) { | |
172 | std::cerr << "Couldn't find medium " << medName << std::endl; | |
173 | return; | |
174 | } | |
175 | int medNo = medium->GetMaterial()->GetUniqueID(); | |
176 | TDatabasePDG* pdgDb = TDatabasePDG::Instance(); | |
177 | fPDG = pdgDb->GetParticle(pdgName.c_str()); | |
178 | if (!fPDG) { | |
179 | std::cerr << "Couldn't find particle " << pdgName << std::endl; | |
180 | return; | |
181 | } | |
182 | int pdgNo = fPDG->PdgCode(); | |
183 | int pidNo = mc->IdFromPDG(pdgNo); | |
184 | ||
185 | std::stringstream vars; | |
186 | vars << "betagamma/F"; | |
187 | ||
188 | size_t nOk = 0; | |
189 | // Loop over defined mechanisms | |
190 | for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) { | |
191 | if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo, fPDG->Mass()))continue; | |
192 | vars << ":" << (*i)->fMech.fName; | |
193 | nOk ++; | |
194 | } | |
195 | ||
196 | std::stringstream tName; | |
197 | tName << medName << "_" << pdgName; | |
198 | TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str()); | |
199 | ||
200 | float_array cache(nOk+1); | |
201 | tree->Branch("xsec", &(cache[0]), vars.str().c_str()); | |
202 | for (size_t i = 0; i < fTKine.size(); i++) { | |
203 | cache[0] = fTKine[i] / fPDG->Mass(); | |
204 | int k = 0; | |
205 | for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) { | |
206 | if (!(*j)->fStatus) continue; | |
207 | cache[++k] = (*j)->fValues[i]; | |
208 | } | |
209 | tree->Fill(); | |
210 | } | |
211 | tree->Write(); | |
212 | } | |
213 | void Draw() | |
214 | { | |
215 | float min = 100000; | |
216 | float max = 0; | |
217 | std::vector<TGraph*> gs; | |
218 | float_array bg(fTKine.size()); | |
219 | for (size_t i = 0; i < fTKine.size(); i++) bg[i] = fTKine[i]/fPDG->Mass(); | |
220 | for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) { | |
221 | if (!(*j)->fStatus) continue; | |
222 | for (size_t i = 0; i < fTKine.size(); i++) { | |
223 | if ((*j)->fValues[i] == 0) continue; | |
224 | min = std::min(min, (*j)->fValues[i]); | |
225 | max = std::max(max, (*j)->fValues[i]); | |
226 | } | |
227 | } | |
228 | ||
229 | TCanvas* c = new TCanvas("c", "C", 700, 700); | |
230 | c->SetFillColor(0); | |
231 | c->SetLogy(); | |
232 | c->SetLogx(); | |
233 | c->SetGridy(); | |
234 | c->SetGridx(); | |
235 | float_array y(101); | |
236 | float ymin = log10(min); | |
237 | float dy = (log10(max)+.5 - log10(min)) / y.size(); | |
238 | for (size_t i = 1; i < y.size(); i++) y[i] = pow(10, ymin + i * dy); | |
239 | TH2* f = new TH2F("x", "X-sec",bg.size()-1,&(bg[0]),y.size()-1,&(y[0])); | |
240 | f->SetXTitle("#beta#gamma"); | |
241 | f->SetDirectory(0); | |
242 | f->SetStats(kFALSE); | |
243 | f->Draw(); | |
244 | TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45); | |
245 | l->SetFillColor(0); | |
246 | // l->SetFillStyle(0); | |
247 | for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) { | |
248 | if (!(*j)->fStatus) continue; | |
249 | TGraph& g = (*j)->Draw(); | |
250 | l->AddEntry(&g, g.GetTitle(), "l"); | |
251 | } | |
252 | l->Draw("same"); | |
253 | } | |
254 | protected: | |
255 | typedef std::vector<float> float_array; | |
256 | typedef std::vector<MechValue*> mech_array; | |
257 | float_array fTKine; | |
258 | float_array fCuts; | |
259 | mech_array fMechs; | |
260 | TParticlePDG* fPDG; | |
261 | }; | |
262 | ||
263 | bool init = false; | |
264 | ||
265 | void XSection(const char* pdgName="pi-") | |
266 | { | |
267 | if (!init) { | |
268 | gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C"); | |
269 | init = true; | |
270 | } | |
271 | TFile* file = TFile::Open("xsec.root", "RECREATE"); | |
272 | XSections xs; | |
273 | xs.Run("FMD_Si$", pdgName); | |
274 | xs.Draw(); | |
275 | file->Close(); | |
276 | } | |
277 | //____________________________________________________________________ | |
278 | // | |
279 | // EOF | |
280 | // |