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

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