1 //____________________________________________________________________
5 // Script to get the various cross sections, energy loss, and ranges
6 // of a particular particle type in a particular medium.
8 // This script should be compiled to speed it up.
10 // It creates a tree on the current output file, with the relevant
13 // Note, that VMC _must_ be the TGeant3TGeo VMC.
19 #include <TVirtualMC.h>
20 #include <TDatabasePDG.h>
22 #include <TGeoManager.h>
23 #include <TGeoMedium.h>
26 #include <TGraphErrors.h>
38 /** @defgroup xsec_script X-section script
41 //____________________________________________________________________
42 /** @ingroup xsec_script
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},
76 //____________________________________________________________________
77 /** @ingroup xsec_script
81 MechValue(const Mech& m, size_t n)
82 : fMech(m), fValues(n), fStatus(0), fELoss(0)
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;
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);
105 bool Get(TGeant3* mc, std::vector<float>& t,
106 std::vector<float>& c, int medNo, int pidNo,
110 mc->Gftmat(medNo, pidNo, fMech.fName, t.size(),
111 &(t[0]), &(fValues[0]), &(c[0]), 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]);
120 fELoss->SetPointError(i, 0, .5 * fValues[i]);
123 if (fELoss) fELoss->Write();
128 if (fELoss) fELoss->DrawClone("4 same");
129 fGraph.DrawClone("l same");
133 std::vector<float> fValues;
136 TGraphErrors* fELoss;
139 //____________________________________________________________________
140 /** @ingroup xsec_script
144 XSections(size_t n=91, float emin=1e-5, float emax=1e4)
148 float dp = 1. / log10(emax/emin);
149 float pmin = log10(emin);
151 for (size_t i = 1; i < fTKine.size(); i++) {
152 float el = pmin + i * dp;
153 fTKine[i] = pow(10, el);
155 for (float_array::iterator i = fCuts.begin(); i != fCuts.end(); ++i)
157 Mech* mech = &(fgMechs[0]);
158 while (mech->fName) {
159 fMechs.push_back(new MechValue(*mech, n));
163 void Run(const std::string& medName, const std::string& pdgName)
165 TGeant3* mc = static_cast<TGeant3*>(gMC);
167 std::cerr << "Couldn't get VMC" << std::endl;
170 TGeoMedium* medium = gGeoManager->GetMedium(medName.c_str());
172 std::cerr << "Couldn't find medium " << medName << std::endl;
175 int medNo = medium->GetMaterial()->GetUniqueID();
176 TDatabasePDG* pdgDb = TDatabasePDG::Instance();
177 fPDG = pdgDb->GetParticle(pdgName.c_str());
179 std::cerr << "Couldn't find particle " << pdgName << std::endl;
182 int pdgNo = fPDG->PdgCode();
183 int pidNo = mc->IdFromPDG(pdgNo);
185 std::stringstream vars;
186 vars << "betagamma/F";
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;
196 std::stringstream tName;
197 tName << medName << "_" << pdgName;
198 TTree* tree = new TTree(tName.str().c_str(), tName.str().c_str());
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();
205 for (mech_array::iterator j = fMechs.begin(); j != fMechs.end(); ++j) {
206 if (!(*j)->fStatus) continue;
207 cache[++k] = (*j)->fValues[i];
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]);
229 TCanvas* c = new TCanvas("c", "C", 700, 700);
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");
244 TLegend* l = new TLegend(0.45, 0.125, 0.90, 0.45);
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");
255 typedef std::vector<float> float_array;
256 typedef std::vector<MechValue*> mech_array;
265 void XSection(const char* pdgName="pi-")
268 gAlice->InitMC("$(ALICE_ROOT)/FMD/Config.C");
271 TFile* file = TFile::Open("xsec.root", "RECREATE");
273 xs.Run("FMD_Si$", pdgName);
277 //____________________________________________________________________