Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
[u/mrichter/AliRoot.git] / FMD / scripts / GetXsection.C
CommitLineData
a1f80595 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//
2b893216 15#ifdef __CINT__
16XSection()
17{
18 gROOT->ProcessLine(".x Compile.C(\"$(ALICE_ROOT)/FMD/scripts/XSection.C\"");
19 gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
20 TFile* file = TFile::Open("xsec.root", "RECREATE");
21 GetXsection("FMD_Si$", "pi+");
22 file->Close();
23}
24#else
a1f80595 25#include <TArrayF.h>
26#include <TTree.h>
27#include <TMath.h>
28#include <iostream>
29#include <TVirtualMC.h>
30#include <TDatabasePDG.h>
31#include <TString.h>
32#include <TGeoManager.h>
33#include <TGeoMedium.h>
34#include <TGeant3.h>
2b893216 35#include <TGraph.h>
36#include <TAxis.h>
37#include <cmath>
38#include <string>
39#include <vector>
40#include <sstream>
41#include <iomanip>
9b48326f 42/** @defgroup xsec_script X-section script
43 @ingroup FMD_script
44*/
bf000c32 45//____________________________________________________________________
9b48326f 46/** @ingroup xsec_script
47 */
a1f80595 48struct Mech
49{
2b893216 50 char* fName;
51 char* fTitle;
52 char* fUnit;
a1f80595 53};
54
2b893216 55Mech fgMechs[] =
56 {{ "HADF","total hadronic x-section according to FLUKA", "cm^{1}"},
57 { "INEF","hadronic inelastic x-section according to FLUKA", "cm^{1}"},
58 { "ELAF","hadronic elastic x-section according to FLUKA", "cm^{1}"},
59 { "HADG","total hadronic x-section according to GHEISHA", "cm^{1}"},
60 { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}"},
61 { "ELAG","hadronic elastic x-section according to GHEISHA", "cm^{1}"},
62 { "FISG","nuclear fission x-section according to GHEISHA", "cm^{1}"},
63 { "CAPG","neutron capture x-section according to GHEISHA", "cm^{1}"},
64 { "LOSS","stopping power", "cm^{1}"},
65 { "PHOT","photoelectric x-section", "MeV/cm"},
66 { "ANNI","positron annihilation x-section", "cm^{1}"},
67 { "COMP","Compton effect x-section", "cm^{1}"},
68 { "BREM","bremsstrahlung x-section", "cm^{1}"},
69 { "PAIR","photon and muon direct- pair x-section", "cm^{1}"},
70 { "DRAY","delta-rays x-section", "cm^{1}"},
71 { "PFIS","photo-fission x-section", "cm^{1}"},
72 { "RAYL","Rayleigh scattering x-section", "cm^{1}"},
73 { "MUNU","muon-nuclear interaction x-section", "cm^{1}"},
74 { "RANG","range", "cm"},
75 { "STEP","maximum step", "cm"},
76 { 0, 0, 0,}};
77
78//____________________________________________________________________
79/** @ingroup xsec_script
80 */
81struct MechValue
82{
83 MechValue(const Mech& m, size_t n)
84 : fMech(m), fValues(n), fStatus(0)
85 {}
86 bool Get(TGeant3* mc, std::vector<float>& t,
87 std::vector<float>& c, int medNo, int pidNo)
88 {
89 TGraph g(t.size());
90 g.SetName(fMech.fName);
91 g.SetTitle(fMech.fTitle);
92 g.GetXaxis()->SetTitle("p [GeV]");
93 g.GetYaxis()->SetTitle(fMech.fUnit);
94 int ixst;
95 mc->Gftmat(medNo, pidNo, fMech.fName, t.size(),
96 &(t[0]), &(fValues[0]), &(c[0]), ixst);
97 fStatus = ixst;
98 if (!fStatus) return false;
99 for (size_t i = 0; i < t.size(); i++) g.SetPoint(i, t[i], fValues[i]);
100 g.Write();
101 return true;
102 }
103 const Mech& fMech;
104 std::vector<float> fValues;
105 bool fStatus;
106};
107
108//____________________________________________________________________
109/** @ingroup xsec_script
110 */
111struct XSection
112{
113 XSection(size_t n=91, float emin=1e-5, float emax=1e4)
114 : fTKine(n),
115 fCuts(5)
116 {
117 float dp = 1. / log10(emax/emin);
118 float pmin = log10(emin);
119 fTKine[0] = emin;
120 for (size_t i = 1; i < fTKine.size(); i++) {
121 float el = pmin + i * dp;
122 fTKine[i] = pow(10, el);
123 }
124 for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i)
125 *i = 1e-4;
126 Mech* mech = &(fgMechs[0]);
127 size_t i = 0;
128 while (mech->fName) {
129 fMechs.push_back(new MechValue(*mech, n));
130 mech++;
131 }
132 }
133 void Run(const std::string& medName, const std::string& pdgName)
134 {
135 TGeant3* mc = static_cast<TGeant3*>(gMC);
136 if (!mc) {
137 std::cerr << "Couldn't get VMC" << std::endl;
138 return;
139 }
140 TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
141 if (!medium) {
142 std::cerr << "Couldn't find medium " << medName << std::endl;
143 return;
144 }
145 int medNo = medium->GetMaterial()->GetUniqueID();
146 TDatabasePDG* pdgDb = TDatabasePDG::Instance();
147 TParticlePDG* pdgP = pdgDb->GetParticle(pdgName.c_str());
148 if (!pdgP) {
149 std::cerr << "Couldn't find particle " << pdgName << std::endl;
150 return;
151 }
152 int pdgNo = pdgP->PdgCode();
153 int pidNo = mc->IdFromPDG(pdgNo);
154
155 std::stringstream vars;
156 vars << "T/F";
157
158 size_t nOk = 0;
159 // Loop over defined mechanisms
160 for (mech_array::iterator i = fMechs.begin(); i != fMechs.end(); ++i) {
161 if (!(*i)->Get(mc, fTKine, fCuts, medNo, pidNo))continue;
162 vars << ":" << (*i)->fMech.fName;
163 nOk ++;
164 }
165
166 std::stringstream tName;
167 tName << medName << "_" << pdgName;
168 TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
169
170 float_array cache(nOk+1);
171 tree->Branch("xsec", &(cache[0]), vars.str().c_str());
172 for (size_t i = 0; i < fTKine.size(); i++) {
173 cache[0] = fTKine[i];
174 int k = 0;
175 for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
176 if (!(*j)->fStatus) continue;
177 cache[++k] = (*j)->fValues[i];
178 }
179 tree->Fill();
180 }
181 tree->Write();
182 }
183protected:
184 typedef std::vector<float> float_array;
185 typedef std::vector<MechValue*> mech_array;
186 float_array fTKine;
187 float_array fCuts;
188 mech_array fMechs;
189};
a1f80595 190
2b893216 191#if 0
bf000c32 192//____________________________________________________________________
9b48326f 193/** @ingroup xsec_script
194 @param medName
195 @param pdgName
196 @param n
197 @param emin
198 @param emax
199*/
a1f80595 200void
201GetXsection(const char* medName, const char* pdgName,
202 Int_t n=91, Float_t emin=1e-5, Float_t emax=1e4)
203{
204 TArrayF tkine(n);
205 Float_t dp = 1/TMath::Log10(emax/emin);
206 Float_t pmin = TMath::Log10(emin);
207 tkine[0] = emin;
208 for (Int_t i=1; i < tkine.fN; i++) {
209 Float_t el = pmin + i * dp;
210 tkine[i] = TMath::Power(10, el);
211 }
212 TArrayF cuts(5);
213 cuts.Reset(1e-4);
214
215 Mech mechs[] =
216 {{ "HADF","total hadronic x-section according to FLUKA","cm^{1}",n,0},
217 { "INEF","hadronic inelastic x-section according to FLUKA","cm^{1}",n,0},
218 { "ELAF","hadronic elastic x-section according to FLUKA","cm^{1}",n,0},
219 { "HADG","total hadronic x-section according to GHEISHA","cm^{1}",n,0},
220 { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}",n,0},
221 { "ELAG","hadronic elastic x-section according to GHEISHA","cm^{1}",n,0},
222 { "FISG","nuclear fission x-section according to GHEISHA","cm^{1}",n,0},
223 { "CAPG","neutron capture x-section according to GHEISHA","cm^{1}",n,0},
224 { "LOSS","stopping power","cm^{1}",n,0},
225 { "PHOT","photoelectric x-section","MeV/cm",n,0},
226 { "ANNI","positron annihilation x-section","cm^{1}",n,0},
227 { "COMP","Compton effect x-section","cm^{1}",n,0},
228 { "BREM","bremsstrahlung x-section","cm^{1}",n,0},
229 { "PAIR","photon and muon direct- pair x-section","cm^{1}",n,0},
230 { "DRAY","delta-rays x-section","cm^{1}",n,0},
231 { "PFIS","photo-fission x-section","cm^{1}",n,0},
232 { "RAYL","Rayleigh scattering x-section","cm^{1}",n,0},
233 { "MUNU","muon-nuclear interaction x-section","cm^{1}",n,0},
234 { "RANG","range","cm",n,0},
235 { "STEP","maximum step","cm",n,0},
236 { 0, 0, 0, 0, 0}};
237 TGeant3* mc = (TGeant3*)gMC;
238 if (!mc) {
239 std::cerr << "Couldn't get VMC" << std::endl;
240 return;
241 }
242 TGeoMedium* medium = gGeoManager->GetMedium(medName);
243 if (!medium) {
244 std::cerr << "Couldn't find medium " << medName << std::endl;
245 return;
246 }
247 Int_t medNo = medium->GetMaterial()->GetUniqueID();
248 TDatabasePDG* pdgDb = TDatabasePDG::Instance();
249 TParticlePDG* pdgP = pdgDb->GetParticle(pdgName);
250 if (!pdgP) {
251 std::cerr << "Couldn't find particle " << pdgName << std::endl;
252 return;
253 }
254 Int_t pdgNo = pdgP->PdgCode();
255 Int_t pidNo = mc->IdFromPDG(pdgNo);
256
257 Mech* mech = &(mechs[0]);
258 Int_t nMech = 0;
259 Int_t nOk = 0;
260 TString vars("T/F");
261 while (mech->name) {
262 cout << mech->name << ": " << mech->title << " ... " << std::flush;
263 nMech++;
264 Int_t ixst;
265 mc->Gftmat(medNo, pidNo, mech->name, n,
266 tkine.fArray, mech->values.fArray, cuts.fArray, ixst);
267 mech->status = ixst;
268 if (ixst) {
269 nOk++;
270 vars.Append(Form(":%s", mech->name));
271 if (!strcmp("LOSS", mech->name)) {
272 for (Int_t i = 0; i < n; i++)
273 std::cout << i << "\t" << tkine[i] << "\t"
274 << mech->values[i] << std::endl;
275 }
276 }
277 std::cout << (ixst ? "ok" : "failed") << std::endl;
278 mech++;
279 }
280 // TFile* file = TFile::Open(Form("xsec-%d.root", pdgNo),
281 // "RECREATE");
282 TArrayF cache(nOk+1);
283 TTree* tree = new TTree(Form("%s_%s", medName, pdgName),
284 Form("%s_%s", medName, pdgName));
285 tree->Branch("xsec", cache.fArray, vars.Data());
286 for (Int_t i = 0; i < n; i++) {
287 cache[0] = tkine[i];
288 Int_t k = 0;
289 for (Int_t j = 0; j < nMech; j++) {
290 if (mechs[j].status) {
291 if (!strcmp(mechs[j].name, "LOSS"))
292 std::cout << tkine[i] << "\t" << mechs[j].values[i] << std::endl;
293 cache[k+1] = mechs[j].values[i];
294 k++;
295 }
296 }
297 std::cout << k << "\t" << (k == nOk) << std::endl;
298 tree->Fill();
299 }
300 tree->Write();
301}
2b893216 302#endif
303#endif
a1f80595 304//____________________________________________________________________
305//
306// EOF
307//