]>
Commit | Line | Data |
---|---|---|
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 | // | |
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> | |
9b48326f | 25 | /** @defgroup xsec_script X-section script |
26 | @ingroup FMD_script | |
27 | */ | |
bf000c32 | 28 | //____________________________________________________________________ |
9b48326f | 29 | /** @ingroup xsec_script |
30 | */ | |
a1f80595 | 31 | struct Mech |
32 | { | |
33 | char* name; | |
34 | char* title; | |
35 | char* unit; | |
36 | TArrayF values; | |
37 | int status; | |
38 | }; | |
39 | ||
40 | ||
bf000c32 | 41 | //____________________________________________________________________ |
9b48326f | 42 | /** @ingroup xsec_script |
43 | @param medName | |
44 | @param pdgName | |
45 | @param n | |
46 | @param emin | |
47 | @param emax | |
48 | */ | |
a1f80595 | 49 | void |
50 | GetXsection(const char* medName, const char* pdgName, | |
51 | Int_t n=91, Float_t emin=1e-5, Float_t emax=1e4) | |
52 | { | |
53 | TArrayF tkine(n); | |
54 | Float_t dp = 1/TMath::Log10(emax/emin); | |
55 | Float_t pmin = TMath::Log10(emin); | |
56 | tkine[0] = emin; | |
57 | for (Int_t i=1; i < tkine.fN; i++) { | |
58 | Float_t el = pmin + i * dp; | |
59 | tkine[i] = TMath::Power(10, el); | |
60 | } | |
61 | TArrayF cuts(5); | |
62 | cuts.Reset(1e-4); | |
63 | ||
64 | Mech mechs[] = | |
65 | {{ "HADF","total hadronic x-section according to FLUKA","cm^{1}",n,0}, | |
66 | { "INEF","hadronic inelastic x-section according to FLUKA","cm^{1}",n,0}, | |
67 | { "ELAF","hadronic elastic x-section according to FLUKA","cm^{1}",n,0}, | |
68 | { "HADG","total hadronic x-section according to GHEISHA","cm^{1}",n,0}, | |
69 | { "INEG","hadronic inelastic x-section according to GHEISHA","cm^{1}",n,0}, | |
70 | { "ELAG","hadronic elastic x-section according to GHEISHA","cm^{1}",n,0}, | |
71 | { "FISG","nuclear fission x-section according to GHEISHA","cm^{1}",n,0}, | |
72 | { "CAPG","neutron capture x-section according to GHEISHA","cm^{1}",n,0}, | |
73 | { "LOSS","stopping power","cm^{1}",n,0}, | |
74 | { "PHOT","photoelectric x-section","MeV/cm",n,0}, | |
75 | { "ANNI","positron annihilation x-section","cm^{1}",n,0}, | |
76 | { "COMP","Compton effect x-section","cm^{1}",n,0}, | |
77 | { "BREM","bremsstrahlung x-section","cm^{1}",n,0}, | |
78 | { "PAIR","photon and muon direct- pair x-section","cm^{1}",n,0}, | |
79 | { "DRAY","delta-rays x-section","cm^{1}",n,0}, | |
80 | { "PFIS","photo-fission x-section","cm^{1}",n,0}, | |
81 | { "RAYL","Rayleigh scattering x-section","cm^{1}",n,0}, | |
82 | { "MUNU","muon-nuclear interaction x-section","cm^{1}",n,0}, | |
83 | { "RANG","range","cm",n,0}, | |
84 | { "STEP","maximum step","cm",n,0}, | |
85 | { 0, 0, 0, 0, 0}}; | |
86 | TGeant3* mc = (TGeant3*)gMC; | |
87 | if (!mc) { | |
88 | std::cerr << "Couldn't get VMC" << std::endl; | |
89 | return; | |
90 | } | |
91 | TGeoMedium* medium = gGeoManager->GetMedium(medName); | |
92 | if (!medium) { | |
93 | std::cerr << "Couldn't find medium " << medName << std::endl; | |
94 | return; | |
95 | } | |
96 | Int_t medNo = medium->GetMaterial()->GetUniqueID(); | |
97 | TDatabasePDG* pdgDb = TDatabasePDG::Instance(); | |
98 | TParticlePDG* pdgP = pdgDb->GetParticle(pdgName); | |
99 | if (!pdgP) { | |
100 | std::cerr << "Couldn't find particle " << pdgName << std::endl; | |
101 | return; | |
102 | } | |
103 | Int_t pdgNo = pdgP->PdgCode(); | |
104 | Int_t pidNo = mc->IdFromPDG(pdgNo); | |
105 | ||
106 | Mech* mech = &(mechs[0]); | |
107 | Int_t nMech = 0; | |
108 | Int_t nOk = 0; | |
109 | TString vars("T/F"); | |
110 | while (mech->name) { | |
111 | cout << mech->name << ": " << mech->title << " ... " << std::flush; | |
112 | nMech++; | |
113 | Int_t ixst; | |
114 | mc->Gftmat(medNo, pidNo, mech->name, n, | |
115 | tkine.fArray, mech->values.fArray, cuts.fArray, ixst); | |
116 | mech->status = ixst; | |
117 | if (ixst) { | |
118 | nOk++; | |
119 | vars.Append(Form(":%s", mech->name)); | |
120 | if (!strcmp("LOSS", mech->name)) { | |
121 | for (Int_t i = 0; i < n; i++) | |
122 | std::cout << i << "\t" << tkine[i] << "\t" | |
123 | << mech->values[i] << std::endl; | |
124 | } | |
125 | } | |
126 | std::cout << (ixst ? "ok" : "failed") << std::endl; | |
127 | mech++; | |
128 | } | |
129 | // TFile* file = TFile::Open(Form("xsec-%d.root", pdgNo), | |
130 | // "RECREATE"); | |
131 | TArrayF cache(nOk+1); | |
132 | TTree* tree = new TTree(Form("%s_%s", medName, pdgName), | |
133 | Form("%s_%s", medName, pdgName)); | |
134 | tree->Branch("xsec", cache.fArray, vars.Data()); | |
135 | for (Int_t i = 0; i < n; i++) { | |
136 | cache[0] = tkine[i]; | |
137 | Int_t k = 0; | |
138 | for (Int_t j = 0; j < nMech; j++) { | |
139 | if (mechs[j].status) { | |
140 | if (!strcmp(mechs[j].name, "LOSS")) | |
141 | std::cout << tkine[i] << "\t" << mechs[j].values[i] << std::endl; | |
142 | cache[k+1] = mechs[j].values[i]; | |
143 | k++; | |
144 | } | |
145 | } | |
146 | std::cout << k << "\t" << (k == nOk) << std::endl; | |
147 | tree->Fill(); | |
148 | } | |
149 | tree->Write(); | |
150 | } | |
151 | //____________________________________________________________________ | |
152 | // | |
153 | // EOF | |
154 | // |