source: releases/Mars.2014.05.26/mtemp/mmpi/macros/data_shape.C@ 18027

Last change on this file since 18027 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 = 36058;
32Int_t calr = 36058;
33//Int_t datar = 12517;
34
35
36
37void data_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(WindowSize,WindowSize);
149
150
151 // timecalc->SetRange(8,14,8,14,6);
152
153 // filter against cosmics
154 MFCosmics cosmics;
155 MContinue cont(&cosmics);
156
157 tlist2.AddToList(&read2);
158 tlist2.AddToList(&geomapl);
159 tlist2.AddToList(&extractor);
160 tlist2.AddToList(&timeext);
161
162 //
163 // In case, you want to skip the cosmics rejection,
164 // uncomment the next line
165 //
166 tlist2.AddToList(&cont);
167 //
168 // In case, you want to skip the somewhat lengthy calculation
169 // of the arrival times using a spline, uncomment the next two lines
170 //
171
172
173 //
174 // Create and setup the eventloop
175 //
176 MEvtLoop evtloop;
177 evtloop.SetParList(&plist2);
178 // evtloop.SetDisplay(display);
179
180 cout << "***************************" << endl;
181 cout << "** COMPUTING CALIBRATION **" << endl;
182 cout << "***************************" << endl;
183
184 //
185 // Execute second analysis
186 //
187
188
189 if (!evtloop.PreProcess())
190 return;
191
192 int count = 1;
193
194 TH2F * shape = new TH2F("shape","shape",30,0,30,1000,-50,250);
195 TH2F * shape_corr = new TH2F("shape_corr","shape_corr",300,0,30,1000,-50,250);
196 TH2F * shape_corr_all = new TH2F("shape_corr_all","shape_corr_all",300,0,30,1000,-50,250);
197 TH1F * htime = new TH1F("htime","htime",150,0,15);
198 TH1F * hsig = new TH1F("hsig","hsig",300,-50,550);
199
200 // cout << " hello " << endl;
201
202 const Float_t ped_mean = pedcam[ipix].GetPedestal();
203 const Float_t pedRMS = pedcam[ipix].GetPedestalRms();
204
205 while (tlist2.Process()){ //loop over the events
206
207 count++;
208
209 if (count%1000==0) cout << "Event #" << count << endl;
210
211 MRawEvtPixelIter pixel(&evtdata);
212
213
214 double value = 0;
215 pedcam.GetPixelContent(value, ipix, geomcam, 0);
216
217
218 double time = 0;
219 timecam.GetPixelContent(time, ipix, geomcam, 0); // was 0)
220
221 double sig = 0;
222 sigcam.GetPixelContent(sig, ipix, geomcam, 0);
223
224 Byte_t sat = (sigcam)[ipix].GetNumHiGainSaturated();
225
226 hsig->Fill(sig);
227
228 if (sig>60 && sat==0){
229
230 htime->Fill(time);
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 Float_t PedMean[2];
247 PedMean[0] = ped_mean + ABoffs;
248 PedMean[1] = ped_mean - ABoffs;
249
250
251
252
253 for (int sample=0; sample<nh; sample++){
254 shape->Fill(sample+0.5,higains[sample]-PedMean[(sample+ABFlag)&0x1]);
255 shape->Fill(sample+15.5,logains[sample]-PedMean[(sample+ABFlag)&0x1]);
256 shape_corr->Fill(sample+0.5+6.-time,higains[sample]-PedMean[(sample+ABFlag)&0x1]);
257 shape_corr->Fill(sample+15.5+6.-time,logains[sample]-PedMean[(sample+nh+ABFlag)&0x1]);
258 shape_corr_all->Fill(sample+0.5+6.-time,(higains[sample]-PedMean[(sample+ABFlag)&0x1])/sig*250);
259 shape_corr_all->Fill(sample+15.5+6.-time,(logains[sample]-PedMean[(sample+nh+ABFlag)&0x1])/sig*250);
260 }
261 }// end if (sig > 50)
262
263 } // end loop over the entries
264
265 tlist2.PrintStatistics();
266
267
268 TCanvas * c_shape = new TCanvas("c_shape","c_shape",600,400);
269 c_shape->SetFillColor(0);
270 c_shape->SetBorderMode(0);
271 c_shape->SetGridx();
272 c_shape->SetGridy();
273 shape->SetStats(0);
274 TString title2 = "Raw Data Signal Shape, Pixel ";
275 title2+=ipix;
276 shape->SetTitle(title2);
277 shape->SetXTitle("FADC sample no.");
278 shape->SetYTitle("signal [FADC counts]");
279 shape->SetMarkerStyle(20);
280 shape->SetMarkerColor(4);
281 shape->SetMarkerSize(0.3);
282 shape->Draw();
283
284
285 TCanvas * c_shape_corr = new TCanvas("c_shape_corr","c_shape_corr",600,400);
286 c_shape_corr->SetFillColor(0);
287 c_shape_corr->SetBorderMode(0);
288 c_shape_corr->SetGridx();
289 c_shape_corr->SetGridy();
290 shape_corr->SetStats(0);
291 TString title3 = "Data Signal Shape, Arrival Time Corrected, Pixel ";
292 title3+=ipix;
293 shape_corr->SetTitle(title3);
294 shape_corr->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
295 shape_corr->SetYTitle("signal [FADC counts]");
296 shape_corr->SetMarkerStyle(20);
297 shape_corr->SetMarkerColor(4);
298 shape_corr->SetMarkerSize(0.3);
299 shape_corr->Draw();
300
301
302 TCanvas * c_shape_corr_all = new TCanvas("c_shape_corr_all","c_shape_corr_all",600,400);
303 c_shape_corr_all->SetFillColor(0);
304 c_shape_corr_all->SetBorderMode(0);
305 c_shape_corr_all->SetGridx();
306 c_shape_corr_all->SetGridy();
307 shape_corr_all->SetStats(0);
308 TString title3 = "Data Signal Shape, Arrival Time + Amplitude Corrected, Pixel ";
309 title3+=ipix;
310 shape_corr_all->SetTitle(title3);
311 shape_corr_all->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
312 shape_corr_all->SetYTitle("signal [a.u.]");
313 shape_corr_all->SetMarkerStyle(20);
314 shape_corr_all->SetMarkerColor(4);
315 shape_corr_all->SetMarkerSize(0.3);
316 shape_corr_all->Draw();
317
318
319
320
321 TCanvas * c_sig = new TCanvas("c_sig","c_sig",600,400);
322 c_sig->SetFillColor(0);
323 c_sig->SetBorderMode(0);
324 c_sig->SetGridx();
325 c_sig->SetGridy();
326
327 hsig->SetStats(0);
328 TString title4 = "Data Extracted Charge, MExtractFixedWindowPeakSearch, Pixel ";
329 title4+=ipix;
330 hsig->SetTitle(title4);
331 hsig->SetXTitle("reco charge [FADC counts]");
332 hsig->SetYTitle("events / 2 FADC counts");
333 hsig->SetLineColor(4);
334 hsig->SetLineWidth(2);
335 hsig->Draw();
336
337 TCanvas * c_time = new TCanvas("c_time","c_time",600,400);
338 c_time->SetFillColor(0);
339 c_time->SetBorderMode(0);
340 c_time->SetGridx();
341 c_time->SetGridy();
342
343 htime->SetStats(0);
344 TString title5 = "Data Arrival Time, MExtractTimeHighestIntegral, Pixel ";
345 title5+=ipix;
346 htime->SetTitle(title5);
347 htime->SetXTitle("arrival time [T_{ADC}]");
348 htime->SetYTitle("events / 0.1 T_{ADC}");
349 htime->SetLineColor(4);
350 htime->SetLineWidth(2);
351 htime->Draw();
352
353 TCanvas * c7 = new TCanvas("c7","c7",600,400);
354 c7->cd();
355 c7->SetGridx();
356 c7->SetGridy();
357 c7->SetFillColor(0);
358 // c7->SetBordermode(0);
359
360 TProfile * shape_corr_all_pfx = (TProfile*)shape_corr_all->ProfileX();
361
362 shape_corr_all_pfx->SetStats(0);
363 TString title6 = "Average Data Signal Shape, Amplitude + Arrival Time Corrected, Pixel";
364 title6+=ipix;
365 shape_corr_all_pfx->SetTitle(title6);
366 shape_corr_all_pfx->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
367 shape_corr_all_pfx->SetYTitle("signal [a.u.]");
368 shape_corr_all_pfx->SetMarkerStyle(20);
369 shape_corr_all_pfx->SetMarkerColor(4);
370 shape_corr_all_pfx->SetMarkerSize(0.5);
371 shape_corr_all_pfx->Draw();
372
373}
374
375
376
377
Note: See TracBrowser for help on using the repository browser.