source: branches/Mars_use_drstimefiles/macros/pedphotcalc.C@ 18459

Last change on this file since 18459 was 4773, checked in by gaug, 20 years ago
*** empty log message ***
File size: 15.8 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! Author(s): Josep Flix, 01/2004 <mailto:jflix@ifae.es>
18! Javier Rico, 01/2004 <mailto:jrico@ifae.es>
19! Markus Gaug, 03/2004 <mailto:markus@ifae.es>
20!
21! (based on bootcampstandardanalysis.C by Javier López)
22!
23!
24! Copyright: MAGIC Software Development, 2000-2004
25!
26!
27\* ======================================================================== */
28/////////////////////////////////////////////////////////////////////////////
29//
30// pedphotcalc.C
31//
32// This macro is an example of the use of MPedPhotCalc. It computes
33// the pedestal mean and rms from pedestal files undergoing the
34// signal extraction + calibration chain, in units of photons. It
35// produces plots of the relevant computed quantities.
36//
37// Needs as arguments the run number of a pedestal file ("*_P_*.root"),
38// one of a calibration file ("*_C_*.root").
39// performs the pedestal calculation, the calibration
40/// constants calculation and the calibration of the pedestals.
41//
42// The TString inpath has to be set correctly.
43//
44// The macro searches for the pulser colour which corresponds to the calibration
45// run number. If the run number is smaller than 20000, pulser colour "CT1"
46// is assumed, otherwise, it searches for the strings "green", "blue", "uv" or
47// "ct1" in the filenames. If no colour or multiple colours are found, the
48// execution is aborted.
49//
50////////////////////////////////////////////////////////////////////////////////////
51#include "MParList.h"
52#include "MTaskList.h"
53#include "MJPedestal.h"
54#include "MJCalibration.h"
55#include "MPedestalCam.h"
56#include "MPedestalPix.h"
57#include "MReadMarsFile.h"
58#include "MGeomApply.h"
59#include "MGeomCamMagic.h"
60#include "MEvtLoop.h"
61#include "MCalibrationCam.h"
62#include "MCalibrationChargeCam.h"
63#include "MCalibrationChargePix.h"
64#include "MCalibrationQECam.h"
65#include "MCalibrationQEPix.h"
66#include "MExtractedSignalCam.h"
67#include "MExtractSlidingWindow.h"
68#include "MExtractFixedWindow.h"
69#include "MCerPhotEvt.h"
70#include "MCalibrateData.h"
71#include "MPedPhotCalc.h"
72#include "MPedPhotCam.h"
73#include "MPedPhotPix.h"
74#include "MHCamera.h"
75#include "MRunIter.h"
76#include "MDirIter.h"
77#include "MStatusDisplay.h"
78#include "MHCamera.h"
79
80#include "TTimer.h"
81#include "TString.h"
82#include "TCanvas.h"
83#include "TStyle.h"
84#include "TF1.h"
85#include "TLegend.h"
86
87#include <iostream.h>
88#include "Getline.h"
89
90const TString inpath = "/mnt/Data/rootdata/CrabNebula/2004_02_10/";
91const Int_t dpedrun = 14607;
92const Int_t dcalrun1 = 14608;
93const Int_t dcalrun2 = 0;
94const Bool_t usedisplay = kTRUE;
95
96
97
98void DrawProjection(MHCamera *obj1, Int_t fit)
99{
100 TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
101 obj2->SetDirectory(0);
102 obj2->Draw();
103 obj2->SetBit(kCanDelete);
104 gPad->SetLogy();
105
106 if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
107 {
108 TArrayI s0(3);
109 s0[0] = 6;
110 s0[1] = 1;
111 s0[2] = 2;
112
113 TArrayI s1(3);
114 s1[0] = 3;
115 s1[1] = 4;
116 s1[2] = 5;
117
118 TArrayI inner(1);
119 inner[0] = 0;
120
121 TArrayI outer(1);
122 outer[0] = 1;
123
124 // Just to get the right (maximum) binning
125 TH1D *half[4];
126 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
127 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
128 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
129 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
130
131 for (int i=0; i<4; i++)
132 {
133 half[i]->SetLineColor(kRed+i);
134 half[i]->SetDirectory(0);
135 half[i]->SetBit(kCanDelete);
136 half[i]->Draw("same");
137 }
138 }
139
140 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
141 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
142 const Double_t integ = obj2->Integral("width")/2.5066283;
143 const Double_t mean = obj2->GetMean();
144 const Double_t rms = obj2->GetRMS();
145 const Double_t width = max-min;
146
147 if (rms==0 || width==0)
148 return;
149
150 const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
151 "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])");
152
153 TF1 *f=0;
154 switch (fit)
155 {
156 // FIXME: MAYBE add function to TH1?
157 case 0:
158 f = new TF1("sgaus", "gaus(0)", min, max);
159 f->SetLineColor(kYellow);
160 f->SetBit(kCanDelete);
161 f->SetParNames("Area", "#mu", "#sigma");
162 f->SetParameters(integ/rms, mean, rms);
163 f->SetParLimits(0, 0, integ);
164 f->SetParLimits(1, min, max);
165 f->SetParLimits(2, 0, width/1.5);
166 obj2->Fit(f, "QLRM");
167 break;
168
169 case 1:
170 f = new TF1("dgaus", dgausformula, min, max);
171 f->SetLineColor(kYellow);
172 f->SetBit(kCanDelete);
173 f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
174 f->SetParameters(integ, (min+mean)/2, width/4,
175 integ/width/2, (max+mean)/2, width/4);
176 // The left-sided Gauss
177 f->SetParLimits(0, integ-1.5, integ+1.5);
178 f->SetParLimits(1, min+(width/10), mean);
179 f->SetParLimits(2, 0, width/2);
180 // The right-sided Gauss
181 f->SetParLimits(3, 0, integ);
182 f->SetParLimits(4, mean, max-(width/10));
183 f->SetParLimits(5, 0, width/2);
184 obj2->Fit(f, "QLRM");
185 break;
186
187 default:
188 obj2->Fit("gaus", "Q");
189 obj2->GetFunction("gaus")->SetLineColor(kYellow);
190 break;
191 }
192}
193
194void CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
195{
196 c.cd(x);
197 gPad->SetBorderMode(0);
198 MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");
199
200 c.cd(x+y);
201 gPad->SetBorderMode(0);
202 obj1->Draw();
203
204 c.cd(x+2*y);
205 gPad->SetBorderMode(0);
206 DrawProjection(obj1, fit);
207}
208
209void pedphotcalc(const Int_t prun=dpedrun, // pedestal file
210 const Int_t crun1=dcalrun1, const Int_t crun2=dcalrun2 // calibration file(s)
211 )
212{
213
214 MExtractFixedWindow extractor;
215
216 MRunIter pruns;
217 MRunIter cruns;
218
219 pruns.AddRun(prun,inpath);
220
221 if (crun2==0)
222 cruns.AddRun(crun1,inpath);
223 else
224 cruns.AddRuns(crun1,crun2,inpath);
225
226 //
227 // Now setup the tasks and tasklist for the pedestals:
228 // ---------------------------------------------------
229 //
230 MGeomCamMagic geomcam;
231 MGeomApply geomapl;
232 MStatusDisplay *display = NULL;
233
234 /************************************/
235 /* FIRST LOOP: PEDESTAL COMPUTATION */
236 /************************************/
237
238 MJPedestal pedloop;
239 pedloop.SetInput(&pruns);
240 if (usedisplay)
241 {
242 display = new MStatusDisplay;
243 display->SetUpdateTime(3000);
244 display->Resize(850,700);
245 display->SetBit(kCanDelete);
246 pedloop.SetDisplay(display);
247 }
248
249 cout << "*************************" << endl;
250 cout << "** COMPUTING PEDESTALS **" << endl;
251 cout << "*************************" << endl;
252
253 if (!pedloop.Process())
254 return;
255
256 MPedestalCam pedcam = pedloop.GetPedestalCam();
257
258 /****************************/
259 /* SECOND LOOP: CALIBRATION */
260 /****************************/
261
262 //
263 // Now setup the new tasks for the calibration:
264 // ---------------------------------------------------
265 //
266 MJCalibration calloop;
267 calloop.SetInput(&cruns);
268 //
269 // Use as signal extractor MExtractSignal:
270 //
271 calloop.SetExtractor(&extractor);
272 //
273 // The next two commands are for the display:
274 //
275 if (usedisplay)
276 {
277 calloop.SetDisplay(display);
278 calloop.SetDataCheck();
279 }
280
281 cout << "***********************************" << endl;
282 cout << "** COMPUTING CALIBRATION FACTORS **" << endl;
283 cout << "***********************************" << endl;
284
285 if (!calloop.Process(pedcam))
286 return;
287
288 MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
289 MCalibrationQECam &qecam = calloop.GetQECam();
290
291 /************************************************************************/
292 /* THIRD LOOP: PEDESTAL COMPUTATION USING EXTRACTED SIGNAL (IN PHOTONS) */
293 /************************************************************************/
294
295 // Create an empty Parameter List and an empty Task List
296 MParList plist3;
297 MTaskList tlist3;
298 plist3.AddToList(&tlist3);
299
300 // containers
301 MCerPhotEvt photcam;
302 MPedPhotCam pedphotcam;
303 MExtractedSignalCam extedsig;
304
305 plist3.AddToList(&geomcam);
306 plist3.AddToList(&pedcam );
307 plist3.AddToList(&calcam );
308 plist3.AddToList(&qecam );
309 plist3.AddToList(&photcam);
310 plist3.AddToList(&extedsig);
311 plist3.AddToList(&pedphotcam);
312
313 //tasks
314 MReadMarsFile read3("Events");
315 read3.DisableAutoScheme();
316 static_cast<MRead&>(read3).AddFiles(pruns);
317
318 MCalibrateData photcalc;
319 photcalc.SetCalibrationMode(MCalibrateData::kFfactor);
320 MPedPhotCalc pedphotcalc;
321
322 tlist3.AddToList(&read3);
323 tlist3.AddToList(&geomapl);
324 tlist3.AddToList(&extractor);
325 tlist3.AddToList(&photcalc);
326 tlist3.AddToList(&pedphotcalc);
327
328 // Create and execute eventloop
329 MEvtLoop evtloop3;
330 evtloop3.SetParList(&plist3);
331
332 cout << "*************************************************************" << endl;
333 cout << "** COMPUTING PEDESTALS USING EXTRACTED SIGNAL (IN PHOTONS) **" << endl;
334 cout << "*************************************************************" << endl;
335
336 if (!evtloop3.Eventloop())
337 return;
338 tlist3.PrintStatistics();
339
340 /**********************/
341 /* PRODUCE NICE PLOTS */
342 /**********************/
343
344 if (usedisplay)
345 {
346
347 MHCamera disp0(geomcam, "MPedPhotCam;ped", "Pedestals");
348 MHCamera disp1(geomcam, "MPedPhotCam;rms", "Sigma Pedestal");
349
350 disp0.SetCamContent(pedphotcam, 0);
351 disp0.SetCamError (pedphotcam, 1);
352
353 disp1.SetCamContent(pedphotcam, 1);
354
355 disp0.SetYTitle("Pedestal signal [photons]");
356 disp1.SetYTitle("Pedestal signal RMS [photons]");
357
358 //
359 // Display data
360 //
361 TCanvas &c1 = display->AddTab("PedPhotCam");
362 c1.Divide(2,3);
363
364 CamDraw(c1, 1, 2, disp0, 0);
365 CamDraw(c1, 2, 2, disp1, 1);
366
367 }
368
369
370#if 0
371 const UShort_t npix = 577;
372
373 // declare histograms
374 // pedestals
375 TH1F* pedhist = new TH1F("pedhist","Pedestal",npix,0,npix);
376 TH1F* pedrmshist = new TH1F("pedrmshist","Pedestal RMS",npix,0,npix);
377 TH1F* peddist = new TH1F("peddist","Pedestal",100,0,20);
378 TH1F* pedrmsdist = new TH1F("pedrmsdist","Pedestal RMS",100,0,15);
379
380 // calibration factors
381 TH1F* calhist = new TH1F("calhist","Calibration factors",npix,0,npix);
382 TH1F* caldist = new TH1F("caldist","Calibration factors",100,0,1);
383 TH1F* qehist = new TH1F("qehist", "Quantrum efficiencies",npix,0,npix);
384 TH1F* qedist = new TH1F("qedist", "Quantrum efficiencies",100,0,1);
385
386 // pedestal signals
387 TH1F* pedphothist = new TH1F("pedphothist","Pedestal Signal",npix,0,npix);
388 TH1F* pedphotrmshist = new TH1F("pedphotrmshist","Pedestal Signal RMS",npix,0,npix);
389 TH1F* pedphotdist = new TH1F("pedphotdist","Pedestal Signal",100,-0.4,0.4);
390 TH1F* pedphotrmsdist = new TH1F("pedphotrmsdist","Pedestal Signal RMS",100,0,25);
391
392 // fill histograms
393 for(int i=0;i<npix;i++)
394 {
395 MCalibrationChargePix &calpix = (MCalibrationChargePix&)calcam[i];
396 MCalibrationQEPix &qepix = (MCalibrationQEPix&) qecam[i];
397
398 const Float_t ped = pedcam[i].GetPedestal();
399 const Float_t pedrms = pedcam[i].GetPedestalRms();
400 const Float_t cal = calpix.GetMeanConvFADC2Phe();
401 const Float_t qe = qepix .GetQECascadesFFactor();
402 const Float_t pedphot = pedphotcam[i].GetMean();
403 const Float_t pedphotrms = pedphotcam[i].GetRms();
404
405 pedhist->Fill(i,ped);
406 peddist->Fill(ped);
407 pedrmshist->Fill(i,pedrms);
408 pedrmsdist->Fill(pedrms);
409
410 calhist->Fill(i,cal);
411 caldist->Fill(cal);
412 qehist->Fill(i,qe);
413 qedist->Fill(qe);
414
415 pedphothist->Fill(i,pedphot);
416 pedphotdist->Fill(pedphot);
417 pedphotrmshist->Fill(i,pedphotrms);
418 pedphotrmsdist->Fill(pedphotrms);
419 }
420
421 // Draw
422 gROOT->Reset();
423 gStyle->SetCanvasColor(0);
424 TCanvas* canvas = new TCanvas("canvas","pedphotcalc.C", 0, 100, 650, 800);
425 canvas->SetBorderMode(0); // Delete the Canvas' border line
426 canvas->cd();
427
428 canvas->Divide(2,5);
429
430 // draw pedestal histo
431 canvas->cd(1);
432 gPad->cd();
433 gPad->SetBorderMode(0);
434
435 pedhist->SetStats(kFALSE);
436 pedhist->GetXaxis()->SetTitle("Pixel SW Id");
437 pedhist->GetYaxis()->SetTitle("Pedestal (ADC counts)");
438 pedrmshist->SetStats(kFALSE);
439 pedrmshist->SetLineColor(2);
440 pedhist->Draw();
441 pedrmshist->Draw("same");
442
443 TLegend* leg1 = new TLegend(.14,.68,.39,.88);
444 leg1->SetHeader("");
445 leg1->AddEntry(pedhist,"Pedestal","L");
446 leg1->AddEntry(pedrmshist,"Pedestal RMS","L");
447 leg1->SetFillColor(0);
448 leg1->SetLineColor(0);
449 leg1->SetBorderSize(0);
450 leg1->Draw();
451
452 // draw pedestal distribution
453 canvas->cd(2);
454 gPad->cd();
455 gPad->SetBorderMode(0);
456 peddist->GetXaxis()->SetTitle("Pedestal (ADC counts)");
457 pedrmsdist->SetLineColor(2);
458 peddist->Draw();
459 pedrmsdist->Draw("same");
460
461 TLegend* leg2 = new TLegend(.14,.68,.39,.88);
462 leg2->SetHeader("");
463 leg2->AddEntry(pedhist,"Pedestal","L");
464 leg2->AddEntry(pedrmshist,"Pedestal RMS","L");
465 leg2->SetFillColor(0);
466 leg2->SetLineColor(0);
467 leg2->SetBorderSize(0);
468 leg2->Draw();
469
470 // draw calibration histo
471 canvas->cd(3);
472 gPad->cd();
473 gPad->SetBorderMode(0);
474 calhist->GetXaxis()->SetTitle("Pixel SW Id");
475 calhist->SetMaximum(1);
476 calhist->SetMinimum(0);
477 calhist->GetYaxis()->SetTitle("Calibration factor (phe/ADC count)");
478 calhist->SetStats(kFALSE);
479 calhist->Draw();
480
481 // draw calibration distribution
482 canvas->cd(4);
483 gPad->cd();
484 gPad->SetBorderMode(0);
485 caldist->GetXaxis()->SetTitle("Calibration factor (phe/ADC count)");
486 caldist->Draw();
487
488 // draw qe histo
489 canvas->cd(5);
490 gPad->cd();
491 gPad->SetBorderMode(0);
492 qehist->GetXaxis()->SetTitle("Pixel SW Id");
493 qehist->SetMaximum(1);
494 qehist->SetMinimum(0);
495 qehist->GetYaxis()->SetTitle("Quantum efficiency for cascades");
496 qehist->SetStats(kFALSE);
497 qehist->Draw();
498
499 // draw qe distribution
500 canvas->cd(6);
501 gPad->cd();
502 gPad->SetBorderMode(0);
503 qedist->GetXaxis()->SetTitle("Quantum efficiency for cascades");
504 qedist->Draw();
505
506 // draw pedestal signal histo
507 canvas->cd(7);
508 gPad->cd();
509 gPad->SetBorderMode(0);
510 pedphothist->GetXaxis()->SetTitle("Pixel SW Id");
511 pedphothist->GetYaxis()->SetTitle("Pedestal signal (photons)");
512 pedphothist->SetStats(kFALSE);
513 pedphothist->Draw();
514
515 // draw pedestal signal distribution
516 canvas->cd(8);
517 gPad->cd();
518 gPad->SetBorderMode(0);
519 pedphotdist->GetXaxis()->SetTitle("Pedestal signal (photons)");
520 pedphotdist->Draw();
521
522 // draw pedestal signal rms histo
523 canvas->cd(9);
524 gPad->cd();
525 gPad->SetBorderMode(0);
526 pedphotrmshist->GetXaxis()->SetTitle("Pixel SW Id");
527 pedphotrmshist->GetYaxis()->SetTitle("Pedestal signal rms (photons)");
528 pedphotrmshist->SetStats(kFALSE);
529 pedphotrmshist->Draw();
530
531 // draw pedestal signal rms distribution
532 canvas->cd(10);
533 gPad->cd();
534 gPad->SetBorderMode(0);
535 pedphotrmsdist->GetXaxis()->SetTitle("Pedestal signal rms (photons)");
536 pedphotrmsdist->Draw();
537
538 canvas->SaveAs("pedphotcalc.root");
539
540#endif
541}
Note: See TracBrowser for help on using the repository browser.