source: trunk/MagicSoft/Mars/macros/calibrate_data.C@ 3838

Last change on this file since 3838 was 3832, checked in by gaug, 21 years ago
*** empty log message ***
File size: 12.0 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
59 MStatusDisplay *display = new MStatusDisplay;
60 display->SetUpdateTime(3000);
61 display->Resize(850,700);
62
63 gStyle->SetOptStat(1111);
64 gStyle->SetOptFit();
65
66 /************************************/
67 /* FIRST LOOP: PEDESTAL COMPUTATION */
68 /************************************/
69
70 MParList plist1;
71 MTaskList tlist1;
72 plist1.AddToList(&tlist1);
73
74 // containers
75 MPedestalCam pedcam;
76 MBadPixelsCam badcam;
77 //
78 // for excluding pixels from the beginning:
79 //
80 // badcam.AsciiRead("badpixels.dat");
81
82
83 plist1.AddToList(&pedcam);
84 plist1.AddToList(&badcam);
85
86 //tasks
87 MReadMarsFile read("Events");
88 read.DisableAutoScheme();
89 static_cast<MRead&>(read).AddFiles(pruns);
90
91 MGeomApply geomapl;
92 MPedCalcPedRun pedcalc;
93 MGeomCamMagic geomcam;
94
95 tlist1.AddToList(&read);
96 tlist1.AddToList(&geomapl);
97 tlist1.AddToList(&pedcalc);
98
99 // Create and execute the event looper
100 MEvtLoop pedloop;
101 pedloop.SetParList(&plist1);
102 pedloop.SetDisplay(display);
103
104 cout << "*************************" << endl;
105 cout << "** COMPUTING PEDESTALS **" << endl;
106 cout << "*************************" << endl;
107
108 if (!pedloop.Eventloop())
109 return;
110
111 tlist1.PrintStatistics();
112
113 //
114 // Now the short version:
115 //
116 //
117 // Now setup the new tasks for the calibration:
118 // ---------------------------------------------------
119 //
120 MCalibrationQECam qecam;
121 MJCalibration calloop;
122 calloop.SetInput(&cruns);
123 // calloop.SetFullDisplay();
124 //
125 // Use as signal extractor MExtractSignal2:
126 //
127 calloop.SetExtractorLevel(2);
128 //
129 // Set the corr. cams:
130 //
131 calloop.SetQECam(qecam);
132 calloop.SetBadPixels(badcam);
133 //
134 // The next two commands are for the display:
135 //
136 calloop.SetDisplay(display);
137
138 //
139 // Apply rel. time calibration:
140 //
141 calloop.SetRelTimeCalibration();
142 //
143 // Use as arrival time extractor MArrivalTimeCalc2:
144 //
145 calloop.SetArrivalTimeLevel(2);
146
147 //
148 // Do the event-loop:
149 //
150 cout << "***************************" << endl;
151 cout << "** COMPUTING CALIBRATION **" << endl;
152 cout << "***************************" << endl;
153
154 if (!calloop.Process(pedcam))
155 return;
156
157 badcam.Print();
158
159 MBadPixelsCam &badbad = calloop.GetBadPixels();
160 MCalibrationChargeCam &calcam = calloop.GetCalibrationCam();
161 MCalibrationRelTimeCam &timecam = calloop.GetRelTimeCam();
162
163 badbad.Print();
164
165 /************************************************************************/
166 /* THIRD LOOP: DATA CALIBRATION INTO PHOTONS */
167 /************************************************************************/
168
169 // Create an empty Parameter List and an empty Task List
170 MParList plist3;
171 MTaskList tlist3;
172 plist3.AddToList(&tlist3);
173
174 // containers
175 MCerPhotEvt photevt;
176 MPedPhotCam pedphotcam;
177 MSrcPosCam srccam;
178 MRawRunHeader runhead;
179 MExtractedSignalCam sigcam;
180
181 plist3.AddToList(&geomcam );
182 plist3.AddToList(&pedcam );
183 plist3.AddToList(&calcam );
184 plist3.AddToList(&qecam );
185 plist3.AddToList(&badbad );
186 plist3.AddToList(&timecam );
187 plist3.AddToList(&sigcam );
188 plist3.AddToList(&photevt);
189 plist3.AddToList(&pedphotcam);
190 plist3.AddToList(&srccam);
191 plist3.AddToList(&runhead);
192
193 //tasks
194 MReadMarsFile read3("Events");
195 read3.DisableAutoScheme();
196 static_cast<MRead&>(read3).AddFiles(druns);
197
198 MExtractSignal2 sigcalc;
199 MArrivalTimeCalc2 timecalc;
200 MCalibrateData photcalc;
201 photcalc.SetCalibrationMode(MCalibrateData::kFfactor); // !!! was only MCalibrate
202 // MPedPhotCalc pedphotcalc; // already done by MCalibrate Data
203 // MCerPhotCalc cerphotcalc; // already done by MCalibrate Data
204
205 tlist3.AddToList(&read3);
206 tlist3.AddToList(&geomapl);
207 tlist3.AddToList(&sigcalc);
208 tlist3.AddToList(&timecalc);
209 // tlist3.AddToList(&cerphotcalc); // already done by MCalibrate Data
210 tlist3.AddToList(&photcalc);
211 // tlist3.AddToList(&pedphotcalc); // already done by MCalibrate Data
212
213 MWriteRootFile write(resname);
214
215 write.AddContainer("MGeomCam" , "RunHeaders");
216 write.AddContainer("MRawRunHeader" , "RunHeaders");
217 write.AddContainer("MSrcPosCam" , "RunHeaders");
218 write.AddContainer("MCalibrationChargeCam" , "RunHeaders");
219 write.AddContainer("MCalibrationQECam" , "RunHeaders");
220 // write.AddContainer("MPedPhotCam","RunHeaders"); // Attention, was in Events - Tree!!
221 write.AddContainer("MPedestalCam" , "RunHeaders");
222 write.AddContainer("MCalibrationRelTimeCam", "RunHeaders");
223
224 write.AddContainer("MCerPhotEvt" , "Events");
225 write.AddContainer("MRawEvtHeader" , "Events");
226 write.AddContainer("MBadPixelsCam" , "Events");
227 write.AddContainer("MPedPhotCam" , "Events");
228
229 tlist3.AddToList(&write);
230
231 // Create and execute eventloop
232 MEvtLoop evtloop3;
233 evtloop3.SetParList(&plist3);
234
235 cout << "*************************************************************" << endl;
236 cout << "*** COMPUTING DATA USING EXTRACTED SIGNAL (IN PHOTONS) ***" << endl;
237 cout << "*************************************************************" << endl;
238
239 if (!evtloop3.Eventloop())
240 return;
241 tlist3.PrintStatistics();
242
243}
244
245void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
246{
247
248 c.cd(i);
249 gPad->SetBorderMode(0);
250 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
251 // obj1->AddNotify(evt);
252
253 c.cd(i+j);
254 gPad->SetBorderMode(0);
255 obj1->Draw();
256 ((MHCamera*)obj1)->SetPrettyPalette();
257
258 if (fit != 0)
259 {
260 c.cd(i+2*j);
261 gPad->SetBorderMode(0);
262 TH1D *obj2 = (TH1D*)obj1->Projection(obj1.GetName());
263
264// obj2->Sumw2();
265 obj2->Draw();
266 obj2->SetBit(kCanDelete);
267
268 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
269 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
270 const Double_t integ = obj2->Integral("width")/2.5066283;
271 const Double_t mean = obj2->GetMean();
272 const Double_t rms = obj2->GetRMS();
273 const Double_t width = max-min;
274
275 if (rms == 0. || width == 0. )
276 return;
277
278 switch (fit)
279 {
280 case 1:
281 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
282 sgaus->SetBit(kCanDelete);
283 sgaus->SetParNames("Area","#mu","#sigma");
284 sgaus->SetParameters(integ/rms,mean,rms);
285 sgaus->SetParLimits(0,0.,integ);
286 sgaus->SetParLimits(1,min,max);
287 sgaus->SetParLimits(2,0,width/1.5);
288 obj2->Fit("sgaus","QLR");
289 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
290 break;
291
292 case 2:
293 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
294 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
295 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
296 dgaus->SetBit(kCanDelete);
297 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
298 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
299 integ/width/2.,(max+mean)/2.,width/4.);
300 // The left-sided Gauss
301 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
302 dgaus->SetParLimits(1,min+(width/10.),mean);
303 dgaus->SetParLimits(2,0,width/2.);
304 // The right-sided Gauss
305 dgaus->SetParLimits(3,0,integ);
306 dgaus->SetParLimits(4,mean,max-(width/10.));
307 dgaus->SetParLimits(5,0,width/2.);
308 obj2->Fit("dgaus","QLRM");
309 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
310 break;
311
312 case 3:
313 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
314 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
315 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
316 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
317 tgaus->SetBit(kCanDelete);
318 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
319 "A_{2}","#mu_{2}","#sigma_{2}",
320 "A_{3}","#mu_{3}","#sigma_{3}");
321 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
322 integ/width/3.,(max+mean)/2.,width/4.,
323 integ/width/3.,mean,width/2.);
324 // The left-sided Gauss
325 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
326 tgaus->SetParLimits(1,min+(width/10.),mean);
327 tgaus->SetParLimits(2,width/15.,width/2.);
328 // The right-sided Gauss
329 tgaus->SetParLimits(3,0.,integ);
330 tgaus->SetParLimits(4,mean,max-(width/10.));
331 tgaus->SetParLimits(5,width/15.,width/2.);
332 // The Gauss describing the outliers
333 tgaus->SetParLimits(6,0.,integ);
334 tgaus->SetParLimits(7,min,max);
335 tgaus->SetParLimits(8,width/4.,width/1.5);
336 obj2->Fit("tgaus","QLRM");
337 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
338 break;
339 case 4:
340 obj2->Fit("pol0","Q");
341 obj2->GetFunction("pol0")->SetLineColor(kYellow);
342 break;
343 case 9:
344 break;
345 default:
346 obj2->Fit("gaus","Q");
347 obj2->GetFunction("gaus")->SetLineColor(kYellow);
348 break;
349 }
350
351 TArrayI s0(3);
352 s0[0] = 6;
353 s0[1] = 1;
354 s0[2] = 2;
355
356 TArrayI s1(3);
357 s1[0] = 3;
358 s1[1] = 4;
359 s1[2] = 5;
360
361 TArrayI inner(1);
362 inner[0] = 0;
363
364 TArrayI outer(1);
365 outer[0] = 1;
366
367 // Just to get the right (maximum) binning
368 TH1D *half[4];
369 half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
370 half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
371 half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
372 half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");
373
374 for (int i=0; i<4; i++)
375 {
376 half[i]->SetLineColor(kRed+i);
377 half[i]->SetDirectory(0);
378 half[i]->SetBit(kCanDelete);
379 half[i]->Draw("same");
380 }
381
382 gPad->Modified();
383 gPad->Update();
384
385 }
386}
387
388
Note: See TracBrowser for help on using the repository browser.