Adding scripts directory including:
[u/mrichter/AliRoot.git] / test / vmctest / scripts / cuts.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 // $Id$
17
18 // Program to generate the plots with energy cuts and corresponding
19 // cuts in range for all materials per detector.
20 // The program uses a text file output from geant4_vmc
21 // (in development version only).
22 //
23 // By I.Hrivnacova, IPN Orsay
24
25 #include <Rtypes.h>
26 #include <TPaveText.h>
27 #include <TString.h>
28 #include <TGraph.h>
29 #include <TCanvas.h>
30 #include <TAxis.h>
31 #include <TLegend.h>
32
33 #include <map>
34 #include <vector>
35 #include <string>
36 #include <fstream>
37 #include <iostream>
38
39 using namespace std;
40
41 struct MaterialData {
42   MaterialData(string name, 
43                Double_t x1, Double_t x2, Double_t x3, 
44                Double_t x4, Double_t x5, Double_t x6)
45     : matName(name), 
46       rangeGam(x1), 
47       rangeEle(x2), 
48       cutGam(x3), 
49       cutEle(x4),
50       cutGamLimit(x5), 
51       cutEleLimit(x6) {}
52   string matName;
53   Double_t rangeGam;
54   Double_t rangeEle;
55   Double_t cutGam;
56   Double_t cutEle;
57   Double_t cutGamLimit;
58   Double_t cutEleLimit;
59 };  
60   
61 map<string, vector<MaterialData> > materialDataMap;  
62
63 void readCuts() 
64 {
65   std::ifstream input("cuts.txt");
66   while ( ! input.eof() ) {
67     string matName;
68     Int_t number;
69     Double_t x1, x2, x3, x4, x5, x6;
70     input >> number;
71     if ( input.eof() ) break;
72     input  >> matName >> x1 >> x2 >> x3 >> x4 >> x5 >> x6;
73     // cout << ".. reading " << number << " " << matName << endl;
74     
75     if ( x1 > 1e10 ) x1 = 1e04;
76     if ( x2 > 1e10 ) x2 = 1e04;
77     if ( x3 > 1e10 ) x3 = 1e04;
78     if ( x4 > 1e10 ) x4 = 1e04;
79     if ( x5 > 1e10 ) x5 = 1e04;
80     if ( x6 > 1e10 ) x6 = 1e04;
81     
82     string detName = matName;
83     detName.erase(detName.find('_'), detName.size()-detName.find('_'));
84
85     std::map<string, std::vector<MaterialData> >::iterator it
86       = materialDataMap.find(detName);
87     if ( it == materialDataMap.end() ) {
88       materialDataMap.insert(
89         std::pair<string, vector<MaterialData> >(detName, vector<MaterialData>()));
90     }
91     materialDataMap[detName].
92       push_back(MaterialData(matName, x1, x2, x3, x4, x5, x6));
93   }
94 }        
95   
96 void plotCuts(const string& detName) 
97 {
98   cout << "Processing " << detName << " ..." << endl;
99   
100   //if ( detName == "DIPO" ) return;
101   //if ( detName == "HALL" ) return;
102   //if ( detName == "SHIL" ) return;
103   //if ( detName == "ZDC" ) return;
104
105   TPaveText* paveText = new TPaveText(0.1, 0.1, 0.98, 0.98);
106   paveText->SetTextAlign(11);
107
108   std::vector<MaterialData> materialData 
109     = materialDataMap[detName];
110
111   Double_t matNumber[300];
112   Double_t rangeGam[300];
113   Double_t rangeEle[300];
114   Double_t cutGam[300];
115   Double_t cutEle[300];
116   Double_t cutGamLimit[300];
117   Double_t cutEleLimit[300];
118   std::vector<string> matNames;
119  
120   for (UInt_t i=0; i<materialData.size(); i++ ) {
121   
122     matNumber[i] = i;
123     rangeGam[i] = materialData[i].rangeGam;  
124     rangeEle[i] = materialData[i].rangeEle;
125     cutGam[i] = materialData[i].cutGam;
126     cutEle[i] = materialData[i].cutEle;
127     cutGamLimit[i] = materialData[i].cutGamLimit;
128     cutEleLimit[i] = materialData[i].cutEleLimit;
129    
130     TString legend;
131     legend += i;
132     legend += "  ";
133     legend += materialData[i].matName.data();    
134     paveText->AddText(legend.Data());
135   }
136   
137   TGraph* gr1 = new TGraph(materialData.size(), matNumber, rangeGam); 
138   TGraph* gr2 = new TGraph(materialData.size(), matNumber, cutGam); 
139   TGraph* gr3 = new TGraph(materialData.size(), matNumber, cutGamLimit); 
140   TGraph* gr4 = new TGraph(materialData.size(), matNumber, rangeEle); 
141   TGraph* gr5 = new TGraph(materialData.size(), matNumber, cutEle); 
142   TGraph* gr6 = new TGraph(materialData.size(), matNumber, cutEleLimit); 
143
144   // gamma range cut
145   gr1->SetMarkerColor(kBlue);
146   gr1->SetMarkerStyle(22);
147   gr1->SetTitle("Range cut for gamma");
148   gr1->GetXaxis()->SetTitle("Material number");
149   gr1->GetYaxis()->SetTitle("Range [mm]");
150   gr1->SetMinimum(1e-03);
151   gr1->SetMaximum(1e+04);
152
153   // gamma energy threshold
154   gr2->SetMarkerColor(kBlue);
155   gr2->SetMarkerStyle(22);
156   gr2->SetTitle("Energy threshold for gamma");
157   gr2->GetXaxis()->SetTitle("Material number");
158   gr2->GetYaxis()->SetTitle("Energy [MeV]");
159   gr2->SetMinimum(1e-04);
160   gr2->SetMaximum(1e+04);
161
162   // gamma user limit
163   gr3->SetMarkerColor(kBlack);
164   gr3->SetMarkerStyle(23);
165   gr3->SetTitle("User limit for gamma (GUTGAM)");
166   gr3->GetXaxis()->SetTitle("Material number");
167   gr3->GetYaxis()->SetTitle("Energy [MeV]");
168   gr3->SetMinimum(1e-04);
169   gr3->SetMaximum(1e+04);
170
171   // e- range cut
172   gr4->SetMarkerColor(kRed);
173   gr4->SetMarkerStyle(22);
174   gr4->SetTitle("Range cut for e-");
175   gr4->GetXaxis()->SetTitle("Material number");
176   gr4->GetYaxis()->SetTitle("Range [mm]");
177   gr4->GetYaxis()->SetRange(1e-03, 1e+04);
178   gr4->SetMinimum(1e-03);
179   gr4->SetMaximum(1e+04);
180
181   // e- energy threshold
182   gr5->SetMarkerColor(kRed);
183   gr5->SetMarkerStyle(22);
184   gr5->SetTitle("Energy threshold for e-");
185   gr5->GetXaxis()->SetTitle("Material number");
186   gr5->GetYaxis()->SetTitle("Energy [MeV]");
187   gr5->SetMinimum(1e-04);
188   gr5->SetMaximum(1e+04);
189
190   // e- user limit
191   gr6->SetMarkerColor(kBlack);
192   gr6->SetMarkerStyle(23);
193   gr6->SetTitle("User limit for e- (CUTELE)");
194   gr6->GetXaxis()->SetTitle("Material number");
195   gr6->GetYaxis()->SetTitle("Energy [MeV]");
196   gr6->SetMinimum(1e-04);
197   gr6->SetMaximum(1e+04);
198
199   //TLegend* leg = new TLegend(0.1,0.7,0.48,0.9);
200   TLegend* leg = new TLegend(0.1, 0.1, 0.98, 0.98);
201   leg->SetHeader("Energy threshold/User cuts");
202   leg->AddEntry(gr2,"Energy threshold for gamma","P");
203   leg->AddEntry(gr3,"User cut for gamma (CUGAM)","P");
204   leg->AddEntry(gr5,"Energy threshold for e-","P");
205   leg->AddEntry(gr6,"User cut for e- (CUTELE)","P");
206
207   TCanvas* canvas = new TCanvas("Range cuts", "Range cuts", 1200, 800);
208   canvas->Divide(3,2);
209   
210   // Draw graph using the function Draw()
211   canvas->cd(1);
212   canvas->cd(1)->SetLogy();
213   gr1->Draw("AP");
214  
215   canvas->cd(2);
216   canvas->cd(2)->SetLogy();
217   gr2->Draw("AP");
218
219   //canvas->cd(3);
220   //canvas->cd(3)->SetLogy();
221   gr3->Draw("P");
222  
223   canvas->cd(3);
224   leg->Draw();
225
226   canvas->cd(4);
227   canvas->cd(4)->SetLogy();
228   gr4->Draw("AP");
229
230   canvas->cd(5);
231   canvas->cd(5)->SetLogy();
232   gr5->Draw("AP");
233  
234   //canvas->cd(6);
235   //canvas->cd(6)->SetLogy();
236   gr6->Draw("P");
237
238   canvas->cd(6);
239
240   //TCanvas* canvas2 = new TCanvas("Materials", "Materials", 200, 800);
241   //canvas2->cd();
242   paveText->Draw();
243   
244   string outName = detName;
245   outName += "_cuts.gif";
246   canvas->SaveAs(outName.data());
247   
248   //string outName2 = detName;
249   //outName2 += "_cuts_legend.gif";
250   //canvas2->SaveAs(outName2.data());
251 }
252
253 int main() {
254   
255   readCuts();
256   
257   std::map<string, std::vector<MaterialData> >::iterator it;
258   for ( it= materialDataMap.begin(); it != materialDataMap.end(); it++ )
259     plotCuts(it->first);
260     
261   return 0;  
262 }  
263
264
265