source: tags/Mars-V0.8.4/macros/calibrate_data.C

Last change on this file was 3762, checked in by gaug, 21 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): Hendrik Bartko, 03/2004 <mailto:hbartko@mppmu.mpg.de>
19! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26const TString defpath = "/mnt/Data/rootdata/CrabNebula/2004_01_27/";
27const TString defrout = "output_test.root";
28
29const Int_t defpedr [] = {12461};
30const Int_t defcalr [] = {12526,12527,12528,12529};
31const Int_t defdatar[] = {12517,12518};
32
33void calibrate_data(const TString inpath=defpath,
34 const Int_t psize=1, const Int_t pedruns[]=defpedr,
35 const Int_t csize=4, const Int_t calruns[]=defcalr,
36 const Int_t dsize=2, const Int_t dataruns[]=defdatar,
37 const TString resname=defrout)
38
39{
40
41 MRunIter pruns;
42 MRunIter cruns;
43 MRunIter druns;
44
45 for (Int_t i=0;i<psize;i++) {
46 cout << "Adding pedestal run: " << pedruns[i] << endl;
47 pruns.AddRun(pedruns[i],inpath);
48 }
49 for (Int_t i=0;i<csize;i++) {
50 cout << "Adding calibration run: " << calruns[i] << endl;
51 cruns.AddRun(calruns[i],inpath);
52 }
53 for (Int_t i=0;i<dsize;i++) {
54 cout << "Adding data run: " << dataruns[i] << endl;
55 druns.AddRun(dataruns[i],inpath);
56 }
57
58 MCalibrationCam::PulserColor_t color;
59
60 if (calruns[0] < 20000)
61 color = MCalibrationCam::kCT1;
62 else
63 color = FindColor((MDirIter*)&cruns);
64
65 if (color == MCalibrationCam::kNONE)
66 {
67 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
68
69 while (1)
70 {
71 timer.TurnOn();
72 TString input = Getline("Could not find the correct colour: Type 'q' to exit, "
73 "green, blue, uv or ct1 to go on: ");
74 timer.TurnOff();
75
76 if (input=="q\n")
77 return ;
78
79 if (input=="green")
80 color = MCalibrationCam::kGREEN;
81 if (input=="blue")
82 color = MCalibrationCam::kBLUE;
83 if (input=="uv")
84 color = MCalibrationCam::kUV;
85 if (input=="ct1")
86 color = MCalibrationCam::kCT1;
87 }
88 }
89
90 MStatusDisplay *display = new MStatusDisplay;
91 display->SetUpdateTime(3000);
92 display->Resize(850,700);
93
94 gStyle->SetOptStat(1111);
95 gStyle->SetOptFit();
96
97 /************************************/
98 /* FIRST LOOP: PEDESTAL COMPUTATION */
99 /************************************/
100
101 MParList plist1;
102 MTaskList tlist1;
103 plist1.AddToList(&tlist1);
104
105 // containers
106 MPedestalCam pedcam;
107 MBadPixelsCam badcam;
108 //
109 // for excluding pixels from the beginning:
110 //
111 // badcam.AsciiRead("badpixels.dat");
112
113 plist1.AddToList(&pedcam);
114 plist1.AddToList(&badcam);
115
116 //tasks
117 MReadMarsFile read("Events");
118 read.DisableAutoScheme();
119 static_cast<MRead&>(read).AddFiles(pruns);
120
121 MGeomApply geomapl;
122 MPedCalcPedRun pedcalc;
123 MGeomCamMagic geomcam;
124
125 tlist1.AddToList(&read);
126 tlist1.AddToList(&geomapl);
127 tlist1.AddToList(&pedcalc);
128
129 // Create and execute the event looper
130 MEvtLoop pedloop;
131 pedloop.SetParList(&plist1);
132 pedloop.SetDisplay(display);
133
134 cout << "*************************" << endl;
135 cout << "** COMPUTING PEDESTALS **" << endl;
136 cout << "*************************" << endl;
137
138 if (!pedloop.Eventloop())
139 return;
140
141 tlist1.PrintStatistics();
142
143 //
144 // Now the short version:
145 //
146 //
147 // Now setup the new tasks for the calibration:
148 // ---------------------------------------------------
149 //
150 MCalibrationQECam qecam;
151 MJCalibration calloop;
152 calloop.SetColor(color);
153 calloop.SetInput(&cruns);
154 // calloop.SetFullDisplay();
155 //
156 // Use as signal extractor MExtractSignal2:
157 //
158 calloop.SetExtractorLevel(2);
159 //
160 // Set the corr. cams:
161 //
162 calloop.SetQECam(qecam);
163 calloop.SetBadPixels(badcam);
164 //
165 // The next two commands are for the display:
166 //
167 calloop.SetDisplay(display);
168
169 //
170 // Apply rel. time calibration:
171 //
172 calloop.SetRelTimeCalibration();
173 //
174 // Use as arrival time extractor MArrivalTimeCalc2:
175 //
176 calloop.SetArrivalTimeLevel(2);
177
178 //
179 // Do the event-loop:
180 //
181 cout << "***************************" << endl;
182 cout << "** COMPUTING CALIBRATION **" << endl;
183 cout << "***************************" << endl;
184
185 if (!calloop.Process(pedcam))
186 return;
187
188 MBadPixelsCam &badbad = calloop.GetBadPixels();
189 MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
190 MCalibrationRelTimeCam &timecam = calloop.GetRelTimeCam();
191
192 /************************************************************************/
193 /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */
194 /************************************************************************/
195
196 // Create an empty Parameter List and an empty Task List
197 MParList plist3;
198 MTaskList tlist3;
199 plist3.AddToList(&tlist3);
200
201 // containers
202 MCerPhotEvt photevt;
203 MPedPhotCam pedphotcam;
204 MSrcPosCam srccam;
205 MRawRunHeader runhead;
206 MExtractedSignalCam sigcam;
207
208 plist3.AddToList(&geomcam );
209 plist3.AddToList(&pedcam );
210 plist3.AddToList(&calcam );
211 plist3.AddToList(&qecam );
212 plist3.AddToList(&badbad );
213 plist3.AddToList(&timecam );
214 plist3.AddToList(&sigcam );
215 plist3.AddToList(&photevt);
216 plist3.AddToList(&pedphotcam);
217 plist3.AddToList(&srccam);
218 plist3.AddToList(&runhead);
219
220 //tasks
221 MReadMarsFile read3("Events");
222 read3.DisableAutoScheme();
223 static_cast<MRead&>(read3).AddFiles(druns);
224
225 MExtractSignal2 sigcalc;
226 MArrivalTimeCalc2 timecalc;
227 MCalibrateData photcalc;
228 photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate
229 // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data
230 // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data
231
232 tlist3.AddToList(&read3);
233 tlist3.AddToList(&geomapl);
234 tlist3.AddToList(&sigcalc);
235 tlist3.AddToList(&timecalc);
236 // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data
237 tlist3.AddToList(&photcalc);
238 // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data
239
240 MWriteRootFile write(resname);
241
242 write.AddContainer("MGeomCam" , "RunHeaders");
243 write.AddContainer("MRawRunHeader" , "RunHeaders");
244 write.AddContainer("MSrcPosCam" , "RunHeaders");
245 write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
246 write.AddContainer("MCalibrationQECam" , "RunHeaders");
247 // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!!
248 write.AddContainer("MPedestalCam" , "RunHeaders");
249 write.AddContainer("MCalibrationRelTimeCam", "RunHeaders");
250
251 write.AddContainer("MCerPhotEvt" , "Events");
252 write.AddContainer("MRawEvtHeader" , "Events");
253 write.AddContainer("MBadPixelsCam" , "Events");
254 write.AddContainer("MPedPhotCam" , "Events");
255
256 tlist3.AddToList(&write);
257
258 // Create and execute eventloop
259 MEvtLoop evtloop3;
260 evtloop3.SetParList(&plist3);
261
262 cout << "*************************************************************" << endl;
263 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
264 cout << "*************************************************************" << endl;
265
266 if (!evtloop3.Eventloop())
267 return;
268 tlist3.PrintStatistics();
269
270}
271
272void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
273{
274
275 c.cd(i);
276 gPad->SetBorderMode(0);
277 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
278 // obj1->AddNotify(evt);
279
280 c.cd(i+j);
281 gPad->SetBorderMode(0);
282 obj1->Draw();
283 ((MHCamera*)obj1)->SetPrettyPalette();
284
285 if (fit != 0)
286 {
287 c.cd(i+2*j);
288 gPad->SetBorderMode(0);
289 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
290
291// obj2->Sumw2();
292 obj2->Draw();
293 obj2->SetBit(kCanDelete);
294
295 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
296 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
297 const Double_t integ = obj2->Integral("width")/2.5066283;
298 const Double_t mean = obj2->GetMean();
299 const Double_t rms = obj2->GetRMS();
300 const Double_t width = max-min;
301
302 if (rms == 0. || width == 0. )
303 return;
304
305 switch (fit)
306 {
307 case 1:
308 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
309 sgaus->SetBit(kCanDelete);
310 sgaus->SetParNames("Area","#mu","#sigma");
311 sgaus->SetParameters(integ/rms,mean,rms);
312 sgaus->SetParLimits(0,0.,integ);
313 sgaus->SetParLimits(1,min,max);
314 sgaus->SetParLimits(2,0,width/1.5);
315 obj2->Fit("sgaus","QLR");
316 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
317 break;
318
319 case 2:
320 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
321 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
322 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
323 dgaus->SetBit(kCanDelete);
324 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
325 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
326 integ/width/2.,(max+mean)/2.,width/4.);
327 // The left-sided Gauss
328 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
329 dgaus->SetParLimits(1,min+(width/10.),mean);
330 dgaus->SetParLimits(2,0,width/2.);
331 // The right-sided Gauss
332 dgaus->SetParLimits(3,0,integ);
333 dgaus->SetParLimits(4,mean,max-(width/10.));
334 dgaus->SetParLimits(5,0,width/2.);
335 obj2->Fit("dgaus","QLRM");
336 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
337 break;
338
339 case 3:
340 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
341 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
342 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
343 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
344 tgaus->SetBit(kCanDelete);
345 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
346 "A_{2}","#mu_{2}","#sigma_{2}",
347 "A_{3}","#mu_{3}","#sigma_{3}");
348 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
349 integ/width/3.,(max+mean)/2.,width/4.,
350 integ/width/3.,mean,width/2.);
351 // The left-sided Gauss
352 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
353 tgaus->SetParLimits(1,min+(width/10.),mean);
354 tgaus->SetParLimits(2,width/15.,width/2.);
355 // The right-sided Gauss
356 tgaus->SetParLimits(3,0.,integ);
357 tgaus->SetParLimits(4,mean,max-(width/10.));
358 tgaus->SetParLimits(5,width/15.,width/2.);
359 // The Gauss describing the outliers
360 tgaus->SetParLimits(6,0.,integ);
361 tgaus->SetParLimits(7,min,max);
362 tgaus->SetParLimits(8,width/4.,width/1.5);
363 obj2->Fit("tgaus","QLRM");
364 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
365 break;
366 case 4:
367 obj2->Fit("pol0","Q");
368 obj2->GetFunction("pol0")->SetLineColor(kYellow);
369 break;
370 case 9:
371 break;
372 default:
373 obj2->Fit("gaus","Q");
374 obj2->GetFunction("gaus")->SetLineColor(kYellow);
375 break;
376 }
377
378 TArrayI s0(3);
379 s0[0] = 6;
380 s0[1] = 1;
381 s0[2] = 2;
382
383 TArrayI s1(3);
384 s1[0] = 3;
385 s1[1] = 4;
386 s1[2] = 5;
387
388 TArrayI inner(1);
389 inner[0] = 0;
390
391 TArrayI outer(1);
392 outer[0] = 1;
393
394 // Just to get the right (maximum) binning
395 TH1D *half[4];
396 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
397 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
398 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
399 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
400
401 for (int i=0; i<4; i++)
402 {
403 half[i]->SetLineColor(kRed+i);
404 half[i]->SetDirectory(0);
405 half[i]->SetBit(kCanDelete);
406 half[i]->Draw("same");
407 }
408
409 gPad->Modified();
410 gPad->Update();
411
412 }
413}
414
415
416MCalibrationCam::PulserColor_t FindColor(MDirIter* run)
417{
418
419 MCalibrationCam::PulserColor_t col = MCalibrationCam::kNONE;
420
421 TString filenames;
422
423 while (!(filenames=run->Next()).IsNull())
424 {
425
426 filenames.ToLower();
427
428 if (filenames.Contains("green"))
429 if (col == MCalibrationCam::kNONE)
430 {
431 cout << "Found colour: Green in " << filenames << endl;
432 col = MCalibrationCam::kGREEN;
433 }
434 else if (col != MCalibrationCam::kGREEN)
435 {
436 cout << "Different colour found in " << filenames << "... abort" << endl;
437 return MCalibrationCam::kNONE;
438 }
439
440 if (filenames.Contains("blue"))
441 if (col == MCalibrationCam::kNONE)
442 {
443 cout << "Found colour: Blue in " << filenames << endl;
444 col = MCalibrationCam::kBLUE;
445 }
446 else if (col != MCalibrationCam::kBLUE)
447 {
448 cout << "Different colour found in " << filenames << "... abort" << endl;
449 return MCalibrationCam::kNONE;
450 }
451
452 if (filenames.Contains("uv"))
453 if (col == MCalibrationCam::kNONE)
454 {
455 cout << "Found colour: Uv in " << filenames << endl;
456 col = MCalibrationCam::kUV;
457 }
458 else if (col != MCalibrationCam::kUV)
459 {
460 cout << "Different colour found in " << filenames << "... abort" << endl;
461 return MCalibrationCam::kNONE;
462 }
463
464 if (filenames.Contains("ct1"))
465 if (col == MCalibrationCam::kNONE)
466 {
467 cout << "Found colour: Ct1 in " << filenames << endl;
468 col = MCalibrationCam::kCT1;
469 }
470 else if (col != MCalibrationCam::kCT1)
471 {
472 cout << "Different colour found in " << filenames << "... abort" << endl;
473 return MCalibrationCam::kNONE;
474 }
475
476 }
477
478
479
480 if (col == MCalibrationCam::kNONE)
481 cout << "No colour found in filenames of runs: " << ((MRunIter*)run)->GetRunsAsString()
482 << "... abort" << endl;
483
484 return col;
485}
Note: See TracBrowser for help on using the repository browser.