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

Last change on this file since 18656 was 9389, checked in by tbretz, 16 years ago
*** empty log message ***
File size: 14.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, 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-2008
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 <TFile.h>
62#include <TH1.h>
63#include <TH2.h>
64#include <TGraph.h>
65
66#include "MSQLMagic.h"
67
68#include "MHCamera.h"
69#include "MHMuonPar.h"
70#include "MSequence.h"
71#include "MStatusArray.h"
72#include "MBadPixelsCam.h"
73#include "MHEffectiveOnTime.h"
74
75using namespace std;
76
77bool CheckGraph(const TGraph *g)
78{
79 return g && g->GetN()>0 && !(g->GetN()==1 && g->GetX()[0]==0);
80}
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
92 MStatusArray arr;
93 if (arr.Read()<=0)
94 {
95 cout << "ERROR - Reading of MStatusDisplay failed." << endl;
96 return 2;
97 }
98
99 MHCamera *hsparks = (MHCamera*)arr.FindObjectInCanvas("Sparks;avg", "MHCamera", "Sparks");
100 if (!hsparks)
101 {
102 cout << "ERROR - Reading of Sparks failed." << endl;
103 return 2;
104 }
105
106 TH2F *hcog = (TH2F*)arr.FindObjectInCanvas("Center", "TH2F", "MHHillas");
107 if (!hcog)
108 {
109 cout << "ERROR - Reading of MHHillas failed." << endl;
110 return 2;
111 }
112
113 MHMuonPar *hmuon = (MHMuonPar*)arr.FindObjectInCanvas("MHMuonPar", "MHMuonPar", "MHMuonPar");
114 if (!hmuon)
115 {
116 cout << "ERROR - Reading of MHMuon failed." << endl;
117 return 2;
118 }
119
120 Double_t val[6];
121 for (int x=1; x<hcog->GetNbinsX(); x++)
122 for (int y=1; y<hcog->GetNbinsY(); y++)
123 {
124 Stat_t px = hcog->GetXaxis()->GetBinCenter(x);
125 Stat_t py = hcog->GetYaxis()->GetBinCenter(y);
126 Int_t i = (TMath::Nint(3*TMath::ATan2(px,py)/TMath::Pi())+6)%6;
127 val[i] += hcog->GetBinContent(x, y);
128 }
129
130 Float_t inhom = TMath::RMS(6, val)*6/hcog->GetEntries()*100;
131 TString inhomogen = Form("%5.1f", inhom);
132
133 Float_t mw = hmuon->GetMeanWidth();
134 Float_t psf = 70.205*mw - 28.055;
135 TString PSF = Form("%5.1f", (psf<0?0:psf));
136 Int_t num = (int)hmuon->GetEntries();
137
138 Float_t ratiodatamc = (hmuon->GetMeanSize()/7216)*100;
139 TString ratio = Form("%5.1f", ratiodatamc);
140
141 TH1 *h = (TH1*)arr.FindObjectInCanvas("Islands", "TH1F", "MHImagePar");
142 if (!h)
143 {
144 cout << "ERROR - Reading of Islands failed." << endl;
145 return 2;
146 }
147
148 TString islands = Form("%6.2f", h->GetMean());
149
150 h = (TH1*)arr.FindObjectInCanvas("EffOnTheta", "TH1D", "EffOnTime");
151
152 Float_t effon = h ? h->Integral() : -1;
153
154 h = (TH1*)arr.FindObjectInCanvas("Cleaned;avg", "TH1D", "Cleaned");
155 if (!h)
156 {
157 cout << "ERROR - Reading of Cleaned;avg failed." << endl;
158 return 2;
159 }
160
161 if (effon<0)
162 effon = h->GetEntries()/100;
163
164 Float_t mrate = num/effon;
165
166 TString muonrate = mrate<0 || effon==0 ? "NULL" : Form("%6.2f", mrate);
167 TString effontime = Form("%.1f", effon);
168
169 Int_t numsparks = (int)hsparks->GetEntries();
170 Int_t numevents = (int)h->GetEntries() - numsparks;
171
172 TString datarate = effon==0 ? "NULL" : Form("%.0f", numevents/effon);
173 TString sparkrate = Form("%5.2f", numsparks/effon);
174 if (sparkrate.Contains("inf") || sparkrate.Contains("nan"))
175 sparkrate="NULL";
176
177 MHEffectiveOnTime *ontm = (MHEffectiveOnTime*)arr.FindObjectInCanvas("MHEffectiveOnTime", "MHEffectiveOnTime", "EffOnTime");
178 TString totontime = ontm ? Form("%d", TMath::Nint(ontm->GetTotalTime())) : "NULL";
179 TString relontime = ontm ? Form("%.2f", effon/ontm->GetTotalTime()*100) : "NULL";
180
181 TGraph *g = (TGraph*)arr.FindObjectInCanvas("MCamEvent", "TGraph", "Currents");
182
183 TString maxdc = CheckGraph(g) ? Form("%5.2f", TMath::MaxElement(g->GetN(), g->GetY())) : "NULL";
184 TString mindc = CheckGraph(g) ? Form("%5.2f", TMath::MinElement(g->GetN(), g->GetY())) : "NULL";
185 TString meddc = CheckGraph(g) ? Form("%5.2f", TMath::Median(g->GetN(), g->GetY())) : "NULL";
186
187 g = (TGraph*)arr.FindObjectInCanvas("Humidity", "TGraph", "MHWeather");
188
189 TString maxhum = CheckGraph(g) ? Form("%5.1f", TMath::MaxElement(g->GetN(), g->GetY())) : "NULL";
190 TString avghum = CheckGraph(g) ? Form("%5.1f", g->GetMean(2)) : "NULL";
191
192 g = (TGraph*)arr.FindObjectInCanvas("Temperature", "TGraph", "MHWeather");
193
194 TString avgtemp = CheckGraph(g) ? Form("%5.1f", g->GetMean(2)) : "NULL";
195
196 g = (TGraph*)arr.FindObjectInCanvas("WindSpeed", "TGraph", "MHWeather");
197
198 TString avgwind = CheckGraph(g) ? Form("%5.1f", g->GetMean(2)) : "NULL";
199
200 g = (TGraph*)arr.FindObjectInCanvas("Cloudiness", "TGraph", "MHWeather");
201 //if (!g)
202 // cout << "WARNING - Reading of Cloudiness failed." << endl;
203 TString avgclouds = CheckGraph(g) ? Form("%5.1f", g->GetMean(2)) : "NULL";
204 TString rmsclouds = CheckGraph(g) ? Form("%5.1f", g->GetRMS(2)) : "NULL";
205
206 g = (TGraph*)arr.FindObjectInCanvas("TempSky", "TGraph", "MHWeather");
207 //if (!g)
208 // cout << "WARNING - Reading of TempSky failed." << endl;
209 TString avgsky = CheckGraph(g) ? Form("%5.1f", g->GetMean(2)+200) : "NULL";
210
211
212 g = (TGraph*)arr.FindObjectInCanvas("NumStars", "TGraph", "MHPointing");
213 //if (!g)
214 // cout << "WARNING - Reading of NumStars failed." << endl;
215 TString numstarsmed = CheckGraph(g) ? Form("%5.1f", TMath::Median(g->GetN(), g->GetY())) : "NULL";
216 TString numstarsrms = CheckGraph(g) ? Form("%5.1f", g->GetRMS(2)) : "NULL";
217
218 g = (TGraph*)arr.FindObjectInCanvas("NumStarsCor", "TGraph", "MHPointing");
219 //if (!g)
220 // cout << "WARNING - Reading of NumStarsCor failed." << endl;
221 TString numcorsmed = CheckGraph(g) ? Form("%5.1f", TMath::Median(g->GetN(), g->GetY())) : "NULL";
222 TString numcorsrms = CheckGraph(g) ? Form("%5.1f", g->GetRMS(2)) : "NULL";
223
224 g = (TGraph*)arr.FindObjectInCanvas("Brightness", "TGraph", "MHPointing");
225 //if (!g)
226 // cout << "WARNING - Reading of SkyBrightness failed." << endl;
227 TString skybrightnessmed = CheckGraph(g) ? Form("%5.1f", TMath::Median(g->GetN(), g->GetY())) : "NULL";
228 TString skybrightnessrms = CheckGraph(g) ? Form("%5.1f", g->GetRMS(2)) : "NULL";
229
230 MSequence seq;
231 if (seq.Read("sequence[0-9]{8}[.]txt|MSequence")<=0)
232 {
233 cout << "ERROR - Could not find sequence in file: " << fname << endl;
234 return 2;
235 }
236 if (!seq.IsValid())
237 {
238 cout << "ERROR - Sequence read from file inavlid: " << fname << endl;
239 return 2;
240 }
241
242 cout << "Sequence M" << seq.GetTelescope() << ":" << seq.GetSequence() << endl;
243 cout << " Inhomogeneity " << inhomogen << endl;
244 cout << " PSF [mm] " << PSF << endl;
245 cout << " Island Quality " << islands << endl;
246 cout << " Ratio [%] " << ratio << endl;
247 cout << " Muon Number " << num << endl;
248 cout << " Eff. OnTime [s] " << effontime << endl;
249 cout << " Tot. OnTime [s] " << totontime << endl;
250 cout << " Rel. OnTime [%] " << relontime << endl;
251 cout << " Muon Rate [Hz] " << muonrate << endl;
252 cout << " # of Events (w/o sparks) " << numevents << endl;
253 cout << " # of Sparks " << numsparks << endl;
254 cout << " Rate after ImgCl [Hz] " << datarate << endl;
255 cout << " Rate of sparks [Hz] " << sparkrate << endl;
256 cout << " Minimum DC current [nA] " << mindc << endl;
257 cout << " Median DC current [nA] " << meddc << endl;
258 cout << " Maximum DC current [nA] " << maxdc << endl;
259 cout << " Maximum Humidity [%] " << maxhum << endl;
260 cout << " Average Humidity [%] " << avghum << endl;
261 cout << " Average WindSpeed [km/h] " << avgwind << endl;
262 cout << " Average Temp [" << UTF8::kDeg << "C] " << avgtemp << endl;
263 cout << " Average Sky Temp [K] " << avgsky << endl;
264 cout << " Cloundiness [%] " << avgclouds << " +/- " << rmsclouds << endl;
265 cout << " Number of Stars " << numstarsmed << " +/- " << numstarsrms << endl;
266 cout << " Number of cor. Stars " << numcorsmed << " +/- " << numcorsrms << endl;
267 cout << " Skybrightness " << skybrightnessmed << " +/- " << skybrightnessrms << endl;
268
269 TString vars = Form(" fSequenceFirst=%d, "
270 " fTelescopeNumber=%d, "
271 " fMeanNumberIslands=%s,"
272 " fRatio=%s,"
273 " fMuonNumber=%d,"
274 " fEffOnTime=%s,"
275 " fTotOnTime=%s,"
276 " fMuonRate=%s,"
277 " fPSF=%s,"
278 " fDataRate=%s,"
279 " fSparkRate=%s,"
280 " fMinCurrents=%s, "
281 " fMedCurrents=%s, "
282 " fMaxCurrents=%s, "
283 " fMaxHumidity=%s,"
284 " fAvgHumidity=%s, "
285 " fAvgTemperature=%s,"
286 " fAvgWindSpeed=%s,"
287 " fAvgTempSky=%s,"
288 " fAvgCloudiness=%s, "
289 " fRmsCloudiness=%s, "
290 " fNumStarsMed=%s,"
291 " fNumStarsRMS=%s,"
292 " fNumStarsCorMed=%s,"
293 " fNumStarsCorRMS=%s,"
294 " fBrightnessMed=%s,"
295 " fBrightnessRMS=%s,"
296 " fInhomogeneity=%s ",
297 seq.GetSequence(), seq.GetTelescope(),
298 islands.Data(), ratio.Data(),
299 num, effontime.Data(), totontime.Data(),
300 muonrate.Data(), PSF.Data(),
301 datarate.Data(), sparkrate.Data(),
302 mindc.Data(), meddc.Data(), maxdc.Data(),
303 maxhum.Data(), avghum.Data(),
304 avgtemp.Data(), avgwind.Data(),
305 avgsky.Data(), avgclouds.Data(), rmsclouds.Data(),
306 numstarsmed.Data(), numstarsrms.Data(),
307 numcorsmed.Data(), numcorsrms.Data(),
308 skybrightnessmed.Data(), skybrightnessrms.Data(),
309 inhomogen.Data());
310
311 TString where = Form("fSequenceFirst=%d AND fTelescopeNumber=%d",
312 seq.GetSequence(), seq.GetTelescope());
313
314 if (!serv.InsertUpdate("Star", vars, where))
315 return 2;
316
317 cout << endl;
318
319 h = (TH1*)arr.FindObjectInCanvas("Rate", "TH2D", "Rate");
320 if (!h)
321 return 1;
322
323 h->ResetBit(TH1::kCanRebin);
324
325 Int_t itrig = h->GetYaxis()->FindBin("Trig");
326 Int_t isum = h->GetYaxis()->FindBin("Sum");
327 Int_t inull = h->GetYaxis()->FindBin("0");
328
329 for (int i=0; i<h->GetNbinsX(); i++)
330 {
331 Int_t id = atoi(h->GetXaxis()->GetBinLabel(i+1));
332
333 Int_t run = seq.GetSequence()<1000000 ? id : id/1000 + 1000000;
334 Int_t file = seq.GetSequence()<1000000 ? 0 : id%1000;
335
336 const char *rtrig = itrig<0 ? "NULL" : Form("%8.1f", h->GetBinContent(i+1, itrig));
337 const char *rsum = isum <0 ? "NULL" : Form("%8.1f", h->GetBinContent(i+1, isum));
338 const char *rnull = inull<0 ? "NULL" : Form("%8.1f", h->GetBinContent(i+1, inull));
339
340 cout << " M" << seq.GetTelescope() << ":" << run << "/" << file << " " <<rtrig << " " << rsum << " " << rnull << endl;
341
342 TString vars = Form(" fTelescopeNumber=%d, "
343 " fRunNumber=%d, "
344 " fFileNumber=%d, "
345 " fRateCleanedTrig=%s,"
346 " fRateCleanedSum=%s,"
347 " fRateCleanedNull=%s",
348 seq.GetTelescope(), run, file,
349 rtrig, rsum, rnull);
350
351 TString where = Form("fTelescopeNumber=%d AND fRunNumber=%d AND fFileNumber=%d",
352 seq.GetTelescope(), run, file);
353
354 if (!serv.InsertUpdate("RunDataCheck", vars, where))
355 return 2;
356 }
357
358 return 1;
359}
360
361int fillstar(TString fname, Bool_t dummy=kTRUE)
362{
363 MSQLMagic serv("sql.rc");
364 if (!serv.IsConnected())
365 {
366 cout << "ERROR - Connection to database failed." << endl;
367 return 0;
368 }
369
370 cout << "fillstar" << endl;
371 cout << "--------" << endl;
372 cout << endl;
373 cout << "Connected to " << serv.GetName() << endl;
374 cout << "File: " << fname << endl;
375 cout << endl;
376
377 serv.SetIsDummy(dummy);
378
379 return Process(serv, fname);
380}
Note: See TracBrowser for help on using the repository browser.