source: branches/Corsika7500Compatibility/macros/fact/ratescan.C@ 18590

Last change on this file since 18590 was 13248, checked in by tbretz, 13 years ago
Plot the results of a ratescan from a slow control file.
File size: 7.9 KB
Line 
1MStatusDisplay *d = new MStatusDisplay;
2TList list[201];
3
4void Fit(TGraph &g, double rate_start=3e6, double rate_kink=100, double rate_end=10)
5{
6
7 double begin = -1;
8 double end = -1;
9
10 double kinkx = -1;
11 double kinky = 0;
12
13 TGraph gl;
14 for (int j=0; j<g.GetN(); j++)
15 {
16 gl.SetPoint(gl.GetN(), g.GetX()[j], log10(g.GetY()[j]));
17 if (g.GetY()[j]<rate_start && begin<0)
18 begin = g.GetX()[j];
19 if (g.GetY()[j]<rate_end && end<0)
20 end = g.GetX()[j];
21 if (g.GetY()[j]<rate_kink && kinkx<0)
22 {
23 kinkx = g.GetX()[j];
24 kinky = g.GetY()[j];
25 }
26 }
27
28 TF1 func2("f2", "log10("
29 " [1]*exp([2]*(x-[0]))"
30 "+ [1]*exp([3]*(x-[0]))"
31 //"+ [4]*exp([5]*(x-700))"
32 ")", begin, end);
33
34 Double_t v00 = g.Eval(kinkx-150, 0, "S");
35 Double_t v01 = g.Eval(kinkx- 50, 0, "S");
36
37 Double_t v10 = g.Eval(kinkx+100, 0, "S");
38 Double_t v11 = g.Eval(kinkx+200, 0, "S");
39
40 //Double_t v20 = g.Eval(700-50, 0, "S");
41 //Double_t v21 = g.Eval(700, 0, "S");
42 //Double_t v22 = g.Eval(700+50, 0, "S");
43
44 func2.SetParameter(0, kinkx);
45 func2.SetParameter(1, kinky);
46 func2.SetParameter(2, log(v01/v00)/100);
47 func2.SetParameter(3, log(v11/v10)/100);
48
49 func2.FixParameter(4, 0);
50 func2.FixParameter(5, 0);
51 func2.SetLineWidth(1);
52 func2.SetLineColor(kBlue);
53
54 gl.Fit(&func2, "N0", "", begin, end);
55
56 TF1 func("f1",
57 "[0]*exp([1]*(x-[6])) + "
58 "[2]*exp([3]*(x-[6])) + "
59 "[4]*exp([5]*(x-700))", begin, end);
60 func.FixParameter(0, func2.GetParameter(1));
61 func.FixParameter(1, func2.GetParameter(2));
62 func.FixParameter(2, func2.GetParameter(1));
63 func.FixParameter(3, func2.GetParameter(3));
64 func.FixParameter(4, func2.GetParameter(4));
65 func.FixParameter(5, func2.GetParameter(5));
66 func.FixParameter(6, func2.GetParameter(0));
67 func.SetLineWidth(2);
68 func.SetLineColor(kBlue);
69 func.DrawClone("same");
70 func.SetLineWidth(1);
71
72 func.SetLineStyle(kDashed);
73 func.SetRange(0, 1000);
74
75 func.FixParameter(0, 0);
76 func.FixParameter(2, 0);
77 func.FixParameter(4, func2.GetParameter(4));
78 func.DrawClone("same");
79
80 func.FixParameter(0, func2.GetParameter(1));
81 func.FixParameter(2, 0);
82 func.FixParameter(4, 0);
83 func.DrawClone("same");
84
85 func.FixParameter(0, 0);
86 func.FixParameter(2, func2.GetParameter(1));
87 func.FixParameter(4, 0);
88 func.DrawClone("same");
89
90
91 TLatex txt;
92 txt.SetTextColor(kBlue);
93 txt.SetTextAlign(33);
94 txt.DrawLatex(func2.GetParameter(0), func2.GetParameter(1),
95 Form("%.1fHz @ %.0f", func2.GetParameter(1), func2.GetParameter(0)));
96
97 // [1]*exp([2]*(x-[0])) / exp(1) = [1]*exp([3]*(x-[0]))
98 // exp([2]*(x-[0])) / exp(1) = exp([3]*(x-[0]))
99 // exp([2]*(x-[0])) / exp([3]*(x-[0])) = exp(1)
100 // exp([2]*(x-[0]) - [3]*(x-[0])) = exp(1)
101 // [2]*(x-[0]) - [3]*(x-[0]) = 1
102 //
103 // x = 1/([2]-[3]) + [0]
104
105
106 Double_t th = 1. / (func2.GetParameter(3)-func2.GetParameter(2)) + func2.GetParameter(0);
107 Double_t rate = pow(10, func2.Eval(th));
108
109 txt.SetTextAlign(11);
110 txt.DrawLatex(th+10, rate, Form("%.1fHz @ %.0f", rate, th));
111
112 txt.SetTextAlign(12);
113 txt.DrawLatex(180, 3e6, Form("%.1f", 1./func2.GetParameter(2)));
114 txt.DrawLatex(350, 3, Form("%.0f", 1./func2.GetParameter(3)));
115
116 TLine line;
117 line.SetLineStyle(kDashed);
118 line.SetLineColor(kBlue);
119 line.DrawLine(th, 1e-2, th, 1e8);
120
121}
122
123void plot(const char *tab, const char *title="")
124{
125 gROOT->SetSelectedPad(0);
126 d->AddTab(tab);
127
128 gPad->SetGrid();
129 gPad->SetLogy();
130
131 TH1D h0("frame", title, 1000, 0, 1000);
132 h0.SetStats(kFALSE);
133 h0.SetXTitle("Threshold [dac]");
134 h0.SetYTitle("Trigger rate [Hz]");
135 h0.SetDirectory(0);
136 h0.SetMaximum(1e8);
137 h0.SetMinimum(0.01);
138 h0.DrawCopy();
139}
140
141
142void same(const char *id, const char *fmt, bool fit=true, bool boards=false, bool patches=false)
143{
144 for (int i=0; i<160; i++)
145 {
146 TGraph *g = (TGraph*)list[i+41].FindObject(id);
147 if (!g)
148 {
149 cout << "==W==> " << id << " not found in i=" << i+41 << endl;
150 return;
151 }
152
153 g->SetMarkerStyle(kFullDotSmall);
154 g->SetMarkerColor(kGray+1);
155 g->SetLineColor(kGray+1);
156 if (patches)
157 g->DrawClone(fmt);
158
159 if (g->Eval(400)>10)
160 cout << "Patch " << setw(3) << i << ": C" << i/40 << " B" << (i%40)/4 << " P" << (i%40)%4 << ": " << g->Eval(400) << "Hz @ th=400" << endl;
161 }
162
163 for (int i=0; i<40; i++)
164 {
165 TGraph *g = (TGraph*)list[i+1].FindObject(id);
166 if (!g)
167 {
168 cout << "==W==> " << id << " not found in i=" << i+1 << endl;
169 return;
170 }
171
172 g->SetMarkerStyle(kFullDotSmall);
173 g->SetMarkerColor(kGray+2);
174 g->SetLineColor(kGray+2);
175 if (boards)
176 g->DrawClone(fmt);
177
178 if (g->Eval(400)>10)
179 cout << "Board " << setw(3) << i << ": C" << i/10 << " B" << i%10 << ": " << g->Eval(400) << "Hz @ th=400" << endl;
180 }
181
182
183 TGraph *g001 = (TGraph*)list[0].FindObject(id);
184 if (!g001)
185 {
186 cout << "==W==> " << id << " not found in i=0" << endl;
187 return;
188 }
189
190 g001->SetMarkerStyle(boards||patches?kFullDotMedium:kFullDotSmall);
191 g001->DrawClone(fmt);
192
193 if (fit)
194 Fit(*g001);
195}
196
197void read(const char *filename)
198{
199 fits fin(filename);
200 if (!fin)
201 return;
202
203 cout << '\n' << filename << '\n';
204 cout << setfill('-') <<setw(strlen(filename)) << '-' << setfill(' ') << endl;
205
206 //fin.PrintColumns();
207
208 ULong64_t start;
209 UInt_t threshold;
210 Float_t time;
211 Float_t ontime;
212 Float_t Rc;
213 Float_t Rb[40];
214 Float_t Rp[160];
215
216 fin.SetPtrAddress("Data0", &start);
217 fin.SetPtrAddress("Data1", &threshold);
218 fin.SetPtrAddress("Data2", &time);
219 fin.SetPtrAddress("Data3", &ontime);
220 fin.SetPtrAddress("Data4", &Rc);
221 fin.SetPtrAddress("Data5", Rb);
222 fin.SetPtrAddress("Data6", Rp);
223
224 ULong_t first = 0;
225
226 TGraph *g[201];
227
228 while (fin.GetNextRow())
229 {
230 if (first==0 || first!=start)
231 {
232 first = start;
233
234 MTime t;
235 t.SetUnixTime(start);
236
237 TString name = t.GetString();
238 name = name(0, 19);
239
240 cout << name << endl;
241
242 for (int i=0; i<201; i++)
243 {
244 g[i] = new TGraph;
245 g[i]->SetMarkerStyle(kFullDotMedium);
246 g[i]->SetName(name);
247 g[i]->SetTitle(t.GetStringFmt("%H:%M:%S"));
248
249 list[i].Add(g[i]);
250 }
251 }
252
253 g[0]->SetPoint(g[0]->GetN(), threshold, Rc/ontime);
254 for (int i=0; i<40; i++)
255 g[i+1]->SetPoint(g[i+1]->GetN(), threshold, Rb[i]);
256 for (int i=0; i<160; i++)
257 g[i+41]->SetPoint(g[i+41]->GetN(), threshold, Rp[i]);
258 }
259}
260
261void ratescan()
262{
263 cout << setprecision(4);
264
265 // Read three files
266 read("/loc_data/aux/2012/02/18/20120218.RATE_SCAN_DATA.fits");
267 read("/loc_data/aux/2012/02/19/20120219.RATE_SCAN_DATA.fits");
268 read("/loc_data/aux/2012/02/21/20120221.RATE_SCAN_DATA.fits");
269
270 // add a new tab
271 plot("0.0V", "Temp control 19.02./22.02. (+0V)");
272
273 // Plot a curve, do not fit
274 same("19.02.2012 20:21:26", "LP", false);
275
276 // plot a curve, fit and plot boards
277 same("22.02.2012 04:19:15", "LP", true, true);
278
279 // raed another file
280 read("/loc_data/aux/2012/03/25/20120325.RATE_SCAN_DATA.fits");
281
282 // add another tab
283 plot("25/03", "Current control 25.03. (+0V)");
284
285 // plot curve, fit and plot boards and patches
286 same("26.03.2012 03:01:14", "LP", true, true, true);
287
288 // plot curve, do not fit, do not draw boards and patches
289 same("26.03.2012 03:21:47", "LP", false);
290
291 // write result to output file
292 d->SaveAs("/home/tbretz/ratescan.root");
293}
Note: See TracBrowser for help on using the repository browser.