source: trunk/MagicSoft/Mars/datacenter/macros/fillcamera.C@ 8737

Last change on this file since 8737 was 8606, checked in by hoehne, 18 years ago
*** empty log message ***
File size: 14.9 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, 06/2007 <mailto:tbretz@astro.uni-wuerzburg.de>
19! Author(s): Daniela Dorner, 11/2005 <mailto:dorner@astro.uni-wuerzburg.de>
20! Author(s): Daniel Hoehne, 06/2007 <mailto:hoehne@astro.uni-wuerzburg.de>
21!
22! Copyright: MAGIC Software Development, 2000-2007
23!
24!
25\* ======================================================================== */
26
27/////////////////////////////////////////////////////////////////////////////
28//
29// fillcamera.C
30// ============
31//
32// This macro is used to read the camera-output files.
33//
34// Make sure, that database and password are corretly set in a resource
35// file called sql.rc and the resource file is found.
36//
37// Returns 0 in case of failure and 1 in case of success.
38//
39//
40/////////////////////////////////////////////////////////////////////////////
41#include <iostream>
42#include <iomanip>
43
44#include <TEnv.h>
45#include <TRegexp.h>
46
47#include <TFile.h>
48#include <TTree.h>
49#include <TBranch.h>
50
51#include <TH1.h>
52
53#include <TSQLResult.h>
54#include <TSQLRow.h>
55
56#include "MSQLServer.h"
57#include "MSQLMagic.h"
58
59#include "MMcRunHeader.hxx"
60#include "MMcConfigRunHeader.h"
61#include "MMcCorsikaRunHeader.h"
62#include "MMcFadcHeader.hxx"
63#include "MMcEvtBasic.h"
64#include "MGeomCamMagic.h"
65#include "MRawRunHeader.h"
66
67#include <math.h>
68#include <MBinning.h>
69
70using namespace std;
71
72int Process(MSQLMagic &serv, TString fname, Bool_t dummy)
73{
74 TFile file(fname, "READ");
75 if (!file.IsOpen())
76 {
77 cout << "ERROR - Could not find file " << fname << endl;
78 return 2;
79 }
80
81 //
82 // Get tree RunHeaders from file
83 //
84 TTree *tree = dynamic_cast<TTree*>(file.Get("RunHeaders"));
85 if (!tree)
86 {
87 cout << "ERROR - Tree RunHeaders not found in file " << fname << endl;
88 return 2;
89 }
90
91 //
92 // Get branch MMcCorsikaRunHeader from tree RunHeaders
93 //
94 TBranch *b1 = tree->GetBranch("MMcCorsikaRunHeader.");
95 if (!b1)
96 {
97 cout << "ERROR - Branch MMcCorsikaRunHeader. not found in file " << fname << endl;
98 return 2;
99 }
100
101 MMcCorsikaRunHeader *runheader1 = 0;
102 b1->SetAddress(&runheader1);
103
104 //
105 // Get branch MMcConfigRunHeader from tree RunHeaders
106 //
107 TBranch *b2 = tree->GetBranch("MMcConfigRunHeader.");
108 if (!b2)
109 {
110 cout << "ERROR - Branch MMcConfigRunHeader. not found in file " << fname << endl;
111 return 2;
112 }
113
114 MMcConfigRunHeader *runheader2 = 0;
115 b2->SetAddress(&runheader2);
116
117 //
118 // Get branch MMcRunHeader from tree RunHeaders
119 //
120 TBranch *b3 = tree->GetBranch("MMcRunHeader.");
121 if (!b3)
122 {
123 cout << "ERROR - Branch MMcRunHeader. not found in file " << fname << endl;
124 return 2;
125 }
126
127 MMcRunHeader *runheader3 = 0;
128 b3->SetAddress(&runheader3);
129
130 //
131 // Get branch MMcFadcRunHeader from tree RunHeaders
132 //
133 TBranch *b4 = tree->GetBranch("MMcFadcHeader.");
134 if (!b4)
135 {
136 cout << "ERROR - Branch MMcFadcHeader. not found in file " << fname << endl;
137 return 2;
138 }
139
140 MMcFadcHeader *fadcheader = 0;
141 b4->SetAddress(&fadcheader);
142
143 //
144 // Get branch MRawRunHearder from tree RunHeaders
145 //
146 TBranch *b5 = tree->GetBranch("MRawRunHeader.");
147 if (!b5)
148 {
149 cout << "ERROR - Branch MRawRunHeader. not found in file " << fname << endl;
150 return 2;
151 }
152
153 MRawRunHeader *rawheader = 0;
154 b5->SetAddress(&rawheader);
155
156 tree->GetEvent(0);
157
158 //
159 // Get tree Events from file
160 //
161 TTree *tree2 = dynamic_cast<TTree*>(file.Get("Events"));
162 if (!tree2)
163 {
164 cout << "ERROR - Tree Events not found in file " << fname << endl;
165 return 2;
166 }
167
168 //
169 // Get branch MMcEvtBasic from tree Events
170 //
171 TBranch *b6 = tree2->GetBranch("MMcEvtBasic.");
172 if (!b6)
173 {
174 cout << "ERROR - Branch MMcEvtBasic. not found in file " << fname << endl;
175 return 2;
176 }
177
178 MMcEvtBasic *evtbasic = 0;
179 b6->SetAddress(&evtbasic);
180
181 tree2->GetEvent(0);
182
183 Double_t misptx = runheader2->GetMissPointingX();
184 Double_t mispty = runheader2->GetMissPointingY();
185 Double_t pointspreadx = runheader2->GetPointSpreadX();
186 Double_t tmin = runheader3->GetShowerThetaMin();
187 Double_t tmax = runheader3->GetShowerThetaMax();
188 UInt_t reflvers = runheader3->GetReflVersion();
189 UInt_t camvers = runheader3->GetCamVersion();
190 UInt_t numsimshow = runheader3->GetNumSimulatedShowers();
191 UInt_t numevents = tree2->GetEntries();
192
193 TString elow = Form("%5.1f", runheader1->GetELowLim());
194 TString eupp = Form("%5.1f", runheader1->GetEUppLim());
195 TString slope = Form("%5.1f", runheader1->GetSlopeSpec());
196 TString wobble = Form("%5.0f", runheader1->GetWobbleMode());
197 TString corsika1 = Form("%5.0f", runheader1->GetCorsikaVersion());
198 TString psf = Form("%5.1f", runheader2->GetPointSpread());
199 TString psfx = Form("%5.2f", pointspreadx);
200 TString psfy = Form("%5.2f", runheader2->GetPointSpreadY());
201 TString psfadd = Form("%5.2f", TMath::Hypot(runheader2->GetPointSpreadX(), runheader2->GetPointSpread()));
202 TString mirrfrac = Form("%5.2f", runheader2->GetMirrorFraction());
203 TString misx = Form("%5.2f", misptx);
204 TString misy = Form("%5.2f", mispty);
205// TString reflector = Form("%5.0f", runheader3->GetReflVersion());
206 TString reflector = Form("%5.0i", reflvers);
207 TString camera = Form("%5.0i", camvers);
208// TString camera = Form("%5.0f", runheader3->GetCamVersion());
209 TString imax = Form("%5.1f", runheader3->GetImpactMax());
210 TString numphe = Form("%5.1f", runheader3->GetNumPheFromDNSB());
211 TString pmin = Form("%5.1f", runheader3->GetShowerPhiMin());
212 TString pmax = Form("%5.1f", runheader3->GetShowerPhiMax());
213 TString numss = Form("%7.0i", numsimshow);
214 TString thetamin = Form("%5.1f", tmin);
215 TString thetamax = Form("%5.1f", tmax);
216 TString ped = Form("%5.1f", fadcheader->GetPedestal(1));
217 TString low2high = Form("%5.1f", fadcheader->GetLow2HighGain());
218 TString amplfadc = Form("%5.1f", fadcheader->GetAmplitud());
219 TString amplfadco = Form("%5.1f", fadcheader->GetAmplitudOuter());
220 TString enoise = Form("%5.1f", fadcheader->GetElecNoise(1));
221 TString dnoise = Form("%5.1f", fadcheader->GetDigitalNoise(1));
222
223 TH1I h("myhist", "", 1, -0.5, 0.5);
224 tree2->Draw("MRawEvtData.GetNumPixels()>>myhist", "", "goff");
225 h.SetDirectory(0);
226
227 TString numtrig = Form("%7.0i", TMath::Nint(h.GetBinContent(2)));
228 TString numevt = Form("%7.0i", numevents);
229 TString partid = Form("%5.0f", evtbasic->GetPartId());
230 TString partname = evtbasic->GetParticleName();
231
232 // Bestimmung von fakewobble aus file
233 // Kombination aus Wobble(0,1) und MissPoint
234 TString wobblemod="Wobble";
235
236 // Fake Wobble: als ObservationModeKEY 4 einfügen?
237 if (wobblemod != 1)
238 wobblemod = misptx == 0 && mispty == 0 ? "On" : "Fake Wobble";
239
240 Float_t pointsfx=TMath::Floor(pointspreadx*10);
241 pointsfx=TMath::Nint(pointsfx);
242 TString pointsfuncx=Form("%5.1f",pointsfx);
243
244 Float_t zBin=TMath::Nint((1-((TMath::Cos(tmin*TMath::Pi()/180)+TMath::Cos(tmax*TMath::Pi()/180))/2))*100);
245 zBin=TMath::Nint(zBin);
246
247 // folgende Werte werden aus dem Pfadnamen gewonnen
248 // RunNumber
249 TRegexp reg2("_[0-9][0-9][0-9][0-9][0-9][0-9][0-9][0-9]_");
250 TString Run = fname(reg2);
251 Int_t RunNum = atoi(Run.Data()+1);
252 if (RunNum < 1 || RunNum > 99999)
253 {
254 cout << "ERROR - RunNumber wrong value" << endl;
255 return 2;
256 }
257
258 // PointSpreadFunction
259 TRegexp reg4("/[12][0-9]/");
260 TString pointsf = fname(reg4);
261 Int_t Point = atoi(pointsf.Data()+1);
262 if (Point < 0 || Point > 99)
263 {
264 cout << "ERROR - PointSpreadFunction wrong value" << endl;
265 return 2;
266 }
267
268 // zbin
269 TRegexp reg1("/19[0-9][0-9]/");
270 TString zbin = fname(reg1);
271 Int_t ZBin = atoi(zbin.Data()+3);
272 if (ZBin < 0 || ZBin > 99)
273 {
274 cout << "ERROR - zbin wrong value" << endl;
275 return 2;
276 }
277
278 // WobbleMode
279 TRegexp reg3("/0[0-9]/"); /*extrahiert '/0X' mit X= 1-8 */
280 TString wm = fname(reg3); /* weist WM den extrahierten Wert zu */
281 Int_t wob = atoi(wm.Data()+1); /* schneidet fuehrenden '/' ab */
282
283
284 // Umrechnung von WobbleModus Bezeichnung in 'wobble', 'nowobble', 'fakewobble'
285 TString wobmode; // dieser Wert wird in 'MCDB' Tabelle 'WobbleMode' eingetragen
286 switch (wob)
287 {
288 case 1:
289 case 3:
290 wobmode = "Wobble";
291 break;
292
293 case 2:
294 case 4:
295 case 7:
296 case 8:
297 wobmode = "On";
298 break;
299
300 case 5:
301 case 6:
302 wobmode = "Fake Wobble";
303 break;
304
305 default:
306 cout << "ERROR - ObservationMode wrong value" << endl;
307 return 2;
308 }
309
310
311 cout << endl;
312 cout << endl;
313
314 cout << "--- From File ---" << endl;
315
316 cout << endl;
317 cout << "Energy Range " << elow << " < E < " << eupp << endl;
318 cout << "SpectralIndex " << slope << endl;
319 cout << "ObservationMode " << wobblemod << endl;
320 cout << "CorsikaVer " << corsika1 << endl;
321 cout << "ReflVer " << reflector << endl;
322 cout << "CamVer " << camera << endl;
323 cout << "ParticleId " << partid << endl;
324 cout << "ParticleName " << partname << endl;
325 cout << "PointSpread " << psf << endl;
326 cout << "PointSpreadXY " << psfx << " /" << psfy << endl;
327 cout << "AdditionPSF " << psfadd << endl;
328 cout << "MispointingXY " << misx << " /" << misy <<endl;
329 cout << "NumSimShowers " << numss << endl;
330 cout << "ImpactMax " << imax << endl;
331 cout << "NumEvents " << numevt << endl;
332 cout << "NumTriggers " << numtrig << endl;
333 cout << "NumPheFromDNSB " << numphe << endl;
334 cout << "Pedestal " << ped << endl;
335 cout << "Low2HighGain " << low2high << endl;
336 cout << "AmplitudeFADC " << amplfadc << endl;
337 cout << "AmplFADCOuter " << amplfadco << endl;
338 cout << "ElecNoise " << enoise << endl;
339 cout << "DigiNoise " << dnoise << endl;
340 cout << "PhiMin " << pmin << endl;
341 cout << "PhiMax " << pmax << endl;
342 cout << "ThetaMin " << thetamin << endl;
343 cout << "ThetaMax " << thetamax << endl;
344 cout << "Zenith range " << tmin << " to " << tmax << endl;
345 cout << "MirrorFraction " << mirrfrac << endl;
346 cout << endl;
347
348 cout << endl;
349 cout << "--- key's from mcdb tables ---" << endl;
350 cout << endl;
351
352 Int_t corsikakey = serv.QueryKeyOfName("CorsikaVersion", corsika1.Data());
353 Int_t reflectorkey = serv.QueryKeyOfName("ReflectorVersion", reflector.Data());
354 Int_t camerakey = serv.QueryKeyOfName("CameraVersion", camera.Data());
355 Int_t observationkey = serv.QueryKeyOfName("ObservationMode", wobblemod.Data());
356 Int_t particlekey = serv.QueryKeyOfName("MCParticle", partname.Data());
357
358 cout << "corsikakey: " << corsikakey << endl;
359 cout << "reflectorkey: " << reflectorkey << endl;
360 cout << "camerakey: " << camerakey << endl;
361 cout << "observationkey: " << observationkey << endl;
362 cout << "particlekey: " << particlekey << endl;
363
364 cout << endl;
365 cout << endl;
366 cout << "--- From File ---" << endl;
367 cout << endl;
368 cout << "WobbleMode " << wobblemod << endl;
369 cout << "PSF " << pointsfuncx << endl;
370 cout << "zBin " << zBin << endl;
371 cout << endl;
372 cout << "--- From FileName ---" << endl;
373 cout << endl;
374 cout << "WobbleMode " << wobmode << endl;
375 cout << "RunNum " << RunNum << endl;
376 cout << "PSF " << Point << endl;
377 cout << "ZBin " << ZBin << endl;
378 cout << "WobbleMode(dir) " << wob << endl;
379 cout << endl;
380
381
382 if (wobblemod!=wobmode)
383 {
384 cout << "Error, WobbleMode in file and filename are not the same" << endl;
385 return 2;
386 }
387 if (pointsfx!=Point)
388 {
389 cout << "Error, PSF in file and filename are not the same" << endl;
390 return 2;
391 }
392 if (zBin!=ZBin)
393 {
394 cout << "Error, ZBin in file and filename are not the same" << endl;
395 return 2;
396 }
397
398 TString vars =
399 Form("fELowLim=%s, fEUppLim=%s, fSlopeSpec=%s, "
400 "fImpactMax=%s, fNumSimulatedShowers=%d, fNumEvents=%d, "
401 "fNumPheFromDNSB=%s, fZBin=%d, fThetaMin=%s, "
402 "fThetaMax=%s, fPointSpread=%s, fPointSpreadX=%s, "
403 "fPointSpreadY=%s, fPedesMean=%s, fLow2HighGain=%s, ",
404 elow.Data(), eupp.Data(), slope.Data(), imax.Data(),
405 numsimshow, numevents, numphe.Data(),
406 zBin, thetamin.Data(), thetamax.Data(), psf.Data(),
407 psfx.Data(), psfy.Data(), ped.Data(), low2high.Data());
408 vars +=
409 Form("fAmplFadc=%s, fAmplFadcOuter=%s, fElectricNoise=%s, "
410 "fDigitalNoise=%s, fRunNumber=%d, fMisspointingX=%s, "
411 "fMissPointingY=%s, fCorsikaVersionKEY =%d, "
412 "fReflectorVersionKEY=%d, fCameraVersionKEY=%d, "
413 "fObservationModeKEY=%d, fMCParticleKEY=%d, "
414 "fSequenceFirst=0, fNumTriggers=%s ",
415 amplfadc.Data(), amplfadco.Data(), enoise.Data(),
416 dnoise.Data(), RunNum, misx.Data(), misy.Data(), corsikakey,
417 reflectorkey, camerakey, observationkey, particlekey, numtrig.Data());
418
419 // Comming soon (-1: dummy, 0: failed, 1: succeeded)
420 // return serv.InsertUpdate("MCRunData", "fRunNumber", RunData()+1, vars);
421
422 if (!serv.ExistStr("fRunNumber", "MCRunData", Run.Data()+1))
423 if (!serv.Insert("MCRunData", vars)==kFALSE)
424 return 2;
425 else
426 if (!serv.Update("MCRunData", vars, Form("fRunNumber=%d", Run))==kFALSE)
427 return 2;
428
429 return 1;
430}
431
432int fillcamera(TString fname, Bool_t dummy=kTRUE)
433{
434 TEnv env("sql.rc");
435
436 MSQLMagic serv(env);
437 if (!serv.IsConnected())
438 {
439 cout << "ERROR - Connection to database failed." << endl;
440 return 0;
441 }
442
443 serv.SetIsDummy(dummy);
444
445 cout << endl;
446 cout << "fillcamera" << endl;
447 cout << "----------" << endl;
448 cout << endl;
449 cout << "Connected to " << serv.GetName() << endl;
450 cout << "File: " << fname << endl;
451 cout << endl;
452
453 return Process(serv, fname, dummy);
454}
Note: See TracBrowser for help on using the repository browser.