Script to get # of dead channels from OCDB
[u/mrichter/AliRoot.git] / FMD / scripts / XSection.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 #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>
38 /** @defgroup xsec_script X-section script 
39     @ingroup FMD_script 
40 */
41 //____________________________________________________________________
42 /** @ingroup xsec_script
43  */
44 struct Mech 
45 {
46   char* fName;
47   char* fTitle;
48   char* fUnit;
49   int   fColor;
50   int   fStyle;
51 };
52
53 Mech fgMechs[] = 
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 //____________________________________________________________________
77 /** @ingroup xsec_script
78  */
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 //____________________________________________________________________
140 /** @ingroup xsec_script
141  */
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;
257   float_array   fTKine;
258   float_array   fCuts;
259   mech_array    fMechs;
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 //