source: branches/Corsika7500Compatibility/datacenter/macros/fillmcsignal.C

Last change on this file was 9225, checked in by hoehne, 16 years ago
*** empty log message ***
File size: 10.6 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! Author(s): Daniel Hoehne-Moench, 01/2009 <mailto:hoehne@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2009
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// fillmcsignal.C
30// ============
31//
32// This macro is used to read the calibration-/callisto-output files
33// signal00000.root.
34//
35// From this file the mean pedrms, the mean signal and the pulse position
36// for the inner and outer camera is extracted and inserted into the database
37// in the table MCCalibration, where the results of callisto are stored.
38// The sequence number is extracted from the filename.
39//
40// Usage:
41// .x fillmcsignal.C("/magic/montecarlo/callisto/0002/00028270/signal00028270.root", kTRUE)
42//
43// The second argument is the 'dummy-mode'. If it is kTRUE dummy-mode is
44// switched on and nothing will be written into the database. This is usefull
45// for tests.
46//
47// The macro can also be run without ACLiC but this is a lot slower...
48//
49// Remark: Running it from the commandline looks like this:
50// root -q -l -b fillmcsignal.C+\(\"filename\"\,kFALSE\) 2>&1 | tee fillsignal.log
51//
52// Make sure, that database and password are corretly set in a resource
53// file called sql.rc and the resource file is found.
54//
55// Returns 2 in case of failure, 1 in case of success and 0 if the connection
56// to the database is not working.
57//
58/////////////////////////////////////////////////////////////////////////////
59#include <iostream>
60#include <iomanip>
61
62#include <TRegexp.h>
63
64#include <TFile.h>
65#include <TGraph.h>
66#include <TSQLResult.h>
67#include <TSQLRow.h>
68
69#include "MSQLServer.h"
70#include "MSQLMagic.h"
71
72#include "MStatusArray.h"
73#include "MSequence.h"
74#include "MHCamera.h"
75#include "MHVsTime.h"
76
77#include "MCalibrationPulseTimeCam.h"
78#include "MCalibrationPix.h"
79
80using namespace std;
81
82int Process(MSQLMagic &serv, TString fname)
83{
84 TFile file(fname, "READ");
85 if (!file.IsOpen())
86 {
87 cout << "ERROR - Could not find file " << fname << endl;
88 return 2;
89 }
90
91 TString meanextpulpos("NULL");
92 TString rmsextpulpos("NULL");
93
94 MCalibrationPulseTimeCam *pt;
95 file.GetObject("MCalibrationPulseTimeCam", pt);
96 if (pt)
97 {
98 Double_t meanextpul = pt->GetAverageArea(0).GetHiGainMean();
99 Double_t rmsextpul = pt->GetAverageArea(0).GetHiGainRms();
100
101 if (meanextpul>=0 || rmsextpulpos>=0)
102 {
103 meanextpulpos.Form("%6.2f", meanextpul);
104 rmsextpulpos.Form("%6.2f", rmsextpul);
105 }
106 }
107
108
109 MStatusArray arr;
110 if (arr.Read()<=0)
111 {
112 cout << "ERROR - Reading of MStatusDisplay failed." << endl;
113 return 2;
114 }
115
116 MHCamera *cam = (MHCamera*)arr.FindObjectInCanvas("PedRMS;avg", "MHCamera", "PedRMS");
117 if (!cam)
118 {
119 cout << "ERROR - Reading of PedRMS;avg failed." << endl;
120 return 2;
121 }
122
123 MHCamera *cal = (MHCamera*)arr.FindObjectInCanvas("CalPos;avg", "MHCamera", "CalPos");
124 MHCamera *pul = (MHCamera*)arr.FindObjectInCanvas("PulsePos;avg", "MHCamera", "PulsePos");
125 if (!pul)
126 {
127 cout << "ERROR - Reading of PulsePos;avg failed." << endl;
128 return 2;
129 }
130
131 MSequence seq;
132 if (seq.Read("sequence[0-9]{8}[.]txt|MSequence")<=0)
133 {
134 cout << "ERROR - Could not find sequence in file: " << fname << endl;
135 return 2;
136 }
137 if (!seq.IsValid())
138 {
139 cout << "ERROR - Sequence read from file inavlid: " << fname << endl;
140 return 2;
141 }
142
143 TString medpuloff("NULL");
144 TString devpuloff("NULL");
145 TString medhilocal("NULL");
146 TString devhilocal("NULL");
147
148 TString query(Form("SELECT MCSequences.fAmplFadcKEY FROM MCSequences LEFT JOIN "
149 "AmplFadc on MCSequences.fAmplFadcKEY=AmplFadc.fAmplFadcKEY "
150 "where fSequenceFirst=%d;",seq.GetSequence()));
151
152 TSQLResult *res = serv.Query(query);
153 if (!res)
154 {
155 cout << "ERROR - Query failed: " << query << endl;
156 return 2;
157 }
158
159 TSQLRow *row = res->Next();
160 if (!row)
161 {
162 cout << "ERROR - Query failed: " << query << endl;
163 return 2;
164 }
165 TString epochkey = (*row)[0];
166 Int_t epoch = atoi(epochkey.Data());
167
168 delete res;
169
170 if (epoch!=1)
171 {
172 MHCamera *hilooff = (MHCamera*)arr.FindObjectInCanvas("HiLoOff;avg", "MHCamera", "HiLoOff");
173 if (!hilooff)
174 {
175 cout << "ERROR - Reading of HiLoOff failed." << endl;
176 return 2;
177 }
178
179 MHCamera *hilocal = (MHCamera*)arr.FindObjectInCanvas("HiLoCal;avg", "MHCamera", "HiLoCal");
180 if (!hilocal)
181 {
182 cout << "ERROR - Reading of HiLoCal failed." << endl;
183 return 2;
184 }
185
186 medpuloff.Form("%7.4f", hilooff->GetMedian());
187 devpuloff.Form("%7.4f", hilooff->GetDev());
188 medhilocal.Form("%6.2f", hilocal->GetMedian());
189 devhilocal.Form("%6.2f", hilocal->GetDev());
190 }
191
192 TArrayI inner(1);
193 inner[0] = 0;
194
195 TArrayI outer(1);
196 outer[0] = 1;
197
198 Int_t s0[] = { 1, 2, 3, 4, 5, 6 };
199
200 Stat_t meanrmsi = cam->GetMeanSectors(TArrayI(6, s0), inner);
201 Stat_t meanrmso = cam->GetMeanSectors(TArrayI(6, s0), outer);
202
203 if (meanrmsi<0 || meanrmso<0)
204 {
205 cout << "ERROR - MeanPedRMS inner or outer < 0 " << endl;
206 cout << "MeanPedRMS inner " << meanrmsi << endl;
207 cout << "MeanPedRMS outer " << meanrmso << endl;
208 return 2;
209 }
210
211 TString meanrmsinner=Form("%6.2f", meanrmsi);
212 TString meanrmsouter=Form("%6.2f", meanrmso);
213
214 cam = (MHCamera*)arr.FindObjectInCanvas("Interp'd;avg", "MHCamera", "Interp'd");
215 if (!cam)
216 {
217 cout << "ERROR - Reading of Interp'd;avg failed." << endl;
218 return 2;
219 }
220
221 Stat_t meansigi = cam->GetMeanSectors(TArrayI(6, s0), inner);
222 Stat_t meansigo = cam->GetMeanSectors(TArrayI(6, s0), outer);
223
224 if (meansigi<0 || meansigo<0)
225 {
226 cout << "ERROR - MeanInterp'd inner or outer < 0 " << endl;
227 cout << "MeanInterp'd inner " << meansigi << endl;
228 cout << "MeanInterp'd outer " << meansigo << endl;
229 return 2;
230 }
231
232 TString meansiginner =Form("%6.2f", meansigi);
233 TString meansigouter =Form("%6.2f", meansigo);
234
235 TString calpos = cal ? Form("%5.1f", cal->GetMean()) : "NULL";
236
237 if (pul->GetMean()<0 || pul->GetRMS()<0)
238 {
239 cout << "ERROR - PulsePos'd mean or rms < 0 " << endl;
240 cout << "PulsePos'd mean " << pul->GetMean() << endl;
241 cout << "PulsePos'd rms " << pul->GetRMS() << endl;
242 return 2;
243 }
244
245 TString meanpulpos = Form("%6.2f", pul->GetMean());
246 TString rmspulpos = Form("%6.2f", pul->GetRMS());
247
248 cam = (MHCamera*)arr.FindObjectInCanvas("Unsuitable;avg", "MHCamera", "Unsuitable");
249 if (!cam)
250 {
251 cout << "ERROR - Reading of Unsuitable;avg failed." << endl;
252 return 2;
253 }
254
255 Int_t unsuitable50 = cam->GetNumBinsAboveThreshold(0.50);
256 Int_t unsuitable01 = cam->GetNumBinsAboveThreshold(0.01);
257
258 // *****************************************************
259
260 cout << "Sequence #" << seq.GetSequence() << endl;
261 cout << " Mean Ped RMS inner [phe] " << meanrmsinner << endl;
262 cout << " Mean Ped RMS outer [phe] " << meanrmsouter << endl;
263 cout << " Mean Signal inner [phe] " << meansiginner << endl;
264 cout << " Mean Signal outer [phe] " << meansigouter << endl;
265 cout << " Mean extracted PulsePos " << meanextpulpos << " +- " << rmsextpulpos << endl;
266 cout << " Mean calibrated PulsePos " << meanpulpos << " +- " << rmspulpos << endl;
267 cout << " Mean calib pulse pos " << calpos << endl;
268 cout << " Lo-Hi gain offset: " << medpuloff << " +- " << devpuloff << endl;
269 cout << " Hi/Lo gain ratio: " << medhilocal << " +- " << devhilocal << endl;
270 cout << " Unsuitable > 50%: " << setw(6) << unsuitable50 << endl;
271 cout << " Unsuitable > 1%: " << setw(6) << unsuitable01 << endl;
272 cout << endl;
273
274 //build query and insert information into the database
275 // here only a update query is built, as this macro is exexuted after
276 // the macro fillmccalib.C in the script fillmccallisto
277 // and so the table MCCalibration is always updated
278 TString vars = Form(" fMeanPedRmsInner=%s, fMeanPedRmsOuter=%s, "
279 " fMeanSignalInner=%s, fMeanSignalOuter=%s, "
280 " fPulsePosMean=%s, fPulsePosRms=%s, "
281 " fPulsePosCheckMean=%s, fPulsePosCheckRms=%s, "
282 " fPulsePosOffMed=%s, fPulsePosOffDev=%s, "
283 " fHiLoGainRatioMed=%s, fHiLoGainRatioDev=%s, "
284 " fUnsuitable50=%d, fUnsuitable01=%d, "
285 " fPulsePosCalib=%s ",
286 meanrmsinner.Data(), meanrmsouter.Data(),
287 meansiginner.Data(), meansigouter.Data(),
288 meanpulpos.Data(), rmspulpos.Data(),
289 meanextpulpos.Data(), rmsextpulpos.Data(),
290 medpuloff.Data(), devpuloff.Data(),
291 medhilocal.Data(), devhilocal.Data(),
292 unsuitable50, unsuitable01,
293 calpos.Data());
294
295 TString where = Form("fSequenceFirst=%d", seq.GetSequence());
296 return serv.Update("MCCalibration", vars, where) ? 1 : 2;
297
298}
299
300int fillmcsignal(TString fname, Bool_t dummy=kTRUE)
301{
302 MSQLMagic serv("sql.rc");
303 if (!serv.IsConnected())
304 {
305 cout << "ERROR - Connection to database failed." << endl;
306 return 0;
307 }
308
309 cout << "fillmcsignal" << endl;
310 cout << "----------" << endl;
311 cout << endl;
312 cout << "Connected to " << serv.GetName() << endl;
313 cout << "File: " << fname << endl;
314 cout << endl;
315
316 serv.SetIsDummy(dummy);
317
318 //process file
319 return Process(serv, fname);
320}
Note: See TracBrowser for help on using the repository browser.