source: tags/Mars-V0.9.2/macros/calibration.C

Last change on this file was 6394, checked in by gaug, 20 years ago
*** empty log message ***
File size: 20.5 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): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// calibration.C
27//
28// Needs as arguments the run number of a calibration file ("*_C_*.root") and
29// the run number of the corresponding pedestal file ("*_P_*.root").
30//
31// The TString inpath has to be set correctly.
32//
33// The macro searches for the pulser colour which corresponds to the calibration
34// run number. If the run number is smaller than 20000, pulser colour "CT1"
35// is assumed, otherwise, it searches for the strings "green", "blue", "uv" or
36// "ct1" in the filenames. If no colour or multiple colours are found, the
37// execution is aborted.
38//
39// The container MBadPixelsCam is created and followed during the execution of the
40// rest of the macro.
41//
42// A first loop over the pedestal file is performed using the class MJPedestal
43//
44// The container MCalibrationQECam is created and followed during the execution of the
45// rest of the macro.
46//
47// A loop over the calibration files is performed using the class MJCalibration.
48// The results are displayed using the MJCalibration::SetNormalDisplay() facility,
49// but other displays can easily be uncommented.
50// The call to MJCalibration skips the relative time calibration, which can be
51// uncommented as well.
52//
53// Last, a third loop is performed over the calibration file again in order to
54// "calibrate" it and test the resulting outcome.
55//
56//////////////////////////////////////////////////////////////////////////////////////////
57#include <TStyle.h>
58#include <TObject.h>
59#include <TObjectTable.h>
60#include <TCanvas.h>
61#include <TPad.h>
62#include <TH1.h>
63#include <TPaveStats.h>
64#include <TApplication.h>
65#include <TClass.h>
66
67#include "getExtractor.C"
68
69#include "MJPedestal.h"
70#include "MJCalibration.h"
71#include "MJCalibTest.h"
72#include "MJCalibrateSignal.h"
73#include "MRunIter.h"
74#include "MStatusDisplay.h"
75#include "MCalibrationQECamMagic.h"
76#include "MCalibrationQECam.h"
77#include "MBadPixelsCam.h"
78#include "MArgs.h"
79#include "MArray.h"
80#include "MParContainer.h"
81#include "MGeomCamMagic.h"
82
83using namespace std;
84
85static TString outpath = ".";
86//static TString inpath = "/home/rootdata/Calib/2004_07_06";
87//static TString inpath = "/home/rootdata/Calib/2004_11_11";
88//static TString inpath = "/home/rootdata/Calib/2004_12_18";
89static TString inpath = "/home/rootdata/Calib/2004_09_21";
90//static TString inpath = "/home/rootdata/Calib/1999_01_01";
91static TString badfile;
92//static TString badfile = "badpixels_only0_388_559.dat";
93//
94// the default pedestal run for the calibration
95//
96//static const Int_t pedrun = 43915;
97static const Int_t pedrun = 38995;
98//
99// the default start calibration run
100//
101static const Int_t calrun1 = 38997;
102//
103// the default last calibration run (if 0, only one run is taken, otherwise consecutive runs
104// between calrun1 and calrun2)
105//
106static const Int_t calrun2 = 0;
107//
108// the default start data run
109//
110static const Int_t datrun1 = 39000;
111//
112// the default last calibration run (if 0, only one run is taken, otherwise consecutive runs
113// between calrun1 and calrun2)
114//
115static const Int_t datrun2 = 0;
116//
117// A switch to output debugging information about Objects use
118//
119static Bool_t debug = kFALSE;
120//
121// A switch to use the Blind Pixel
122//
123static Bool_t blindpix = kTRUE;
124//
125// A switch to use the PIN Diode
126//
127static Bool_t pindiode = kFALSE;
128//
129// Tell if you want to use the display:
130//
131static Bool_t useDisplay = kTRUE;
132//
133// Tell if you want to test the result afterwards
134//
135static Bool_t useTest = kFALSE;
136//
137// Tell if you want to test the intensity calibration
138//
139static Bool_t useIntensity = kFALSE;
140//
141// Tell if you to calibrated interlaced calibration events
142//
143static Bool_t useInterlaced = kTRUE;
144//
145// Tell if you want to store and read the F0 and F1- files
146//
147static Bool_t useStorage = kTRUE;
148//
149// Tell which extractor you want to use. The flags are counted according
150// to the extractor-TDAS
151//
152static Int_t extractorflag = 35;
153//
154Int_t calibration(const UInt_t extflag=extractorflag, const Int_t prun=pedrun,
155 const Int_t crun1=calrun1, const Int_t crun2=calrun2,
156 const Int_t drun1=datrun1, const Int_t drun2=datrun2)
157{
158
159 cout << extractorflag << endl;
160
161 MExtractor *extractor = getExtractor(extflag==35 ? 32 : extflag);
162
163 if (!extractor)
164 return 99;
165
166 extractor->SetName(Form("%s_Run_%05d",extractor->GetName(),prun));
167 const Bool_t timeandcharge = extractor->InheritsFrom("MExtractTimeAndCharge");
168
169 gStyle->SetOptStat(1111);
170 gStyle->SetOptFit();
171 gStyle->SetTitleSize(0.35,"u");
172 gStyle->SetTitleFontSize(0.9);
173 gStyle->SetTitleH(0.12);
174 gStyle->SetTitleW(0.95);
175 gStyle->SetLineWidth(1);
176
177 if (debug)
178 TObject::SetObjectStat(kTRUE);
179
180 MRunIter pruns;
181 MRunIter cruns;
182
183 pruns.AddRun(prun,inpath);
184
185 if (crun1!=0)
186 if (crun2==0)
187 cruns.AddRun(crun1,inpath);
188 else
189 cruns.AddRuns(crun1,crun2,inpath);
190
191 MStatusDisplay *display = NULL;
192
193 if (useDisplay)
194 {
195 display = new MStatusDisplay;
196 display->SetUpdateTime(3000);
197 display->Resize(850,700);
198 }
199
200 /*****************************************/
201 /* FIRST LOOP: PURE PEDESTAL COMPUTATION */
202 /*****************************************/
203 //
204 // Hand over to the jobs a QE Cam with Cornings initialized
205 //
206 MGeomCamMagic geomcam;
207 MCalibrationQECamMagic qecam;
208 MBadPixelsCam badcam;
209 badcam.InitSize(geomcam.GetNumPixels());
210 //
211 // If you want to exclude pixels from the beginning, read
212 // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
213 //
214 if (!badfile.IsNull())
215 {
216 ifstream f(badfile.Data());
217 badcam.AsciiRead((istream&)f);
218 f.close();
219 }
220
221 MJPedestal pedloop1;
222 pedloop1.SetNoStorage();
223 pedloop1.SetEnvDebug(debug);
224 pedloop1.SetExtractor(extractor);
225 pedloop1.SetExtractionFundamental();
226 pedloop1.SetInput(&pruns);
227 pedloop1.SetPathOut(outpath.Data());
228 if (useDisplay)
229 {
230 pedloop1.SetDisplay(display);
231 pedloop1.SetDataCheckDisplay();
232 }
233 pedloop1.SetBadPixels(badcam);
234
235 if (!pedloop1.Process())
236 return 1;
237
238 /****************************************/
239 /* SECOND LOOP: CALIBRATION COMPUTATION */
240 /****************************************/
241
242 MJPedestal pedloop2;
243 MJCalibration calloop;
244
245 if (timeandcharge)
246 {
247 /***********************************************************/
248 /* NEEDED FOR SECOND LOOP: EXTRACTOR RESOLUTION COMPUTATION */
249 /***********************************************************/
250
251 pedloop2.SetUseData();
252 pedloop2.SetNoStorage();
253 pedloop2.SetEnvDebug(debug);
254 pedloop2.SetExtractor(extractor);
255 pedloop2.SetExtractionWithExtractorRndm();
256 pedloop2.SetPedestals(pedloop1.GetPedestalCam());
257 pedloop2.SetInput(&pruns);
258 pedloop2.SetPathOut(outpath.Data());
259 if (useDisplay)
260 {
261 pedloop2.SetDisplay(display);
262 pedloop2.SetDataCheckDisplay();
263 }
264 pedloop2.SetBadPixels(badcam);
265
266 if (!pedloop2.Process())
267 return 1;
268
269 calloop.SetExtractorCam(pedloop2.GetPedestalCam());
270 }
271
272 if (crun1 == 0)
273 return 0;
274
275 MPedestalCam &pedcam = pedloop1.GetPedestalCam();
276
277 if (debug)
278 calloop.SetDebug();
279 calloop.SetEnvDebug(debug);
280 if (useIntensity)
281 calloop.SetIntensity();
282 // calloop.SetHistsStorage();
283 calloop.SetNoStorage(!useStorage);
284 calloop.SetRelTimeCalibration(kTRUE);
285 //
286 // If you want to set a colour explicitely from outside (not recommanded!)
287 // calloop.SetColor(MCalibrationCam::kUV);
288 //
289 // If you want to run the data-check on RAW DATA!!!, choose:
290 // calloop.SetDataCheck();
291 //
292 // If you want to see the data-check plots only, choose:
293 calloop.SetDataCheckDisplay();
294 // calloop.SetNormalDisplay();
295 //
296 // For everything, you have ever dreamed of, choose:
297 // calloop.SetFullDisplay();
298
299 //
300 // If you want to calibrate the times as well, choose:
301 //
302 calloop.SetExtractor(extractor);
303 calloop.SetInput(&cruns);
304 calloop.SetPathOut(outpath.Data());
305 if (useDisplay)
306 calloop.SetDisplay(display);
307 calloop.SetUseBlindPixel(blindpix);
308 calloop.SetUsePINDiode(pindiode);
309 calloop.SetQECam(qecam);
310 calloop.SetBadPixels(badcam);
311
312 if (!calloop.Process(pedcam))
313 return 2;
314
315 //
316 // The next lines are the use the Print() function and have
317 // all the results as ascii-tables:
318 //
319 if (debug)
320 {
321 MCalibrationQECam &nqecam = calloop.GetQECam();
322 MBadPixelsCam &badbad = calloop.GetBadPixels();
323 MCalibrationChargeCam &chargecam = calloop.GetCalibrationCam();
324 chargecam.Print();
325 nqecam.Print();
326 badbad.Print();
327 }
328
329 /********************************************************************/
330 /* THIRD LOOP: APPLY CALIBRATION TO THE CALIBRATION FILES FOR TESTS */
331 /********************************************************************/
332
333 if (useTest)
334 {
335
336 MJCalibTest testloop;
337
338 // If you want to see the data-check plots only, choose:
339 testloop.SetDataCheckDisplay();
340
341 testloop.SetExtractor(extractor);
342 testloop.SetInput(&cruns);
343 testloop.SetPathOut(outpath);
344 if (useDisplay)
345 testloop.SetDisplay(display);
346 testloop.SetBadPixels(calloop.GetBadPixels());
347
348 if (!testloop.ProcessFile(pedloop1.GetPedestalCam()))
349 return 3;
350
351 }
352
353 if (drun1 == 0)
354 return 4;
355
356 MRunIter druns;
357
358 if (drun2==0)
359 druns.AddRun(drun1,inpath);
360 else
361 druns.AddRuns(drun1,drun2,inpath);
362
363
364 /********************************************************************/
365 /* FOURTH LOOP: APPLY CALIBRATION TO THE PEDESTAL FILES */
366 /********************************************************************/
367
368 if (extflag == 35)
369 {
370 delete extractor;
371 extractor = getExtractor(28);
372 }
373
374 MJCalibrateSignal calibloop;
375
376 // calibloop.SetExtractor(extractor);
377 calibloop.SetInputCal(&cruns);
378 calibloop.SetInput(&druns);
379 calibloop.SetPathIn(outpath);
380 if (useDisplay)
381 calibloop.SetDisplay(display);
382 // calibloop.SetBadPixels(calloop.GetBadPixels());
383 // calibloop.SetNoStorage(!useStorage);
384 calibloop.SetInterlaced(useInterlaced);
385
386 if (!calibloop.ProcessFile(pedloop1.GetPedestalCam(),
387 timeandcharge ? pedloop2.GetPedestalCam() : pedloop1.GetPedestalCam(),
388 timeandcharge ? pedloop2.GetPedestalCam() : pedloop1.GetPedestalCam()))
389 return 5;
390
391 if (debug)
392 TObject::SetObjectStat(kFALSE);
393
394 //
395 // Debugging at the end:
396 //
397 if (debug)
398 gObjectTable->Print();
399
400 return 0;
401}
402
403static void Usage()
404{
405 gLog << endl;
406 gLog << "Usage:" << endl;
407 gLog << endl;
408 gLog << " calibration [ped.run nr.] [first cal.run nr.] [last cal.run nr.] [first dat.run nr.] [last dat.run nr.]" << endl ;
409 gLog << endl;
410 gLog << " ped.run.nr: Run number of the pedestal file." << endl;
411 gLog << " first cal.run nr.: Run number of the first calibration file." << endl;
412 gLog << " last cal.run nr.: Run number of the last calibration file." << endl;
413 gLog << " first dat.run nr.: Run number of the first data file to be calibrated." << endl;
414 gLog << " last dat.run nr.: Run number of the last data file to be calibrated." << endl;
415 gLog << endl;
416 gLog << "All calibration runs between (first cal.run nr.) and (last cal.run nr.) will be used" << endl;
417 gLog << "If last.cal.run.nr is 0 (default), only one calibration run is taken" << endl;
418 gLog << endl;
419 gLog << "All data runs between (first dat.run nr.) and (last dat.run nr.) will be used" << endl;
420 gLog << "If last.dat.run.nr is 0 (default), only one data run is taken" << endl;
421 gLog << endl;
422 gLog << "Additional Options: " << endl;
423 gLog << " --extractor=# Choose one of the following possible extractors (integer)" << endl;
424 gLog << " (default: Nr. 33) " << endl;
425 gLog << endl;
426 gLog << " Nr. Extractor Parameters " << endl;
427 gLog << endl;
428 gLog << " MExtractFixedWindow: " << endl;
429 gLog << " with the following parameters: " << endl;
430 gLog << " 1: SetRange(4,7,6,9) " << endl;
431 gLog << " 2: SetRange(4,7,5,10)" << endl;
432 gLog << " 3: SetRange(3,8,5,10)" << endl;
433 gLog << " 4: SetRange(2,9,4,11)" << endl;
434 gLog << " 5: SetRange(0,13,4,13)" << endl;
435 gLog << " MExtractFixedWindowSpline: " << endl;
436 gLog << " 6: SetRange(3,7,5,9)" << endl;
437 gLog << " 7: SetRange(3,7,5,11)" << endl;
438 gLog << " 8: SetRange(3,9,5,11)" << endl;
439 gLog << " 9: SetRange(2,10,4,12)" << endl;
440 gLog << " 10: SetRange(0,14,4,14)" << endl;
441 gLog << " MExtractFixedWindowPeakSearch: " << endl;
442 gLog << " SetRange(0,18,2,14) and the following parameters:" << endl;
443 gLog << " 11: SetWindows(2,2,2)" << endl;
444 gLog << " SetOffsetFromWindow(0)" << endl;
445 gLog << " 12: SetWindows(4,4,2)" << endl;
446 gLog << " SetOffsetFromWindow(1)" << endl;
447 gLog << " 13: SetWindows(4,6,4)" << endl;
448 gLog << " SetOffsetFromWindow(0)" << endl;
449 gLog << " 14: SetWindows(6,6,4)" << endl;
450 gLog << " SetOffsetFromWindow(1)" << endl;
451 gLog << " 15: SetWindows(8,8,4)" << endl;
452 gLog << " SetOffsetFromWindow(1)" << endl;
453 gLog << " 16: SetWindows(14,10,4)" << endl;
454 gLog << " SetOffsetFromWindow(2)" << endl;
455 gLog << " MExtractTimeAndChargeSlidingWindow:" << endl;
456 gLog << " SetRange(0,18,2,14) and the following parameters:" << endl;
457 gLog << " 17: SetWindowSize(2,2)" << endl;
458 gLog << " 18: SetWindowSize(4,4)" << endl;
459 gLog << " 19: SetWindowSize(4,6)" << endl;
460 gLog << " 20: SetWindowSize(6,6)" << endl;
461 gLog << " 21: SetWindowSize(8,8)" << endl;
462 gLog << " 22: SetWindowSize(14,10)" << endl;
463 gLog << " MExtractTimeAndChargeSpline: " << endl;
464 gLog << " 23: SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) and:" << endl;
465 gLog << " SetRange(0,13,2,13) " << endl;
466 gLog << " 24: SetChargeType(MExtractTimeAndChargeSpline::kIntegral) and:" << endl;
467 gLog << " SetRiseTime(0.5); SetFallTime(0.5)" << endl;
468 gLog << " 25: SetRiseTime(0.5); SetFallTime(1.5)" << endl;
469 gLog << " 26: SetRiseTime(1.0); SetFallTime(3.0)" << endl;
470 gLog << " 27 SetRiseTime(1.5); SetFallTime(4.5)" << endl;
471 gLog << " MExtractTimeAndChargeDigitalFilter" << endl;
472 gLog << " SetRange(0,18,2,14) and the following parameters:" << endl;
473 gLog << " 28: SetNameWeightsFile('msignal/cosmics_weights.dat')" << endl;
474 gLog << " 29: SetNameWeightsFile('msignal/cosmics_weights4.dat')" << endl;
475 gLog << " 30: SetNameWeightsFile('msignal/cosmics_weights_logain6.dat')" << endl;
476 gLog << " 31: SetNameWeightsFile('msignal/cosmics_weights_logain4.dat')" << endl;
477 gLog << " 32: SetNameWeightsFile('msignal/calibration_weights_UV.dat')" << endl;
478 gLog << " 33: SetNameWeightsFile('msignal/calibration_weights_UV_logain.dat')" << endl;
479 gLog << " 34: SetNameWeightsFile('msignal/MC_weights.dat')" << endl;
480 gLog << " 35: SetNameWeightsFile('msignal/cosmics_weights.dat') (for cosmics) and" << endl;
481 gLog << " SetNameWeightsFile('msignal/calibration_weights_UV.dat') (for cal.)" << endl;
482 gLog << " MExtractTimeAndChargeDigitalFilterPeakSearch" << endl;
483 gLog << " 36: SetNameWeightsFile('msignal/calibration_weights_UV.dat')" << endl;
484 gLog << endl;
485
486 gLog << " --inpath=# Find the data in inpath" << endl;
487 gLog << " --outpath=# Write the output containers to outpath" << endl;
488 gLog << " --badfile=# Use the file # to exclude pixels from the beginning" << endl;
489 gLog << " --debug Use the TObjectTable for debugging " << endl;
490 gLog << " and write out the pixels as ascii tables" << endl;
491 gLog << " --useInterlaced Use the program to calibrate data with the interlaced cal. events" << endl;
492 gLog << " --useTest Use the class MJCalibTest to test the calibration on itself" << endl;
493 gLog << " --skipBlindPix Skip the blind pixel calibration" << endl;
494 gLog << " --skipPINDiode Skip the PIN Diode calibration" << endl;
495}
496
497
498int main(int argc, char **argv)
499{
500
501
502 MArray::Class()->IgnoreTObjectStreamer();
503 MParContainer::Class()->IgnoreTObjectStreamer();
504 //
505 // Evaluate arguments
506 //
507 MArgs arg(argc, argv);
508
509 if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
510 {
511 Usage();
512 return -1;
513 }
514
515 debug = arg.HasOnlyAndRemove("--debug") || arg.HasOnlyAndRemove("-d");
516 useTest = arg.HasOnlyAndRemove("--useTest") || arg.HasOnlyAndRemove("-t");
517 useInterlaced = arg.HasOnlyAndRemove("--useInterlaced") || arg.HasOnlyAndRemove("-i");
518 blindpix = !(arg.HasOnlyAndRemove("--skipBlindPix"));
519 pindiode = !(arg.HasOnlyAndRemove("--skipPINDiode"));
520
521 if (arg.HasOption("--extractor="))
522 extractorflag = arg.GetIntAndRemove("--extractor=");
523
524 if (arg.HasOption("--inpath="))
525 inpath = arg.GetStringAndRemove("--inpath=");
526
527 if (arg.HasOption("--outpath="))
528 outpath = arg.GetStringAndRemove("--outpath=");
529
530 if (arg.HasOption("--badfile="))
531 badfile = arg.GetStringAndRemove("--badfile=");
532
533 if (gSystem->AccessPathName(badfile,kFileExists))
534 {
535 gLog << "WARNING: the bad pixels file '" << badfile.Data() << "' doesn't exist." << endl;
536 badfile = "";
537 }
538
539 // check for the right usage of the program
540 //
541 if (arg.GetNumArguments()>6)
542 {
543 Usage();
544 return -1;
545 }
546
547 //
548 // Initialize Non-GUI (batch) mode
549 //
550 gROOT->SetBatch();
551
552 //
553 // Switch off the display
554 //
555 useDisplay = kTRUE;
556 //
557 // check for the arguments
558 //
559 Int_t pedr = 0;
560 Int_t calr1 = 0;
561 Int_t calr2 = 0;
562 Int_t datr1 = 0;
563 Int_t datr2 = 0;
564
565 const Int_t nargs = arg.GetNumArguments();
566
567 if (nargs>=5)
568 {
569 pedr = arg.GetArgumentInt(0);
570 calr1 = arg.GetArgumentInt(1);
571 calr2 = arg.GetArgumentInt(2);
572 datr1 = arg.GetArgumentInt(3);
573 datr2 = arg.GetArgumentInt(4);
574 return calibration(extractorflag,pedr,calr1,calr2,datr1,datr2);
575 }
576
577 if (nargs>=4)
578 {
579 pedr = arg.GetArgumentInt(0);
580 calr1 = arg.GetArgumentInt(1);
581 calr2 = arg.GetArgumentInt(2);
582 datr1 = arg.GetArgumentInt(3);
583 datr2 = arg.GetArgumentInt(4);
584 return calibration(extractorflag,pedr,calr1,calr2,datr1,0);
585 }
586
587 if (nargs>=3)
588 {
589 pedr = arg.GetArgumentInt(0);
590 calr1 = arg.GetArgumentInt(1);
591 calr2 = arg.GetArgumentInt(2);
592 return calibration(extractorflag,pedr,calr1,calr2,0);
593 }
594
595 if (nargs>=2)
596 {
597 pedr = arg.GetArgumentInt(0);
598 calr1 = arg.GetArgumentInt(1);
599 gLog << "PEDR: " << pedr << " CALR1: " << calr1 << " CALR2 " << calr2 << endl;
600 gLog << "inpath: " << inpath << endl;
601 gLog << "extractor: " << extractorflag << endl;
602 return calibration(extractorflag,pedr,calr1,0,0);
603 }
604
605 if (nargs>=1)
606 {
607 pedr = arg.GetArgumentInt(0);
608 return calibration(extractorflag,pedr,0);
609 }
610
611 return calibration(extractorflag,pedr,calr1,calr2);
612}
613
Note: See TracBrowser for help on using the repository browser.