source: trunk/MagicSoft/Mars/datacenter/macros/fillstar.C@ 8760

Last change on this file since 8760 was 8694, checked in by msmeyer, 18 years ago
*** empty log message ***
File size: 12.2 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 = 70.205*mw - 28.055;
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()/7216)*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 if (sparkrate.Contains("inf"))
214 sparkrate="0.00";
215
216 TGraph *g = (TGraph*)arr.FindObjectInCanvas("Humidity", "TGraph", "MHWeather");
217 if (!g)
218 {
219 cout << "WARNING - Reading of Humidity failed." << endl;
220 return 2;
221 }
222
223 Double_t max = TMath::MaxElement(g->GetN(), g->GetY());
224 TString maxhum = Form("%6.1f", max);
225
226
227 g = (TGraph*)arr.FindObjectInCanvas("NumStars", "TGraph", "MHPointing");
228 if (!g)
229 {
230 cout << "WARNING - Reading of NumStars failed." << endl;
231// return 2;
232 }
233
234 Double_t numstarmed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
235 TString numstarsmed = Form("%5.1f", numstarmed);
236 Double_t numstarrms = g ? g->GetRMS(2) : -1;
237 TString numstarsrms = Form("%5.1f", numstarrms);
238
239 g = (TGraph*)arr.FindObjectInCanvas("NumStarsCor", "TGraph", "MHPointing");
240 if (!g)
241 {
242 cout << "WARNING - Reading of NumStarsCor failed." << endl;
243// return 2;
244 }
245
246 Double_t numcormed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
247 TString numcorsmed = Form("%5.1f", numcormed);
248 Double_t numcorrms = g ? g->GetRMS(2) : -1;
249 TString numcorsrms = Form("%5.1f", numcorrms);
250
251 g = (TGraph*)arr.FindObjectInCanvas("Brightness", "TGraph", "MHPointing");
252 if (!g)
253 {
254 cout << "WARNING - Reading of SkyBrightness failed." << endl;
255// return 2;
256 }
257
258 Double_t brightnessmed = g ? TMath::Median(g->GetN(), g->GetY()) : -1;
259 TString skybrightnessmed = Form("%5.1f", brightnessmed);
260 Double_t brightnessrms = g ? g->GetRMS(2) : -1;
261 TString skybrightnessrms = Form("%5.1f", brightnessrms);
262
263
264 TString sequence = fname(TRegexp("star[0-9]+[.]root$"));
265 if (sequence.IsNull())
266 {
267 cout << "WARNING - Could not determine sequence# from filename: " << fname << endl;
268 return 2;
269 }
270
271 Int_t seq = atoi(sequence.Data()+5);
272
273 cout << "Sequence #" << seq << endl;
274 cout << " Inhomogeneity " << inhomogen << endl;
275 cout << " PSF [mm] " << PSF << endl;
276 cout << " Island Quality " << islands << endl;
277 cout << " Ratio [%] " << ratio << endl;
278 cout << " Muon Number " << num << endl;
279 cout << " EffOnTime [s] " << effontime << endl;
280 cout << " Muon Rate [Hz] " << muonrate << endl;
281 cout << " # of Events (w/o sparks) " << numevents << endl;
282 cout << " # of Sparks " << numsparks << endl;
283 cout << " Rate after ImgCl [Hz] " << datarate << endl;
284 cout << " Rate of sparks [Hz] " << sparkrate << endl;
285 cout << " Maximum Humidity [%] " << maxhum << endl;
286 cout << " Number of Stars " << numstarsmed << " +/- " << numstarsrms << endl;
287 cout << " Number of cor. Stars " << numcormed << " +/- " << numcorrms << endl;
288 cout << " Skybrightness " << skybrightnessmed << " +/- " << skybrightnessrms << endl;
289
290 TString query;
291 if (!ExistStr(serv, "fSequenceFirst", "Star", seq))
292 {
293 query = Form("INSERT Star SET"
294 " fSequenceFirst=%d,"
295 " fMeanNumberIslands=%s, "
296 " fRatio=%s, "
297 " fMuonNumber=%d, "
298 " fEffOnTime=%d, "
299 " fMuonRate=%s, "
300 " fPSF=%s, "
301 " fDataRate=%d, "
302 " fSparkRate=%s, "
303 " fMaxHumidity=%s ,"
304 " fNumStarsMed=%s ,"
305 " fNumStarsRMS=%s ,"
306 " fNumStarsCorMed=%s ,"
307 " fNumStarsCorRMS=%s ,"
308 " fBrightnessMed=%s ,"
309 " fBrightnessRMS=%s ,"
310 " fInhomogeneity=%s ",
311 seq, islands.Data(), ratio.Data(),
312 num, effontime,
313 muonrate.Data(), PSF.Data(),
314 datarate, sparkrate.Data(), maxhum.Data(),
315 numstarsmed.Data(), numstarsrms.Data(),
316 numcorsmed.Data(), numcorsrms.Data(),
317 skybrightnessmed.Data(), skybrightnessrms.Data(),
318 inhomogen.Data());
319 }
320 else
321 {
322 query = Form("UPDATE Star SET"
323 " fMeanNumberIslands=%s, "
324 " fRatio=%s, "
325 " fMuonNumber=%d, "
326 " fEffOnTime=%d, "
327 " fMuonRate=%s, "
328 " fPSF=%s, "
329 " fDataRate=%d, "
330 " fSparkRate=%s, "
331 " fMaxHumidity=%s, "
332 " fNumStarsMed=%s ,"
333 " fNumStarsRMS=%s ,"
334 " fNumStarsCorMed=%s ,"
335 " fNumStarsCorRMS=%s ,"
336 " fBrightnessMed=%s ,"
337 " fBrightnessRMS=%s ,"
338 " fInhomogeneity=%s "
339 " WHERE fSequenceFirst=%d ",
340 islands.Data(), ratio.Data(),
341 num, effontime,
342 muonrate.Data(), PSF.Data(),
343 datarate, sparkrate.Data(), maxhum.Data(),
344 numstarsmed.Data(), numstarsrms.Data(),
345 numcorsmed.Data(), numcorsrms.Data(),
346 skybrightnessmed.Data(), skybrightnessrms.Data(),
347 inhomogen.Data(), seq);
348 }
349
350// cout << "Q: " << query << endl;
351
352 if (dummy)
353 return 1;
354
355 TSQLResult *res = serv.Query(query);
356 if (!res)
357 {
358 cout << "ERROR - Query failed: " << query << endl;
359 return 2;
360 }
361 delete res;
362
363 return 1;
364}
365
366int fillstar(TString fname, Bool_t dummy=kTRUE)
367{
368 TEnv env("sql.rc");
369
370 MSQLServer serv(env);
371 if (!serv.IsConnected())
372 {
373 cout << "ERROR - Connection to database failed." << endl;
374 return 0;
375 }
376
377 cout << "fillstar" << endl;
378 cout << "---------" << endl;
379 cout << endl;
380 cout << "Connected to " << serv.GetName() << endl;
381 cout << "File: " << fname << endl;
382 cout << endl;
383
384 return Process(serv, fname, dummy);
385}
Note: See TracBrowser for help on using the repository browser.