source: trunk/MagicSoft/Mars/macros/pedestalstudies.C@ 3651

Last change on this file since 3651 was 3316, checked in by gaug, 22 years ago
*** empty log message ***
File size: 11.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
25//const TString pedfile = "/remote/home/pc2/operator/Crab20040214/20040215_16743_P_CrabOn_E.root";
26const TString pedfile = "../20040215_16770_P_OffCrab4_E.root";
27
28//const TString pedfile = "/mnt/users/mdoro/Mars/Data/20040201_14418_P_OffMrk421-1_E.root";
29//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_02_10/20040210_14607_P_CrabNebula_E.root";
30//const TString pedfile = "/mnt/Data/rootdata/CrabNebula/2004_01_26/20040125_10412_P_Crab-On_E.root";
31//const TString pedfile = "/mnt/Data/rootdata/Miscellaneous/2003_12_19/20031218_03522_P_Park_E.root";
32
33void pedestalstudies(TString pedname=pedfile)
34{
35
36 gStyle->SetOptStat(1111);
37 gStyle->SetOptFit();
38
39 MStatusDisplay *display = new MStatusDisplay;
40 display->SetUpdateTime(500);
41 display->Resize(850,700);
42
43 //
44 // Create a empty Parameter List and an empty Task List
45 // The tasklist is identified in the eventloop by its name
46 //
47 MParList plist;
48 MTaskList tlist;
49 plist.AddToList(&tlist);
50
51 //
52 // Now setup the tasks and tasklist for the pedestals:
53 // ---------------------------------------------------
54 //
55
56 MReadMarsFile read("Events", pedname);
57 read.DisableAutoScheme();
58
59 MGeomApply geomapl;
60 MExtractSignal sigcalc;
61
62 //
63 // Set the extraction range higher:
64 //
65 //sigcalc.SetRange(1,14,1,14);
66
67 MPedCalcPedRun pedcalc;
68
69 //
70 // Additionally to calculating the pedestals,
71 // you can fill histograms and look at them
72 //
73 MFillH fill("MHPedestalCam", "MExtractedSignalCam");
74
75 tlist.AddToList(&read);
76 tlist.AddToList(&geomapl);
77 tlist.AddToList(&sigcalc);
78 tlist.AddToList(&pedcalc);
79 tlist.AddToList(&fill);
80
81 MGeomCamMagic geomcam;
82 MPedestalCam pedcam;
83 MHPedestalCam hpedcam;
84 plist.AddToList(&geomcam);
85 plist.AddToList(&pedcam);
86 plist.AddToList(&hpedcam);
87
88 //
89 // Create and setup the eventloop
90 //
91 MEvtLoop evtloop;
92
93 evtloop.SetParList(&plist);
94 evtloop.SetDisplay(display);
95
96 //
97 // Execute first analysis
98 //
99 if (!evtloop.Eventloop())
100 return;
101
102 tlist.PrintStatistics();
103
104 //
105 // Look at one specific pixel, after all the histogram manipulations:
106 //
107// hpedcam[9].DrawClone("fourierevents");
108
109
110 MHCamera dispped0 (geomcam, "Ped;Pedestal", "Mean per Slice");
111 MHCamera dispped1 (geomcam, "Ped;PedestalErr", "Mean Error per Slice");
112 MHCamera dispped2 (geomcam, "Ped;PedestalRms", "RMS per Slice");
113 MHCamera dispped3 (geomcam, "Ped;PedestalRmsErr", "RMS Error per Slice");
114
115 MHCamera dispped4 (geomcam, "Ped;Mean", "Fitted Mean per Slice");
116 MHCamera dispped5 (geomcam, "Ped;MeanErr", "Fitted Error of Mean per Slice");
117 MHCamera dispped6 (geomcam, "Ped;Sigma", "Fitted Sigma per Slice");
118 MHCamera dispped7 (geomcam, "Ped;SigmaErr", "Fitted Error of Sigma per Slice");
119 MHCamera dispped8 (geomcam, "Ped;Prob", "Probability of Fit");
120 MHCamera dispped9 (geomcam, "Ped;DeltaPedestalMean", "Rel. Diff. Mean per Slice (Calc.-Fitte)");
121 MHCamera dispped10 (geomcam, "Ped;DeltaPedestalMeanError", "Rel. Diff. Mean Error per Slice (Calc.-Fitted)");
122 MHCamera dispped11 (geomcam, "Ped;DeltaRmsSigma", "Rel. Diff. RMS per Slice (Calc.-Fitted)");
123 MHCamera dispped12 (geomcam, "Ped;DeltaRmsSigmaError", "Rel. Diff. RMS Error per Slice (Calc.-Fitted)");
124 MHCamera dispped13 (geomcam, "Ped;FitOK", "Gaus Fit not OK");
125 MHCamera dispped14 (geomcam, "Ped;FourierOK", "Fourier Analysis not OK");
126
127 dispped0.SetCamContent( pedcam, 0);
128 dispped0.SetCamError( pedcam, 1);
129 dispped1.SetCamContent( pedcam, 1);
130 dispped2.SetCamContent( pedcam, 2);
131 dispped2.SetCamError( pedcam, 3);
132 dispped3.SetCamContent( pedcam, 3);
133
134 dispped4.SetCamContent( hpedcam, 0);
135 dispped4.SetCamError( hpedcam, 1);
136 dispped5.SetCamContent( hpedcam, 1);
137 dispped6.SetCamContent( hpedcam, 2);
138 dispped6.SetCamError( hpedcam, 3);
139 dispped7.SetCamContent( hpedcam, 3);
140 dispped8.SetCamContent( hpedcam, 4);
141 dispped9.SetCamContent( hpedcam, 5);
142 dispped9.SetCamError( hpedcam, 6);
143 dispped10.SetCamContent(hpedcam, 7);
144 dispped11.SetCamContent(hpedcam, 8);
145 dispped11.SetCamError( hpedcam, 9);
146 dispped12.SetCamContent(hpedcam, 10);
147 dispped13.SetCamContent(hpedcam, 11);
148 dispped14.SetCamContent(hpedcam, 12);
149
150 dispped0.SetYTitle("Calc. Pedestal per slice [FADC counts]");
151 dispped1.SetYTitle("Calc. Pedestal Error per slice [FADC counts]");
152 dispped2.SetYTitle("Calc. Pedestal RMS per slice [FADC counts]");
153 dispped3.SetYTitle("Calc. Pedestal RMS Error per slice [FADC counts]");
154 dispped4.SetYTitle("Fitted Mean per slice [FADC counts]");
155 dispped5.SetYTitle("Error of Fitted Mean per slice [FADC counts]");
156 dispped6.SetYTitle("Fitted Sigma per slice [FADC counts]");
157 dispped7.SetYTitle("Error of Fitted Sigma per slice [FADC counts]");
158 dispped8.SetYTitle("Fit Probability [1]");
159 dispped9.SetYTitle("Rel. Diff. Pedestal Calc.-Fitted per slice [1]");
160 dispped10.SetYTitle("Rel. Diff. Pedestal Error Calc.-Fitted per slice [1]");
161 dispped11.SetYTitle("Rel. Diff. Pedestal RMS Calc.-Fitted per slice [1]");
162 dispped12.SetYTitle("Rel. Diff. Pedestal RMS Error Calc.-Fitted per slice [1]");
163 dispped13.SetYTitle("[1]");
164 dispped14.SetYTitle("[1]");
165
166 // Histogram values
167 TCanvas &b1 = display->AddTab("Ped.Calc.");
168 b1.Divide(4,3);
169
170 CamDraw(b1,dispped0,pedcam,1,4,1);
171 CamDraw(b1,dispped1,pedcam,2,4,2);
172 CamDraw(b1,dispped2,pedcam,3,4,2);
173 CamDraw(b1,dispped3,pedcam,4,4,2);
174
175 // Fitted values
176 TCanvas &b2 = display->AddTab("Ped.Fit");
177 b2.Divide(4,3);
178
179 CamDraw(b2,dispped4,hpedcam,1,4,1);
180 CamDraw(b2,dispped5,hpedcam,2,4,2);
181 CamDraw(b2,dispped6,hpedcam,3,4,2);
182 CamDraw(b2,dispped7,hpedcam,4,4,2);
183
184
185 // Fits Probability
186 TCanvas &b3 = display->AddTab("Ped.Fit Prob.");
187 b3.Divide(1,3);
188
189 CamDraw(b3,dispped8,hpedcam,1,1,3);
190
191 // Differences
192 TCanvas &c4 = display->AddTab("Rel.Diff.Calc.-Fit");
193 c4.Divide(4,3);
194
195 CamDraw(c4,dispped9,hpedcam,1,4,1);
196 CamDraw(c4,dispped10,hpedcam,2,4,1);
197 CamDraw(c4,dispped11,hpedcam,3,4,1);
198 CamDraw(c4,dispped12,hpedcam,4,4,1);
199
200 // Defects
201 TCanvas &c5 = display->AddTab("Defects");
202 c5.Divide(2,2);
203
204 CamDraw(c5,dispped13,hpedcam,1,2,0);
205 CamDraw(c5,dispped14,hpedcam,2,2,0);
206
207}
208
209void CamDraw(TCanvas &c, MHCamera &cam, MCamEvent &evt, Int_t i, Int_t j, Int_t fit)
210{
211
212 c.cd(i);
213 gPad->SetBorderMode(0);
214 MHCamera *obj1=(MHCamera*)cam.DrawCopy("hist");
215 // obj1->AddNotify(evt);
216
217 c.cd(i+j);
218 gPad->SetBorderMode(0);
219 obj1->Draw();
220 ((MHCamera*)obj1)->SetPrettyPalette();
221
222 if (fit != 0)
223 {
224 c.cd(i+2*j);
225 gPad->SetBorderMode(0);
226 TH1D *obj2 = (TH1D*)obj1->Projection();
227
228// obj2->Sumw2();
229 obj2->Draw();
230 obj2->SetBit(kCanDelete);
231
232 const Double_t min = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
233 const Double_t max = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
234 const Double_t integ = obj2->Integral("width")/2.5066283;
235 const Double_t mean = obj2->GetMean();
236 const Double_t rms = obj2->GetRMS();
237 const Double_t width = max-min;
238
239 if (rms == 0. || width == 0. )
240 return;
241
242 switch (fit)
243 {
244 case 1:
245 TF1 *sgaus = new TF1("sgaus","gaus(0)",min,max);
246 sgaus->SetBit(kCanDelete);
247 sgaus->SetParNames("Area","#mu","#sigma");
248 sgaus->SetParameters(integ/rms,mean,rms);
249 sgaus->SetParLimits(0,0.,integ);
250 sgaus->SetParLimits(1,min,max);
251 sgaus->SetParLimits(2,0,width/1.5);
252 obj2->Fit("sgaus","QLR");
253 obj2->GetFunction("sgaus")->SetLineColor(kYellow);
254 break;
255
256 case 2:
257 TString dgausform = "([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
258 dgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
259 TF1 *dgaus = new TF1("dgaus",dgausform.Data(),min,max);
260 dgaus->SetBit(kCanDelete);
261 dgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}","A_{2}","#mu_{2}","#sigma_{2}");
262 dgaus->SetParameters(integ,(min+mean)/2.,width/4.,
263 integ/width/2.,(max+mean)/2.,width/4.);
264 // The left-sided Gauss
265 dgaus->SetParLimits(0,integ-1.5,integ+1.5);
266 dgaus->SetParLimits(1,min+(width/10.),mean);
267 dgaus->SetParLimits(2,0,width/2.);
268 // The right-sided Gauss
269 dgaus->SetParLimits(3,0,integ);
270 dgaus->SetParLimits(4,mean,max-(width/10.));
271 dgaus->SetParLimits(5,0,width/2.);
272 obj2->Fit("dgaus","QLRM");
273 obj2->GetFunction("dgaus")->SetLineColor(kYellow);
274 break;
275
276 case 3:
277 TString tgausform = "([0]-[3]-[6])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])";
278 tgausform += "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])";
279 tgausform += "+[6]/[8]*exp(-0.5*(x-[7])*(x-[7])/[8]/[8])";
280 TF1 *tgaus = new TF1("tgaus",tgausform.Data(),min,max);
281 tgaus->SetBit(kCanDelete);
282 tgaus->SetParNames("A_{tot}","#mu_{1}","#sigma_{1}",
283 "A_{2}","#mu_{2}","#sigma_{2}",
284 "A_{3}","#mu_{3}","#sigma_{3}");
285 tgaus->SetParameters(integ,(min+mean)/2,width/4.,
286 integ/width/3.,(max+mean)/2.,width/4.,
287 integ/width/3.,mean,width/2.);
288 // The left-sided Gauss
289 tgaus->SetParLimits(0,integ-1.5,integ+1.5);
290 tgaus->SetParLimits(1,min+(width/10.),mean);
291 tgaus->SetParLimits(2,width/15.,width/2.);
292 // The right-sided Gauss
293 tgaus->SetParLimits(3,0.,integ);
294 tgaus->SetParLimits(4,mean,max-(width/10.));
295 tgaus->SetParLimits(5,width/15.,width/2.);
296 // The Gauss describing the outliers
297 tgaus->SetParLimits(6,0.,integ);
298 tgaus->SetParLimits(7,min,max);
299 tgaus->SetParLimits(8,width/4.,width/1.5);
300 obj2->Fit("tgaus","QLRM");
301 obj2->GetFunction("tgaus")->SetLineColor(kYellow);
302 break;
303 case 4:
304 obj2->Fit("pol0","Q");
305 obj2->GetFunction("pol0")->SetLineColor(kYellow);
306 break;
307 case 9:
308 break;
309 default:
310 obj2->Fit("gaus","Q");
311 obj2->GetFunction("gaus")->SetLineColor(kYellow);
312 break;
313 }
314
315 gPad->Modified();
316 gPad->Update();
317
318 }
319}
Note: See TracBrowser for help on using the repository browser.