source: trunk/MagicSoft/Mars/macros/pedphotcalc.C@ 3842

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