source: branches/Corsika7405Compatibility/mtemp/mmpi/macros/calibration_shape.C

Last change on this file was 5127, checked in by hbartko, 20 years ago
*** empty log message ***
File size: 10.4 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): Markus Gaug, 11/2003 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2003
21!
22!
23\* ======================================================================== */
24#include "MAGIC.h"
25
26
27const TString pedfile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12149_P_Cab-On_E.root";
28const TString calfile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12526_C_Cab-On_E.root";
29const TString datafile="/.magic/magicserv01/MAGIC/rootdata/2004_01_27/20040126_12517_D_Cab-On_E.root";
30
31Int_t pedr = 36038;
32Int_t calr = 36042;
33//Int_t datar = 12517;
34
35
36
37void calibration_shape(const TString inpath = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_06/", Int_t pedruns = pedr, Int_t calruns = calr, Int_t ipix=316)
38
39{
40
41
42 MRunIter pruns;
43 MRunIter cruns;
44
45 pruns.AddRun(pedruns,inpath);
46
47 cruns.AddRun(calruns,inpath);
48
49
50 gStyle->SetOptStat(1111);
51 gStyle->SetOptFit();
52
53 /************************************/
54 /* FIRST LOOP: PEDESTAL COMPUTATION */
55 /************************************/
56
57 MParList plist1;
58 MTaskList tlist1;
59 plist1.AddToList(&tlist1);
60
61 // containers
62 MPedestalCam pedcam;
63 MBadPixelsCam badcam;
64
65
66 plist1.AddToList(&pedcam);
67 plist1.AddToList(&badcam);
68
69 //tasks
70 MReadMarsFile read("Events");
71 read.DisableAutoScheme();
72 static_cast<MRead&>(read).AddFiles(pruns);
73
74 MGeomApply geomapl;
75 // MPedCalcPedRun pedcalc;
76 MPedCalcFromLoGain pedcalc_ped;
77 pedcalc_ped.SetMaxHiGainVar(20);
78 pedcalc_ped.SetRange(0, 11, 1, 14);
79 pedcalc_ped.SetWindowSize(12,14);
80 pedcalc_ped.SetPedestalUpdate(kFALSE);
81
82
83
84 MGeomCamMagic geomcam;
85
86 tlist1.AddToList(&read);
87 tlist1.AddToList(&geomapl);
88 tlist1.AddToList(&pedcalc_ped);
89
90 // Create and execute the event looper
91 MEvtLoop pedloop;
92 pedloop.SetParList(&plist1);
93
94 cout << "*************************" << endl;
95 cout << "** COMPUTING PEDESTALS **" << endl;
96 cout << "*************************" << endl;
97
98 if (!pedloop.Eventloop())
99 return;
100
101 tlist1.PrintStatistics();
102
103
104
105 //
106 // Create a empty Parameter List and an empty Task List
107 //
108 MParList plist2;
109 MTaskList tlist2;
110 plist2.AddToList(&tlist2);
111 plist2.AddToList(&pedcam);
112 plist2.AddToList(&badcam);
113
114
115 MReadMarsFile read2("Events");
116 read2.DisableAutoScheme();
117 static_cast<MRead&>(read2).AddFiles(cruns);
118
119 MGeomCamMagic geomcam;
120 MExtractedSignalCam sigcam;
121 MArrivalTimeCam timecam;
122 MRawEvtData evtdata;
123
124 //
125 // Get the previously created MPedestalCam into the new Parameter List
126 //
127 plist2.AddToList(&geomcam);
128 plist2.AddToList(&sigcam);
129 plist2.AddToList(&timecam);
130 plist2.AddToList(&evtdata);
131
132 //
133 // We saw that the signal jumps between slices,
134 // thus take the sliding window
135 //
136
137 MGeomApply geomapl;
138
139
140 Int_t WindowSize = 6; // size of the chosen integration window for signal and time extraction
141
142 MExtractFixedWindowPeakSearch extractor;
143 extractor.SetRange(0, 14, 0, 14);
144 extractor.SetWindows( WindowSize, WindowSize, 4);
145
146 MExtractTimeHighestIntegral timeext;
147 timeext.SetRange(0, 14, 0, 14);
148 //timeext.SetWindowSize(10,10);
149 timeext.SetWindowSize(WindowSize,WindowSize);
150
151
152 // timecalc->SetRange(8,14,8,14,6);
153
154 // filter against cosmics
155 MFCosmics cosmics;
156 MContinue cont(&cosmics);
157
158 tlist2.AddToList(&read2);
159 tlist2.AddToList(&geomapl);
160 tlist2.AddToList(&extractor);
161 tlist2.AddToList(&timeext);
162
163 //
164 // In case, you want to skip the cosmics rejection,
165 // uncomment the next line
166 //
167 tlist2.AddToList(&cont);
168 //
169 // In case, you want to skip the somewhat lengthy calculation
170 // of the arrival times using a spline, uncomment the next two lines
171 //
172
173
174 //
175 // Create and setup the eventloop
176 //
177 MEvtLoop evtloop;
178 evtloop.SetParList(&plist2);
179 // evtloop.SetDisplay(display);
180
181 cout << "***************************" << endl;
182 cout << "** COMPUTING CALIBRATION **" << endl;
183 cout << "***************************" << endl;
184
185 //
186 // Execute second analysis
187 //
188
189
190 if (!evtloop.PreProcess())
191 return;
192
193 int count = 1;
194
195 TH2F * shape = new TH2F("shape","shape",300,0,30,1000,-20,250);
196 TH2F * shape_corr = new TH2F("shape_corr","shape_corr",300,0,30,1000,-20,250);
197 TH2F * shape_corr_all = new TH2F("shape_corr_all","shape_corr_all",300,0,30,1000,-20,250);
198 TH1F * htime = new TH1F("htime","htime",150,0,15);
199 TH1F * hsig = new TH1F("hsig","hsig",300,-50,550);
200
201 // cout << " hello " << endl;
202
203 const Float_t ped_mean = pedcam[ipix].GetPedestal();
204 const Float_t pedRMS = pedcam[ipix].GetPedestalRms();
205
206 while (tlist2.Process()){ //loop over the events
207
208 count++;
209
210 if (count%1000==0) cout << "Event #" << count << endl;
211
212 MRawEvtPixelIter pixel(&evtdata);
213
214
215 double value = 0;
216 pedcam.GetPixelContent(value, ipix, geomcam, 0);
217
218
219 double time = 0;
220 timecam.GetPixelContent(time, ipix, geomcam, 0); // was 0)
221
222 double sig = 0;
223 sigcam.GetPixelContent(sig, ipix, geomcam, 0);
224
225 Byte_t sat = (sigcam)[ipix].GetNumHiGainSaturated();
226
227
228 htime->Fill(time);
229 hsig->Fill(sig);
230
231
232 pixel.Jump(ipix);
233
234 const Byte_t *higains = pixel.GetHiGainSamples();
235 const Byte_t *logains = pixel.GetLoGainSamples();
236 const Int_t nh = evtdata.GetNumHiGainSamples();
237 const Int_t nl = evtdata.GetNumLoGainSamples();
238
239 Bool_t ABFlag = pixel.HasABFlag();
240
241 const Byte_t *higains = pixel.GetHiGainSamples();
242 const Byte_t *logains = pixel.GetLoGainSamples();
243
244 const Float_t ABoffs = pedcam[ipix].GetPedestalABoffset();
245
246
247 Float_t PedMean[2];
248 PedMean[0] = ped_mean + ABoffs;
249 PedMean[1] = ped_mean - ABoffs;
250
251
252
253
254 for (int sample=0; sample<nh; sample++){
255 shape->Fill(sample+0.5,higains[sample]-PedMean[(sample+ABFlag)&0x1]);
256 shape->Fill(sample+15.5,logains[sample]-PedMean[(sample+ABFlag)&0x1]);
257 shape_corr->Fill(sample+0.5+5.-time,higains[sample]-PedMean[(sample+ABFlag)&0x1]);
258 shape_corr->Fill(sample+15.5+5.-time,logains[sample]-PedMean[(sample+nh+ABFlag)&0x1]);
259 if (sat==0){
260 shape_corr_all->Fill(sample+0.5+5.-time,(higains[sample]-PedMean[(sample+ABFlag)&0x1])/sig*250);
261 shape_corr_all->Fill(sample+15.5+5.-time,(logains[sample]-PedMean[(sample+nh+ABFlag)&0x1])/sig*250);
262 }
263 }
264
265
266 } // end loop over the entries
267
268 tlist2.PrintStatistics();
269
270
271 TCanvas * c_shape = new TCanvas("c_shape","c_shape",600,400);
272 c_shape->SetFillColor(0);
273 c_shape->SetBorderMode(0);
274 c_shape->SetGridx();
275 c_shape->SetGridy();
276 shape->SetStats(0);
277 TString title2 = "Raw Calibration Signal Shape, Pixel ";
278 title2+=ipix;
279 shape->SetTitle(title2);
280 shape->SetXTitle("FADC sample no.");
281 shape->SetYTitle("signal [FADC counts]");
282 shape->SetMarkerStyle(20);
283 shape->SetMarkerColor(4);
284 shape->SetMarkerSize(0.3);
285 shape->Draw();
286
287
288 TCanvas * c_shape_corr = new TCanvas("c_shape_corr","c_shape_corr",600,400);
289 c_shape_corr->SetFillColor(0);
290 c_shape_corr->SetBorderMode(0);
291 c_shape_corr->SetGridx();
292 c_shape_corr->SetGridy();
293 shape_corr->SetStats(0);
294 TString title3 = "Calibration Signal Shape, Arrival Time Corrected, Pixel ";
295 title3+=ipix;
296 shape_corr->SetTitle(title3);
297 shape_corr->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
298 shape_corr->SetYTitle("signal [FADC counts]");
299 shape_corr->SetMarkerStyle(20);
300 shape_corr->SetMarkerColor(4);
301 shape_corr->SetMarkerSize(0.3);
302 shape_corr->Draw();
303
304
305 TCanvas * c_shape_corr_all = new TCanvas("c_shape_corr_all","c_shape_corr_all",600,400);
306 c_shape_corr_all->SetFillColor(0);
307 c_shape_corr_all->SetBorderMode(0);
308 c_shape_corr_all->SetGridx();
309 c_shape_corr_all->SetGridy();
310 shape_corr_all->SetStats(0);
311 TString title3 = "Calibration Signal Shape, Arrival Time + Amplitude Corrected, Pixel ";
312 title3+=ipix;
313 shape_corr_all->SetTitle(title3);
314 shape_corr_all->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
315 shape_corr_all->SetYTitle("signal [a.u.]");
316 shape_corr_all->SetMarkerStyle(20);
317 shape_corr_all->SetMarkerColor(4);
318 shape_corr_all->SetMarkerSize(0.3);
319 shape_corr_all->Draw();
320
321
322
323
324 TCanvas * c_sig = new TCanvas("c_sig","c_sig",600,400);
325 c_sig->SetFillColor(0);
326 c_sig->SetBorderMode(0);
327 c_sig->SetGridx();
328 c_sig->SetGridy();
329
330 hsig->SetStats(0);
331 TString title4 = "Calibration Extracted Charge, MExtractFixedWindowPeakSearch, Pixel ";
332 title4+=ipix;
333 hsig->SetTitle(title4);
334 hsig->SetXTitle("reco charge [FADC counts]");
335 hsig->SetYTitle("events / 2 FADC counts");
336 hsig->SetLineColor(4);
337 hsig->SetLineWidth(2);
338 hsig->Draw();
339
340 TCanvas * c_time = new TCanvas("c_time","c_time",600,400);
341 c_time->SetFillColor(0);
342 c_time->SetBorderMode(0);
343 c_time->SetGridx();
344 c_time->SetGridy();
345
346 htime->SetStats(0);
347 TString title5 = "Calibration Arrival Time, MExtractTimeHighestIntegral, Pixel ";
348 title5+=ipix;
349 htime->SetTitle(title5);
350 htime->SetXTitle("arrival time [T_{ADC}]");
351 htime->SetYTitle("events / 0.1 T_{ADC}");
352 htime->SetLineColor(4);
353 htime->SetLineWidth(2);
354 htime->Draw();
355
356
357 TCanvas * c7 = new TCanvas("c7","c7",600,400);
358 c7->cd();
359 c7->SetGridx();
360 c7->SetGridy();
361 c7->SetFillColor(0);
362 // c7->SetBordermode(0);
363
364 TProfile * shape_corr_all_pfx = (TProfile*)shape_corr_all->ProfileX();
365
366 shape_corr_all_pfx->SetStats(0);
367 TString title6 = "Average Calibration Signal Shape, Amplitude + Arrival Time Corrected, Pixel";
368 title6+=ipix;
369 shape_corr_all_pfx->SetTitle(title6);
370 shape_corr_all_pfx->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
371 shape_corr_all_pfx->SetYTitle("signal [a.u.]");
372 shape_corr_all_pfx->SetMarkerStyle(20);
373 shape_corr_all_pfx->SetMarkerColor(4);
374 shape_corr_all_pfx->SetMarkerSize(0.5);
375 shape_corr_all_pfx->Draw();
376}
377
378
379
380
Note: See TracBrowser for help on using the repository browser.