AliRoot Core  da88d91 (da88d91)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
XSection.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 #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>
41 //____________________________________________________________________
44 struct Mech
45 {
46  char* fName;
47  char* fTitle;
48  char* fUnit;
49  int fColor;
50  int fStyle;
51 };
52 
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 //____________________________________________________________________
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 //____________________________________________________________________
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;
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 //
void Draw()
Definition: XSection.C:213
TFile * Open(const char *filename, Long64_t &nevents)
TGraphErrors * fELoss
Definition: XSection.C:136
void Run(const std::string &medName, const std::string &pdgName)
Definition: XSection.C:163
TFile f("CalibObjects.root")
float_array fTKine
Definition: XSection.C:257
Mech fgMechs[]
Definition: XSection.C:53
mech_array fMechs
Definition: XSection.C:259
char * fUnit
Definition: GetXsection.C:52
TTree * tree
void XSection(const char *pdgName="pi-")
Definition: XSection.C:265
bool init
Definition: XSection.C:263
Double_t t
Definition: AliFMDv1.cxx:90
int fStyle
Definition: XSection.C:50
Bool_t bg
Definition: RunAnaESD.C:6
TParticlePDG * fPDG
Definition: XSection.C:260
std::vector< float > float_array
Definition: XSection.C:255
char * fTitle
Definition: GetXsection.C:51
MechValue(const Mech &m, size_t n)
Definition: XSection.C:81
std::vector< float > fValues
Definition: GetXsection.C:104
std::vector< MechValue * > mech_array
Definition: XSection.C:256
AliRun * gAlice
TGraph fGraph
Definition: XSection.C:135
int fColor
Definition: XSection.C:49
const Mech & fMech
Definition: GetXsection.C:103
TGraph & Draw()
Definition: XSection.C:126
float_array fCuts
Definition: XSection.C:258
XSections(size_t n=91, float emin=1e-5, float emax=1e4)
Definition: XSection.C:144
char * fName
Definition: GetXsection.C:50
bool Get(TGeant3 *mc, std::vector< float > &t, std::vector< float > &c, int medNo, int pidNo, float mass)
Definition: XSection.C:105
bool fStatus
Definition: GetXsection.C:105