source: trunk/Mars/mtemp/mmpi/macros/noise_autocorrelation_AB.C

Last change on this file was 5724, checked in by hbartko, 20 years ago
*** empty log message ***
File size: 9.1 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 09/2004 <mailto:hbartko@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2001
21!
22!
23\* ======================================================================== */
24
25
26const TString defname = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_03/20040903_35698_P_multiplexFADC_E.root";
27const TString defpedname = "/.magic/magicserv01/scratch/hbartko/rootdata/2004_09_03/20040903_35698_P_multiplexFADC_E.root";
28
29// works with one pedestal file to calculate the noise autocorrelation
30
31void noise_autocorrelation_AB(const TString fname = defpedname, const TString pedname = defpedname, Int_t ipix = 67) {
32
33 MParList plist_ped;
34
35 MTaskList tlist_ped;
36 plist_ped.AddToList(&tlist_ped);
37
38 MPedestalCam pedcam;
39 plist_ped.AddToList(&pedcam);
40
41 MReadMarsFile read("Events", pedname);
42 read.DisableAutoScheme();
43 tlist_ped.AddToList(&read);
44
45 MGeomApply geomapl_ped;
46 MGeomCamMagic geomcam;
47 tlist_ped.AddToList(&geomapl_ped);
48
49
50 MPedCalcFromLoGain pedcalc_ped;
51 //pedcalc_ped.SetMaxHiGainVar(20);
52 //pedcalc_ped.SetRange(0, 11, 1, 14);
53 //pedcalc_ped.SetWindowSize(12,14);
54 pedcalc_ped.SetPedestalUpdate(kFALSE);
55 tlist_ped.AddToList(&pedcalc_ped);
56
57 MEvtLoop evtloop_ped;
58 evtloop_ped.SetParList(&plist_ped);
59
60 if (!evtloop_ped.Eventloop())
61 return;
62
63 tlist_ped.PrintStatistics();
64
65
66 // now the event loop for the signal reconstruction with pedestals subtracted
67
68
69 MParList plist;
70 MTaskList tlist;
71 plist.AddToList(&tlist);
72
73 plist.AddToList(&pedcam);
74
75 MRawRunHeader runheader;
76 plist.AddToList(&runheader);
77
78 MRawEvtData evtdata;
79 plist.AddToList(&evtdata);
80
81 MExtractedSignalCam sigcam;
82 plist.AddToList(&sigcam);
83
84 MArrivalTimeCam timecam;
85 plist.AddToList(&timecam);
86
87
88 MReadMarsFile read("Events", fname);
89 read.DisableAutoScheme();
90
91 MGeomApply geomapl;
92
93 Int_t WindowSize = 6; // size of the chosen integration window for signal and time extraction
94
95 MExtractFixedWindowPeakSearch extractor;
96 extractor.SetRange(0, 14, 0, 14);
97 extractor.SetWindows( WindowSize, WindowSize, 4);
98
99 MExtractTimeHighestIntegral timeext;
100 timeext.SetRange(0, 14, 0, 14);
101 timeext.SetWindowSize(WindowSize,WindowSize);
102
103
104 tlist.AddToList(&read);
105 tlist.AddToList(&geomapl);
106 tlist.AddToList(&extractor);
107 tlist.AddToList(&timeext);
108
109
110 MEvtLoop evtloop;
111 evtloop.SetParList(&plist);
112
113 if (!evtloop.PreProcess())
114 return;
115
116 TString title = "Noise Autocorrelation B_{ij} = <b_{i}*b_{j}> - <b_{i}>*<b_{j}>, Pixel ";
117 title += ipix;
118
119
120 Int_t First = 1;
121
122 TH2F * hcorr = new TH2F("hcorr","hcorr",15,0,15,15,0,15);
123
124 TH2F * shape = new TH2F("shape","shape",311,-1.05,30.05,1000,-20,250);
125 TH2F * shape_corr = new TH2F("shape_corr","shape_corr",311,-1.05,30.05,1000,-20,250);
126 TH1F * htime = new TH1F("htime","htime",150,0,15);
127 TH1F * hsig = new TH1F("hsig","hsig",300,-50,550);
128
129
130 Float_t B[15][15];
131 Float_t bi[15];
132 Float_t bj[15];
133
134 for (int i=0; i<15; i++){
135 bi[i]=0;
136 bj[i]=0;
137 for (int j=0; j<15; j++){
138 B[i][j]=0;
139 }
140 }
141
142 const Float_t ped_mean = pedcam[ipix].GetPedestal();
143 const Float_t pedRMS = pedcam[ipix].GetPedestalRms();
144
145 cout << " ped_mean " << ped_mean << " pedRMS " << pedRMS << endl;
146 int count = 0;
147
148 while (tlist.Process()) {
149
150 if (First) {
151 First = 0;
152
153 const Int_t nh = runheader.GetNumSamplesHiGain();
154 const Int_t nl = runheader.GetNumSamplesLoGain();
155 const Int_t nt = nh+nl;
156
157 }
158
159
160 htime->Fill(timecam[ipix].GetArrivalTimeHiGain());
161 hsig->Fill(sigcam[ipix].GetExtractedSignalHiGain());
162
163
164 MRawEvtPixelIter pixel(&evtdata);
165 pixel.Jump(ipix);
166
167 Bool_t ABFlag = pixel.HasABFlag();
168
169
170 count++; // couter for the number of events
171
172 const Byte_t *higains = pixel.GetHiGainSamples();
173 const Byte_t *logains = pixel.GetLoGainSamples();
174
175 const Float_t ABoffs = pedcam[ipix].GetPedestalABoffset();
176
177 Float_t PedMean[2];
178 PedMean[0] = ped_mean + ABoffs;
179 PedMean[1] = ped_mean - ABoffs;
180
181
182 // calculate the noise autocorrelation matrix
183 for (int slice=0; slice<nh; slice++) {
184 bi[slice]+= higains[slice]-PedMean[(slice+ABFlag)&0x1];
185 bj[slice]+= higains[slice]-PedMean[(slice+ABFlag)&0x1];
186 for (int ii = 0; ii< nh; ii++){
187 B[slice][ii]+=(higains[slice]-PedMean[(slice+ABFlag)&0x1])*(higains[ii]-PedMean[(ii+ABFlag)&0x1]);
188 }
189 }
190
191
192 // compute the noise signal shape
193
194 Int_t trigger = 0;
195 Double_t charge = 0;
196 Double_t charge_time = 0;
197
198
199 for (int slice=0; slice<nh; slice++) {
200 shape->Fill(slice,higains[slice]-PedMean[(slice+ABFlag)&0x1]);
201 shape->Fill(slice+nh,logains[slice]-PedMean[(slice+ABFlag)&0x1]);
202 if (trigger==0 && (higains[slice]-PedMean[(slice+ABFlag)&0x1]) >= 4.0)
203 trigger = slice;
204 }
205
206 for (int slice=0; slice<nh; slice++)
207 if (trigger>=1 && slice<=trigger+2){
208 charge += higains[slice]-PedMean[(slice+ABFlag)&0x1];
209 charge_time += slice*(higains[slice]-PedMean[(slice+ABFlag)&0x1]);
210 }
211
212 if (charge>0) charge_time /= charge;
213
214 for (int slice=0; slice<nh; slice++)
215 if (trigger>=1)
216 shape_corr->Fill(slice+3-charge_time,(higains[slice]-PedMean[(slice+ABFlag)&0x1]));
217
218
219
220
221 } // end eventloop.Process()
222
223 evtloop.PostProcess();
224
225
226 for (int i=0; i<15; i++){
227 for (int j=0; j<15; j++){
228 hcorr->Fill(i,j,B[i][j]/count-bi[i]/count*bj[j]/count);
229 // cout << "i " << i << " j " << j << " B[i][j] " << B[i][j] << " bi[i] " << bi[i] << " bj[j] " << bj[j] << endl;
230 }
231 }
232
233 cout << " pedmean " << ped_mean << " pedRMS " << pedRMS << endl;
234
235
236 TCanvas * c_corr = new TCanvas("c_corr","c_corr",600,400);
237 c_corr->SetFillColor(0);
238 c_corr->SetBorderMode(0);
239 c_corr->SetGridx();
240 c_corr->SetGridy();
241 gStyle->SetPalette(1);
242
243 hcorr->SetStats(0);
244 hcorr->SetTitle(title);
245 hcorr->SetXTitle("slice i");
246 hcorr->SetYTitle("slice j");
247 hcorr->Draw("colz");
248
249
250 TCanvas * c_corr_px = new TCanvas("c_corr_px","c_corr_px",600,400);
251 c_corr_px->SetFillColor(0);
252 c_corr_px->SetBorderMode(0);
253 c_corr_px->SetGridx();
254 c_corr_px->SetGridy();
255
256 TH1F * hcorr_px = (TH1F*)hcorr->ProjectionX("hcorr_px",6,6);
257 hcorr_px->SetLineColor(4);
258 hcorr_px->SetLineWidth(2);
259 hcorr_px->SetStats(0);
260 TString title1 = "Noise Autocorrelation, FADC sample 6, Pixel";
261 title1+=ipix;
262 hcorr_px->SetTitle(title);
263 hcorr_px->SetXTitle("slice i");
264 hcorr_px->SetYTitle("B_{ij} = <b_{i}*b_{j}> - <b_{i}>*<b_{j}> [FADC counts]");
265
266 hcorr_px->Draw();
267
268
269
270 TCanvas * c_shape = new TCanvas("c_shape","c_shape",600,400);
271 c_shape->SetFillColor(0);
272 c_shape->SetBorderMode(0);
273 c_shape->SetGridx();
274 c_shape->SetGridy();
275 shape->SetStats(0);
276 TString title2 = "Noise Signal Shape, Signal #geq 3 FADC counts, Pixel ";
277 title2+=ipix;
278 shape->SetTitle(title2);
279 shape->SetXTitle("FADC sample no.");
280 shape->SetYTitle("signal [FADC counts]");
281 shape->SetMarkerStyle(20);
282 shape->SetMarkerColor(4);
283 shape->SetMarkerSize(0.3);
284 shape->Draw();
285
286
287 TCanvas * c_shape_corr = new TCanvas("c_shape_corr","c_shape_corr",600,400);
288 c_shape_corr->SetFillColor(0);
289 c_shape_corr->SetBorderMode(0);
290 c_shape_corr->SetGridx();
291 c_shape_corr->SetGridy();
292 shape_corr->SetStats(0);
293 TString title3 = "Noise Signal Shape, Signal #geq 3 FADC counts, Pixel ";
294 title3+=ipix;
295 shape_corr->SetTitle(title3);
296 shape_corr->SetXTitle("FADC sample no. + (<t_{arrival}> - t_{arrival}) / T_{ADC}");
297 shape_corr->SetYTitle("signal [FADC counts]");
298 shape_corr->SetMarkerStyle(20);
299 shape_corr->SetMarkerColor(4);
300 shape_corr->SetMarkerSize(0.3);
301 shape_corr->Draw();
302
303
304 TCanvas * c_sig = new TCanvas("c_sig","c_sig",600,400);
305 c_sig->SetFillColor(0);
306 c_sig->SetBorderMode(0);
307 c_sig->SetGridx();
308 c_sig->SetGridy();
309
310 hsig->SetStats(0);
311 TString title4 = "Noise Extracted Charge, MExtractFixedWindowPeakSearch, Pixel ";
312 title4+=ipix;
313 hsig->SetTitle(title4);
314 hsig->SetXTitle("reco charge [FADC counts]");
315 hsig->SetYTitle("events / 2 FADC counts");
316 hsig->SetLineColor(4);
317 hsig->SetLineWidth(2);
318 hsig->Draw();
319
320 TCanvas * c_time = new TCanvas("c_time","c_time",600,400);
321 c_time->SetFillColor(0);
322 c_time->SetBorderMode(0);
323 c_time->SetGridx();
324 c_time->SetGridy();
325
326 htime->SetStats(0);
327 TString title5 = "Noise Arrival Time, MExtractTimeHighestIntegral, Pixel ";
328 title5+=ipix;
329 htime->SetTitle(title5);
330 htime->SetXTitle("arrival time [T_{ADC}]");
331 htime->SetYTitle("events / 0.1 T_{ADC}");
332 htime->SetLineColor(4);
333 htime->SetLineWidth(2);
334 htime->Draw();
335
336}
337
338
339
340
341
342
343
Note: See TracBrowser for help on using the repository browser.