Bug fix for reading raw data - thanks Frederic YERMIA, Diego, and Gines MARTINEZ
[u/mrichter/AliRoot.git] / FMD / scripts / GetXsection.C
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>
42 /** @defgroup xsec_script X-section script 
43     @ingroup FMD_script 
44 */
45 //____________________________________________________________________
46 /** @ingroup xsec_script
47  */
48 struct Mech 
49 {
50   char* fName;
51   char* fTitle;
52   char* fUnit;
53 };
54
55 Mech 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  */
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 //____________________________________________________________________
109 /** @ingroup xsec_script
110  */
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;
186   float_array  fTKine;
187   float_array  fCuts;
188   mech_array   fMechs;
189 };
190
191 #if 0
192 //____________________________________________________________________
193 /** @ingroup xsec_script
194     @param medName 
195     @param pdgName 
196     @param n 
197     @param emin 
198     @param emax 
199 */
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 //