source: tags/Mars-V1.1/datacenter/macros/fillstar.C

Last change on this file was 8332, checked in by hoehne, 18 years ago
*** empty log message ***
File size: 12.1 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, 05/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Daniela Dorner, 05/2005 <mailto:dorner@astro.uni-wuerzburg.de>
20!
21! Copyright: MAGIC Software Development, 2000-2006
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// fillstar.C
29// ==========
30//
31// This macro is used to read the star-output files.
32// These files are automatically called star00000000.root.
33// From these files the muon parameters (psf, muon number, ratio, muon rate),
34// the rate, the number of islands, the effective ontime, the maximum
35// humidity and a parameter describing the inhomogeneity are extracted from
36// the status display.
37// The sequence number is extracted from the filename.
38//
39// Usage:
40// .x fillstar.C("/magic/data/star/0004/00047421/star00047421.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 fillstar.C+\(\"filename\"\,kFALSE\) 2>&1 | tee fillstar.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 <TEnv.h>
62#include <TRegexp.h>
63
64#include <TH1.h>
65#include <TH2.h>
66#include <TGraph.h>
67#include <TProfile.h>
68#include <TFile.h>
69#include <TSQLResult.h>
70#include <TSQLRow.h>
71
72#include "MSQLServer.h"
73
74#include "MHCamera.h"
75#include "MHMuonPar.h"
76#include "MStatusArray.h"
77#include "MGeomCamMagic.h"
78#include "MBadPixelsCam.h"
79
80using namespace std;
81
82// --------------------------------------------------------------------------
83//
84// Checks whether an entry is already existing
85//
86Bool_t ExistStr(MSQLServer &serv, const char *column, const char *table, Int_t test)
87{
88 TString query(Form("SELECT %s FROM %s WHERE %s='%d'", column, table, column, test));
89 TSQLResult *res = serv.Query(query);
90 if (!res)
91 return kFALSE;
92
93 TSQLRow *row;
94
95 Bool_t rc = kFALSE;
96 while ((row=res->Next()))
97 {
98 if ((*row)[0])
99 {
100 rc = kTRUE;
101 break;
102 }
103 }
104
105 delete res;
106
107 return rc;
108}
109
110
111int Process(MSQLServer &serv, TString fname, Bool_t dummy)
112{
113 TFile file(fname, "READ");
114 if (!file.IsOpen())
115 {
116 cout << "ERROR - Could not find file " << fname << endl;
117 return 2;
118 }
119
120
121 MStatusArray arr;
122 if (arr.Read()<=0)
123 {
124 cout << "ERROR - Reading of MStatusDisplay failed." << endl;
125 return 2;
126 }
127
128 MHCamera *hsparks = (MHCamera*)arr.FindObjectInCanvas("Sparks;avg", "MHCamera", "Sparks");
129 if (!hsparks)
130 {
131 cout << "WARNING - Reading of Sparks failed." << endl;
132 return 2;
133 }
134
135 TH2F *hcog = (TH2F*)arr.FindObjectInCanvas("Center", "TH2F", "MHHillas");
136 if (!hcog)
137 {
138 cout << "WARNING - Reading of MHHillas failed." << endl;
139 return 2;
140 }
141
142 MHMuonPar *hmuon = (MHMuonPar*)arr.FindObjectInCanvas("MHMuonPar", "MHMuonPar", "MHMuonPar");
143 if (!hmuon)
144 {
145 cout << "WARNING - Reading of MHMuon failed." << endl;
146 return 2;
147 }
148
149 Double_t val[6];
150 for (int x=1; x<hcog->GetNbinsX(); x++)
151 for (int y=1; y<hcog->GetNbinsY(); y++)
152 {
153 Stat_t px = hcog->GetXaxis()->GetBinCenter(x);
154 Stat_t py = hcog->GetYaxis()->GetBinCenter(y);
155 Int_t i = (TMath::Nint(3*TMath::ATan2(px,py)/TMath::Pi())+6)%6;
156 val[i] += hcog->GetBinContent(x, y);
157 }
158
159 Double_t inhom = TMath::RMS(6, val)*6/hcog->GetEntries()*100;
160 inhom = TMath::Nint(inhom*10)/10.;
161 TString inhomogen = Form("%5.1f", inhom);
162
163 Float_t mw = hmuon->GetMeanWidth();
164 Float_t psf = 58.5480*mw - 37.3829;
165 psf = TMath::Nint(psf*10)/10.;
166 if (psf<0)
167 psf=0;
168 TString PSF = Form("%5.1f", psf);
169 Int_t num = (int)hmuon->GetEntries();
170
171 Float_t ratiodatamc = (hmuon->GetMeanSize()/9327.)*100;
172 TString ratio = Form("%5.1f", ratiodatamc);
173
174 TH1 *h = (TH1*)arr.FindObjectInCanvas("Islands", "TH1F", "MHImagePar");
175 if (!h)
176 {
177 cout << "WARNING - Reading of Islands failed." << endl;
178 return 2;
179 }
180
181 Float_t quality = h->GetMean();
182 quality = TMath::Nint(quality*100)/100.;
183 TString islands = Form("%6.2f", quality);
184
185 h = (TH1*)arr.FindObjectInCanvas("EffOnTheta", "TH1D", "EffOnTime");
186 if (!h)
187 {
188 cout << "WARNING - Reading of EffOnTheta failed." << endl;
189 return 2;
190 }
191
192 Float_t effon = h->Integral();
193 Float_t mrate = num/effon;
194 mrate = TMath::Nint(mrate*100)/100.;
195 if (mrate<0)
196 mrate=0;
197 TString muonrate = Form("%6.2f", mrate);
198 Int_t effontime = TMath::Nint(effon);
199
200 h = (TH1*)arr.FindObjectInCanvas("ProjDeltaT", "TH1D", "EffOnTime");
201 if (!h)
202 {
203 cout << "WARNING - Reading of ProjDeltaT failed." << endl;
204 return 2;
205 }
206
207 Int_t numsparks = (int)hsparks->GetEntries();
208 Int_t numevents = (int)h->GetEntries() - numsparks;
209 Int_t datarate = (int)(numevents/effon);
210
211 TString sparkrate = Form("%5.2f", numsparks/effon);
212
213 TGraph *g = (TGraph*)arr.FindObjectInCanvas("Humidity", "TGraph", "MHWeather");
214 if (!g)
215 {
216 cout << "WARNING - Reading of Humidity failed." << endl;
217 return 2;
218 }
219
220 Double_t max = TMath::MaxElement(g->GetN(), g->GetY());
221 TString maxhum = Form("%6.1f", max);
222
223
224 g = (TGraph*)arr.FindObjectInCanvas("NumStars", "TGraph", "MHPointing");
225 if (!g)
226 {
227 cout << "WARNING - Reading of NumStars failed." << endl;
228// return 2;
229 }
230
231 Double_t numstarmed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
232 TString numstarsmed = Form("%5.1f", numstarmed);
233 Double_t numstarrms = g ? g->GetRMS(2) : -1;
234 TString numstarsrms = Form("%5.1f", numstarrms);
235
236 g = (TGraph*)arr.FindObjectInCanvas("NumStarsCor", "TGraph", "MHPointing");
237 if (!g)
238 {
239 cout << "WARNING - Reading of NumStarsCor failed." << endl;
240// return 2;
241 }
242
243 Double_t numcormed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
244 TString numcorsmed = Form("%5.1f", numcormed);
245 Double_t numcorrms = g ? g->GetRMS(2) : -1;
246 TString numcorsrms = Form("%5.1f", numcorrms);
247
248 g = (TGraph*)arr.FindObjectInCanvas("Brightness", "TGraph", "MHPointing");
249 if (!g)
250 {
251 cout << "WARNING - Reading of SkyBrightness failed." << endl;
252// return 2;
253 }
254
255 Double_t brightnessmed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
256 TString skybrightnessmed = Form("%5.1f", brightnessmed);
257 Double_t brightnessrms = g ? g->GetRMS(2) : -1;
258 TString skybrightnessrms = Form("%5.1f", brightnessrms);
259
260
261 TString sequence = fname(TRegexp("star[0-9]+[.]root$"));
262 if (sequence.IsNull())
263 {
264 cout << "WARNING - Could not determine sequence# from filename: " << fname << endl;
265 return 2;
266 }
267
268 Int_t seq = atoi(sequence.Data()+5);
269
270 cout << "Sequence #" << seq << endl;
271 cout << " Inhomogeneity " << inhomogen << endl;
272 cout << " PSF [mm] " << PSF << endl;
273 cout << " Island Quality " << islands << endl;
274 cout << " Ratio [%] " << ratio << endl;
275 cout << " Muon Number " << num << endl;
276 cout << " EffOnTime [s] " << effontime << endl;
277 cout << " Muon Rate [Hz] " << muonrate << endl;
278 cout << " # of Events (w/o sparks) " << numevents << endl;
279 cout << " # of Sparks " << numsparks << endl;
280 cout << " Rate after ImgCl [Hz] " << datarate << endl;
281 cout << " Rate of sparks [Hz] " << sparkrate << endl;
282 cout << " Maximum Humidity [%] " << maxhum << endl;
283 cout << " Number of Stars " << numstarsmed << " +/- " << numstarsrms << endl;
284 cout << " Number of cor. Stars " << numcormed << " +/- " << numcorrms << endl;
285 cout << " Skybrightness " << skybrightnessmed << " +/- " << skybrightnessrms << endl;
286
287 TString query;
288 if (!ExistStr(serv, "fSequenceFirst", "Star", seq))
289 {
290 query = Form("INSERT Star SET"
291 " fSequenceFirst=%d,"
292 " fMeanNumberIslands=%s, "
293 " fRatio=%s, "
294 " fMuonNumber=%d, "
295 " fEffOnTime=%d, "
296 " fMuonRate=%s, "
297 " fPSF=%s, "
298 " fDataRate=%d, "
299 " fSparkRate=%s, "
300 " fMaxHumidity=%s ,"
301 " fNumStarsMed=%s ,"
302 " fNumStarsRMS=%s ,"
303 " fNumStarsCorMed=%s ,"
304 " fNumStarsCorRMS=%s ,"
305 " fBrightnessMed=%s ,"
306 " fBrightnessRMS=%s ,"
307 " fInhomogeneity=%s ",
308 seq, islands.Data(), ratio.Data(),
309 num, effontime,
310 muonrate.Data(), PSF.Data(),
311 datarate, sparkrate.Data(), maxhum.Data(),
312 numstarsmed.Data(), numstarsrms.Data(),
313 numcorsmed.Data(), numcorsrms.Data(),
314 skybrightnessmed.Data(), skybrightnessrms.Data(),
315 inhomogen.Data());
316 }
317 else
318 {
319 query = Form("UPDATE Star SET"
320 " fMeanNumberIslands=%s, "
321 " fRatio=%s, "
322 " fMuonNumber=%d, "
323 " fEffOnTime=%d, "
324 " fMuonRate=%s, "
325 " fPSF=%s, "
326 " fDataRate=%d, "
327 " fSparkRate=%s, "
328 " fMaxHumidity=%s, "
329 " fNumStarsMed=%s ,"
330 " fNumStarsRMS=%s ,"
331 " fNumStarsCorMed=%s ,"
332 " fNumStarsCorRMS=%s ,"
333 " fBrightnessMed=%s ,"
334 " fBrightnessRMS=%s ,"
335 " fInhomogeneity=%s "
336 " WHERE fSequenceFirst=%d ",
337 islands.Data(), ratio.Data(),
338 num, effontime,
339 muonrate.Data(), PSF.Data(),
340 datarate, sparkrate.Data(), maxhum.Data(),
341 numstarsmed.Data(), numstarsrms.Data(),
342 numcorsmed.Data(), numcorsrms.Data(),
343 skybrightnessmed.Data(), skybrightnessrms.Data(),
344 inhomogen.Data(), seq);
345 }
346
347 cout << "Q: " << query << endl;
348
349 if (dummy)
350 return 1;
351
352 TSQLResult *res = serv.Query(query);
353 if (!res)
354 {
355 cout << "ERROR - Query failed: " << query << endl;
356 return 2;
357 }
358 delete res;
359
360 return 1;
361}
362
363int fillstar(TString fname, Bool_t dummy=kTRUE)
364{
365 TEnv env("sql.rc");
366
367 MSQLServer serv(env);
368 if (!serv.IsConnected())
369 {
370 cout << "ERROR - Connection to database failed." << endl;
371 return 0;
372 }
373
374 cout << "fillstar" << endl;
375 cout << "---------" << endl;
376 cout << endl;
377 cout << "Connected to " << serv.GetName() << endl;
378 cout << "File: " << fname << endl;
379 cout << endl;
380
381 return Process(serv, fname, dummy);
382}
Note: See TracBrowser for help on using the repository browser.