source: branches/MarsMoreSimulationTruth/datacenter/macros/fillsignal.C@ 18656

Last change on this file since 18656 was 9062, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 14.8 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz, 04/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Daniela Dorner, 04/2005 <mailto:dorner@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2008
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// fillsignal.C
29// ============
30//
31// This macro is used to read the calibration-/callisto-output files
32// signal00000.root.
33//
34// From this file the mean pedrms, the mean signal and the pulse position
35// for the inner and outer camera is extracted and inserted into the database
36// in the table Calibration, where the results of callisto are stored.
37// The sequence number is extracted from the filename.
38//
39// Usage:
40// .x fillsignal.C("/magic/data/callisto/0004/00047421/signal00047421.root", kTRUE)
41//
42// The second argument is the 'dummy-mode'. If it is kTRUE dummy-mode is
43// switched on and nothing will be written into the database. This is usefull
44// for tests.
45//
46// The macro can also be run without ACLiC but this is a lot slower...
47//
48// Remark: Running it from the commandline looks like this:
49// root -q -l -b fillsignal.C+\(\"filename\"\,kFALSE\) 2>&1 | tee fillsignal.log
50//
51// Make sure, that database and password are corretly set in a resource
52// file called sql.rc and the resource file is found.
53//
54// Returns 2 in case of failure, 1 in case of success and 0 if the connection
55// to the database is not working.
56//
57/////////////////////////////////////////////////////////////////////////////
58#include <iostream>
59#include <iomanip>
60
61#include <TRegexp.h>
62
63#include <TH2.h>
64#include <TFile.h>
65#include <TGraph.h>
66#include <TSQLResult.h>
67
68#include "MSQLMagic.h"
69
70#include "MStatusArray.h"
71#include "MSequence.h"
72#include "MHCamera.h"
73#include "MHVsTime.h"
74
75#include "MCalibrationPulseTimeCam.h"
76#include "MCalibrationPix.h"
77
78using namespace std;
79
80int Process(MSQLMagic &serv, TString fname)
81{
82 TFile file(fname, "READ");
83 if (!file.IsOpen())
84 {
85 cout << "ERROR - Could not find file " << fname << endl;
86 return 2;
87 }
88
89 TString meanextpulpos("NULL");
90 TString rmsextpulpos("NULL");
91
92 MCalibrationPulseTimeCam *pt;
93 file.GetObject("MCalibrationPulseTimeCam", pt);
94 if (pt)
95 {
96 Double_t meanextpul = pt->GetAverageArea(0).GetHiGainMean();
97 Double_t rmsextpul = pt->GetAverageArea(0).GetHiGainRms();
98
99 if (meanextpul>=0 || rmsextpulpos>=0)
100 {
101 meanextpulpos.Form("%6.2f", meanextpul);
102 rmsextpulpos.Form("%6.2f", rmsextpul);
103 }
104 }
105
106
107 MStatusArray arr;
108 if (arr.Read()<=0)
109 {
110 cout << "ERROR - Reading of MStatusDisplay failed." << endl;
111 return 2;
112 }
113
114 MHCamera *cam = (MHCamera*)arr.FindObjectInCanvas("PedRMS;avg", "MHCamera", "PedRMS");
115 if (!cam)
116 {
117 cout << "ERROR - Reading of PedRMS;avg failed." << endl;
118 return 2;
119 }
120
121 MHCamera *cal = (MHCamera*)arr.FindObjectInCanvas("CalPos;avg", "MHCamera", "CalPos");
122 MHCamera *pul = (MHCamera*)arr.FindObjectInCanvas("PulsePos;avg", "MHCamera", "PulsePos");
123 if (!pul)
124 {
125 cout << "ERROR - Reading of PulsePos;avg failed." << endl;
126 return 2;
127 }
128/*
129 MHCamera *difflo = (MHCamera*)arr.FindObjectInCanvas("DiffLo;avg", "MHCamera", "DiffLo");
130 if (!difflo)
131 {
132 cout << "ERROR - Reading of DiffLo;avg failed." << endl;
133 return 2;
134 }
135 MHCamera *diffhi = (MHCamera*)arr.FindObjectInCanvas("DiffHi;avg", "MHCamera", "DiffHi");
136 if (!diffhi)
137 {
138 cout << "ERROR - Reading of DiffHi;avg failed." << endl;
139 return 2;
140 }
141*/
142
143 MSequence seq;
144 if (seq.Read("sequence[0-9]{8}[.]txt|MSequence")<=0)
145 {
146 cout << "ERROR - Could not find sequence in file: " << fname << endl;
147 return 2;
148 }
149 if (!seq.IsValid())
150 {
151 cout << "ERROR - Sequence read from file inavlid: " << fname << endl;
152 return 2;
153 }
154
155 TString medpuloff("NULL");
156 TString devpuloff("NULL");
157 TString medhilocal("NULL");
158 TString devhilocal("NULL");
159
160 if (seq.GetSequence()<200000 && seq.GetTelescope()==1)
161 {
162 MHCamera *hilooff = (MHCamera*)arr.FindObjectInCanvas("HiLoOff;avg", "MHCamera", "HiLoOff");
163 if (!hilooff)
164 {
165 cout << "ERROR - Reading of HiLoOff failed." << endl;
166 return 2;
167 }
168
169 MHCamera *hilocal = (MHCamera*)arr.FindObjectInCanvas("HiLoCal;avg", "MHCamera", "HiLoCal");
170 if (!hilocal)
171 {
172 cout << "ERROR - Reading of HiLoCal failed." << endl;
173 return 2;
174 }
175
176 medpuloff.Form("%7.4f", hilooff->GetMedian());
177 devpuloff.Form("%7.4f", hilooff->GetDev());
178 medhilocal.Form("%6.2f", hilocal->GetMedian());
179 devhilocal.Form("%6.2f", hilocal->GetDev());
180 }
181
182 TArrayI inner(1);
183 inner[0] = 0;
184
185 TArrayI outer(1);
186 outer[0] = 1;
187
188 Int_t s0[] = { 1, 2, 3, 4, 5, 6 };
189
190 Stat_t meanrmsi = cam->GetMeanSectors(TArrayI(6, s0), inner);
191 Stat_t meanrmso = cam->GetMeanSectors(TArrayI(6, s0), outer);
192
193 if (meanrmsi<0 || meanrmso<0)
194 {
195 cout << "ERROR - MeanPedRMS inner or outer < 0 " << endl;
196 cout << "MeanPedRMS inner " << meanrmsi << endl;
197 cout << "MeanPedRMS outer " << meanrmso << endl;
198 return 2;
199 }
200
201 TString meanrmsinner=Form("%6.2f", meanrmsi);
202 TString meanrmsouter=Form("%6.2f", meanrmso);
203
204 cam = (MHCamera*)arr.FindObjectInCanvas("Interp'd;avg", "MHCamera", "Interp'd");
205 if (!cam)
206 {
207 cout << "ERROR - Reading of Interp'd;avg failed." << endl;
208 return 2;
209 }
210
211 Stat_t meansigi = cam->GetMeanSectors(TArrayI(6, s0), inner);
212 Stat_t meansigo = cam->GetMeanSectors(TArrayI(6, s0), outer);
213
214 if (meansigi<0 || meansigo<0)
215 {
216 cout << "ERRROR - MeanInterp'd inner or outer < 0 " << endl;
217 cout << "MeanInterp'd inner " << meansigi << endl;
218 cout << "MeanInterp'd outer " << meansigo << endl;
219 return 2;
220 }
221
222 TString meansiginner =Form("%6.2f", meansigi);
223 TString meansigouter =Form("%6.2f", meansigo);
224
225 TString calpos = cal ? Form("%5.1f", cal->GetMean()) : "NULL";
226
227 if (pul->GetMean()<0 || pul->GetRMS()<0)
228 {
229 cout << "ERROR - PulsePos'd mean or rms < 0 " << endl;
230 cout << "PulsePos'd mean " << pul->GetMean() << endl;
231 cout << "PulsePos'd rms " << pul->GetRMS() << endl;
232 return 2;
233 }
234
235 TString meanpulpos = Form("%6.2f", pul->GetMean());
236 TString rmspulpos = Form("%6.2f", pul->GetRMS());
237
238/*
239 Double_t meanhi = TMath::Nint(pulhi->GetMean()*100.)/100.;
240 Double_t rmshi = TMath::Nint(pulhi->GetRMS() *100.)/100.;
241
242 Double_t meanlo = TMath::Nint(pullo->GetMean()*100.)/100.;
243 Double_t rmslo = TMath::Nint(pullo->GetRMS() *100.)/100.;
244 pullo->Add(pullo, pulhi, 1, -1);
245 pullo->ResetBit(MHCamera::kProfile);
246
247 Double_t meanoff = TMath::Nint(pullo->GetMean()*100.)/100.;
248 Double_t rmsoff = TMath::Nint(pullo->GetRMS() *100.)/100.;
249
250 // USE MEDIAN INSTEAD? GetQuantiles(Int_t nprobSum, Double_t *q, const Double_t *probSum=0);
251
252 TString meanpulhi =Form("%6.2f", meanhi);
253 TString rmspulhi =Form("%6.2f", rmshi);
254
255 TString meanpullo =Form("%6.2f", meanlo);
256 TString rmspullo =Form("%6.2f", rmslo);
257 */
258
259 cam = (MHCamera*)arr.FindObjectInCanvas("Unsuitable;avg", "MHCamera", "Unsuitable");
260 if (!cam)
261 {
262 cout << "ERROR - Reading of Unsuitable;avg failed." << endl;
263 return 2;
264 }
265
266 Int_t unsuitable50 = cam->GetNumBinsAboveThreshold(0.50);
267 Int_t unsuitable01 = cam->GetNumBinsAboveThreshold(0.01);
268
269 TString unsuitablemax = "NULL";
270 TString deadmax = "NULL";
271
272 TGraph *gr = (TGraph*)arr.FindObjectInCanvas("BadPixTm", "TGraph", "BadPixTm");
273 if (gr)
274 {
275 const Int_t p = TMath::FloorNint(gr->GetN()*0.999);
276 unsuitablemax = Form("%d", TMath::Nint(TMath::KOrdStat(gr->GetN(), gr->GetY(), p)));
277 }
278
279 gr = (TGraph*)arr.FindObjectInCanvas("DeadPixTm", "TGraph", "DeadPixTm");
280 if (gr)
281 deadmax = Form("%d", TMath::Nint(TMath::MaxElement(gr->GetN(), gr->GetY())));
282
283 TString rateped = "NULL";
284 TString rateped2 = "NULL";
285 TString ratecal = "NULL";
286 TString ratetrig = "NULL";
287 TString ratesum = "NULL";
288 TString ratena = "NULL";
289 TString ratenull = "NULL";
290
291 TH2D *htp = (TH2D*)arr.FindObjectInCanvas("TrigPat", "TH2D", "TrigPat");
292 if (htp)
293 {
294 htp->ResetBit(TH1::kCanRebin);
295
296 Int_t iped = htp->GetYaxis()->FindBin("Ped");
297 Int_t iped2 = htp->GetYaxis()->FindBin("Ped+Trig");
298 Int_t ical = htp->GetYaxis()->FindBin("Cal");
299 Int_t itrig = htp->GetYaxis()->FindBin("Trig");
300 Int_t isum = htp->GetYaxis()->FindBin("Sum");
301 Int_t inull = htp->GetYaxis()->FindBin("0");
302 Int_t ina = htp->GetYaxis()->FindBin("UNKNOWN");
303
304 Int_t nx = htp->GetNbinsX();
305
306 rateped = iped <0 ? "NULL" : Form("%8.1f", htp->Integral(1, nx, iped, iped) / nx);
307 rateped2 = iped2<0 ? "NULL" : Form("%7.2f", htp->Integral(1, nx, iped2, iped2) / nx);
308 ratecal = ical <0 ? "NULL" : Form("%8.1f", htp->Integral(1, nx, ical, ical) / nx);
309 ratetrig = itrig<0 ? "NULL" : Form("%8.1f", htp->Integral(1, nx, itrig, itrig) / nx);
310 ratesum = isum <0 ? "NULL" : Form("%8.1f", htp->Integral(1, nx, isum, isum) / nx);
311 ratenull = inull<0 ? "NULL" : Form("%8.1f", htp->Integral(1, nx, inull, inull) / nx);
312 ratena = ina <0 ? "NULL" : Form("%7.2f", htp->Integral(1, nx, ina, ina) / nx);
313 }
314
315 // *****************************************************
316
317 cout << "Sequence M" << seq.GetTelescope() << "/" << seq.GetSequence() << endl;
318 cout << " Mean Ped RMS inner [phe] " << meanrmsinner << endl;
319 cout << " Mean Ped RMS outer [phe] " << meanrmsouter << endl;
320 cout << " Mean Signal inner [phe] " << meansiginner << endl;
321 cout << " Mean Signal outer [phe] " << meansigouter << endl;
322 cout << " Mean extracted PulsePos " << meanextpulpos << " +- " << rmsextpulpos << endl;
323 cout << " Mean calibrated PulsePos " << meanpulpos << " +- " << rmspulpos << endl;
324 cout << " Mean calib pulse pos " << calpos << endl;
325// cout << " Mean ext.HiGain PulsePos " << meanpulhi << " +- " << rmspulhi << endl;
326// cout << " Mean ext.LoGain PulsePos " << meanpullo << " +- " << rmspullo << endl;
327 cout << " Lo-Hi gain offset: " << medpuloff << " +- " << devpuloff << endl;
328 cout << " Hi/Lo gain ratio: " << medhilocal << " +- " << devhilocal << endl;
329 cout << " Unsuitable > 50%: " << setw(6) << unsuitable50 << endl;
330 cout << " Unsuitable > 1%: " << setw(6) << unsuitable01 << endl;
331 cout << " UnsuitableMax (99.9%) " << setw(6) << unsuitablemax << endl;
332 cout << " DeadMax " << setw(6) << deadmax << endl;
333 cout << " Rate Trigger [Hz] " << ratetrig << endl;
334 cout << " Rate SUM [Hz] " << ratesum << endl;
335 cout << " Rate Ped+Trigger [Hz] " << rateped2 << endl;
336 cout << " Rate Pedestal [Hz] " << rateped << endl;
337 cout << " Rate Calibration [Hz] " << ratecal << endl;
338 cout << " Rate 0 [Hz] " << ratenull << endl;
339 cout << " Rate UNKNOWN [Hz] " << ratena << endl;
340 cout << endl;
341
342 //build query and insert information into the database
343 // here only a update query is built, as this macro is exexuted after
344 // the macro fillcalib.C in the script fillcallisto
345 // and so the table Calibration is always updated
346 TString vars = Form(" fMeanPedRmsInner=%s, fMeanPedRmsOuter=%s, "
347 " fMeanSignalInner=%s, fMeanSignalOuter=%s, "
348 " fPulsePosMean=%s, fPulsePosRms=%s, "
349 " fPulsePosCheckMean=%s, fPulsePosCheckRms=%s, "
350 " fPulsePosCalib=%s, "
351 //" fPulsePosHiMean=%s, fPulsePosHiRms=%s, "
352 //" fPulsePosLoMean=%s, fPulsePosLoRms=%s, "
353 " fPulsePosOffMed=%s, fPulsePosOffDev=%s, "
354 " fHiLoGainRatioMed=%s, fHiLoGainRatioDev=%s, "
355 " fUnsuitable50=%d, fUnsuitable01=%d, "
356 " fUnsuitableMax=%s, fDeadMax=%s, "
357 " fRateTrigEvts=%s, fRateSumEvts=%s, "
358 " fRatePedTrigEvts=%s, fRatePedEvts=%s, "
359 " fRateCalEvts=%s, fRateNullEvts=%s, "
360 " fRateUnknownEvts=%s ",
361 meanrmsinner.Data(), meanrmsouter.Data(),
362 meansiginner.Data(), meansigouter.Data(),
363 meanpulpos.Data(), rmspulpos.Data(),
364 meanextpulpos.Data(), rmsextpulpos.Data(),
365 calpos.Data(),
366 //meanpulhi.Data(), rmspulhi.Data(),
367 //meanpullo.Data(), rmspullo.Data(),
368 medpuloff.Data(), devpuloff.Data(),
369 medhilocal.Data(), devhilocal.Data(),
370 unsuitable50, unsuitable01,
371 unsuitablemax.Data(), deadmax.Data(),
372 ratetrig.Data(), ratesum.Data(),
373 rateped2.Data(), rateped.Data(),
374 ratecal.Data(), ratenull.Data(),
375 ratena.Data()
376 );
377
378 TString where = Form("fSequenceFirst=%d AND fTelescopeNumber=%d",
379 seq.GetSequence(), seq.GetTelescope());
380 return serv.Update("Calibration", vars, where) ? 1 : 2;
381}
382
383int fillsignal(TString fname, Bool_t dummy=kTRUE)
384{
385 MSQLMagic serv("sql.rc");
386 if (!serv.IsConnected())
387 {
388 cout << "ERROR - Connection to database failed." << endl;
389 return 0;
390 }
391
392 cout << "fillsignal" << endl;
393 cout << "----------" << endl;
394 cout << endl;
395 cout << "Connected to " << serv.GetName() << endl;
396 cout << "File: " << fname << endl;
397 cout << endl;
398
399 serv.SetIsDummy(dummy);
400
401 //process file
402 return Process(serv, fname);
403}
Note: See TracBrowser for help on using the repository browser.