AliRoot Core  da88d91 (da88d91)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GetXsection.C
Go to the documentation of this file.
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 #ifdef __CINT__
16 XSection()
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
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>
35 #include <TGraph.h>
36 #include <TAxis.h>
37 #include <cmath>
38 #include <string>
39 #include <vector>
40 #include <sstream>
41 #include <iomanip>
45 //____________________________________________________________________
48 struct Mech
49 {
50  char* fName;
51  char* fTitle;
52  char* fUnit;
53 };
54 
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 //____________________________________________________________________
81 struct 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 //____________________________________________________________________
111 struct 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  }
183 protected:
184  typedef std::vector<float> float_array;
185  typedef std::vector<MechValue*> mech_array;
189 };
190 
191 #if 0
192 //____________________________________________________________________
200 void
201 GetXsection(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 }
302 #endif
303 #endif
304 //____________________________________________________________________
305 //
306 // EOF
307 //
TFile * Open(const char *filename, Long64_t &nevents)
bool Get(TGeant3 *mc, std::vector< float > &t, std::vector< float > &c, int medNo, int pidNo)
Definition: GetXsection.C:86
TROOT * gROOT
char * fUnit
Definition: GetXsection.C:52
TTree * tree
void XSection(const char *pdgName="pi-")
Definition: XSection.C:265
Double_t t
Definition: AliFMDv1.cxx:90
char * fTitle
Definition: GetXsection.C:51
Mech fgMechs[]
Definition: GetXsection.C:55
MechValue(const Mech &m, size_t n)
Definition: GetXsection.C:83
std::vector< float > fValues
Definition: GetXsection.C:104
mech_array fMechs
Definition: GetXsection.C:188
AliRun * gAlice
float_array fCuts
Definition: GetXsection.C:187
std::vector< float > float_array
Definition: GetXsection.C:184
void Run(const std::string &medName, const std::string &pdgName)
Definition: GetXsection.C:133
XSection(size_t n=91, float emin=1e-5, float emax=1e4)
Definition: GetXsection.C:113
const Mech & fMech
Definition: GetXsection.C:103
char * fName
Definition: GetXsection.C:50
std::vector< MechValue * > mech_array
Definition: GetXsection.C:185
float_array fTKine
Definition: GetXsection.C:186
bool fStatus
Definition: GetXsection.C:105