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 | ! Author(s): Javier Lopez, 12/2003 <mailto:jlopez@ifae.es>
|
---|
18 | ! Author(s): Alex Armada, 1/2004 <mailto:armada@ifae.es>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2003
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 |
|
---|
26 | //------------------------------------------------------------------------- //
|
---|
27 | // //
|
---|
28 | // This macro fits the dc signal of a star using a two dimension gaussian //
|
---|
29 | // for each dc measurement. Then the values of parameters of the fit are //
|
---|
30 | // stored in histograms and shown at the end of the macro. //
|
---|
31 | // //
|
---|
32 | // USAGE: //
|
---|
33 | // It has two arguments, //
|
---|
34 | // 1- The first one is the dc file with the tracked star //
|
---|
35 | // 2- The second one is a continuos light file used to intercalibrate //
|
---|
36 | // the gain of the photomultipliers. //
|
---|
37 | // (It's possible not to use the calibration file and then the gain //
|
---|
38 | // of the pmts are supouse to be the same for all of them. //
|
---|
39 | // 3- The third argument is just the number of dc measurements you want //
|
---|
40 | // analize. If you put 0 it just stops after each fit and show you //
|
---|
41 | // results. //
|
---|
42 | // //
|
---|
43 | //--------------------------------------------------------------------------//
|
---|
44 |
|
---|
45 | const Int_t numPixels = 577;
|
---|
46 | Int_t nPixelsFitted;
|
---|
47 | Bool_t isPixelsFitted[numPixels];
|
---|
48 | Float_t z[numPixels],x[numPixels],y[numPixels],errorz[numPixels];
|
---|
49 | Float_t chisquare;
|
---|
50 |
|
---|
51 | //______________________________________________________________________________
|
---|
52 | //
|
---|
53 | // Function used by Minuit to do the fit
|
---|
54 | //
|
---|
55 | void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
|
---|
56 | {
|
---|
57 | Int_t i;
|
---|
58 |
|
---|
59 | //calculate chisquare
|
---|
60 | Double_t chisq = 0;
|
---|
61 | Double_t delta;
|
---|
62 | nPixelsFitted=0;
|
---|
63 | for (i=1;i<numPixels; i++) {
|
---|
64 | if (isPixelsFitted[i])
|
---|
65 | {
|
---|
66 | if (errorz[i] != 0.0)
|
---|
67 | {
|
---|
68 | delta = (z[i]-func(x[i],y[i],par))/errorz[i];
|
---|
69 | chisq += delta*delta;
|
---|
70 | nPixelsFitted++;
|
---|
71 | }
|
---|
72 | else
|
---|
73 | cout << "This should never happen errorz[" << i << "] " << errorz[i] << endl;
|
---|
74 | }
|
---|
75 | }
|
---|
76 | f = chisq;
|
---|
77 | chisquare = chisq;
|
---|
78 | }
|
---|
79 |
|
---|
80 | //______________________________________________________________________________
|
---|
81 | //
|
---|
82 | // The 2D gaussian fucntion used to fit the spot of the star
|
---|
83 | //
|
---|
84 | Double_t func(float x,float y,Double_t *par)
|
---|
85 | {
|
---|
86 | Double_t value=par[0]*TMath::exp(-(x-par[1])*(x-par[1])/(2*par[2]*par[2]))*TMath::exp(-(y-par[3])*(y-par[3])/(2*par[4]*par[4]));
|
---|
87 | return value;
|
---|
88 | }
|
---|
89 |
|
---|
90 | Bool_t HandleInput()
|
---|
91 | {
|
---|
92 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
---|
93 |
|
---|
94 | while (1)
|
---|
95 | {
|
---|
96 | //
|
---|
97 | // While reading the input process gui events asynchronously
|
---|
98 | //
|
---|
99 | timer.TurnOn();
|
---|
100 | TString input = Getline("Type 'q' to exit, <return> to go on: ");
|
---|
101 | timer.TurnOff();
|
---|
102 |
|
---|
103 | if (input=="q\n")
|
---|
104 | return kFALSE;
|
---|
105 |
|
---|
106 | if (input=="\n")
|
---|
107 | return kTRUE;
|
---|
108 | };
|
---|
109 |
|
---|
110 | return kFALSE;
|
---|
111 | }
|
---|
112 |
|
---|
113 |
|
---|
114 | void pointspreadfunction(TString fname, TString clname = "NULL", Int_t userNumLines = 1000)
|
---|
115 | {
|
---|
116 |
|
---|
117 |
|
---|
118 | Float_t currmean[numPixels];
|
---|
119 | Float_t currsquaremean[numPixels];
|
---|
120 | Float_t currrms[numPixels];
|
---|
121 | Float_t meanofcurrmean = 0;
|
---|
122 |
|
---|
123 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
---|
124 | {
|
---|
125 | currmean[swpix] = 0.;
|
---|
126 | currsquaremean[swpix] = 0.;
|
---|
127 | currrms[swpix] = 0.;
|
---|
128 | }
|
---|
129 |
|
---|
130 | Int_t numLines=0;
|
---|
131 |
|
---|
132 | //containers
|
---|
133 | MGeomCamMagic geomcam;
|
---|
134 | MCameraDC curr;
|
---|
135 | MCameraDC currbis;
|
---|
136 |
|
---|
137 | const Float_t conv4mm2deg = geomcam.GetConvMm2Deg();
|
---|
138 |
|
---|
139 |
|
---|
140 | if (clname != "NULL")
|
---|
141 | {
|
---|
142 |
|
---|
143 | //
|
---|
144 | // First run over continuos light files to have a DC calibration
|
---|
145 | //
|
---|
146 |
|
---|
147 | MParList plist0;
|
---|
148 |
|
---|
149 | MTaskList tlist0;
|
---|
150 | plist0.AddToList(&tlist0);
|
---|
151 |
|
---|
152 | plist0.AddToList(&geomcam);
|
---|
153 | plist0.AddToList(&curr);
|
---|
154 |
|
---|
155 | MReportFileRead read0(clname);
|
---|
156 | read0.SetHasNoHeader();
|
---|
157 | read0.AddToList("MReportCurrents");
|
---|
158 |
|
---|
159 | tlist0.AddToList(&read0);
|
---|
160 |
|
---|
161 | MEvtLoop evtloop0;
|
---|
162 | evtloop0.SetParList(&plist0);
|
---|
163 |
|
---|
164 |
|
---|
165 | if (!evtloop0.PreProcess())
|
---|
166 | return;
|
---|
167 |
|
---|
168 | while (tlist0.Process())
|
---|
169 | {
|
---|
170 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
---|
171 | {
|
---|
172 | meanofcurrmean += curr[swpix];
|
---|
173 | currmean[swpix] += curr[swpix];
|
---|
174 | currsquaremean[swpix] += curr[swpix]*curr[swpix];
|
---|
175 | }
|
---|
176 | numLines++;
|
---|
177 | }
|
---|
178 |
|
---|
179 | evtloop0.PostProcess();
|
---|
180 |
|
---|
181 | meanofcurrmean /= (numLines*numPixels);
|
---|
182 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
---|
183 | {
|
---|
184 |
|
---|
185 | currmean[swpix] /= numLines;
|
---|
186 | currsquaremean[swpix] /= numLines;
|
---|
187 | currrms[swpix] = sqrt(fabs(currsquaremean[swpix] - currmean[swpix]*currmean[swpix]));
|
---|
188 |
|
---|
189 | curr[swpix] = currmean[swpix];
|
---|
190 | currbis[swpix] = currrms[swpix];
|
---|
191 |
|
---|
192 | currmean[swpix] /= meanofcurrmean;
|
---|
193 | currrms[swpix] /= meanofcurrmean;
|
---|
194 |
|
---|
195 | }
|
---|
196 |
|
---|
197 |
|
---|
198 | /* MHCamera display0(geomcam);
|
---|
199 | display0.SetPrettyPalette();
|
---|
200 | display0.Draw();
|
---|
201 |
|
---|
202 | // curr.Print();
|
---|
203 | display0.SetCamContent(currbis);
|
---|
204 | cout << "PSF>> DC mean values drawn" << endl;
|
---|
205 | // Remove the comments if you want to go through the file
|
---|
206 | // event-by-event:
|
---|
207 | if (!HandleInput())
|
---|
208 | break;
|
---|
209 | */
|
---|
210 | }
|
---|
211 | else
|
---|
212 | {
|
---|
213 | // If you don't use the continuous light this is the currrms[] array
|
---|
214 | // is the error associated to the currents TMinuit will use.
|
---|
215 | for (Int_t swpix=0; swpix<numPixels; swpix++)
|
---|
216 | {
|
---|
217 | currmean[swpix] = 1.;
|
---|
218 | currrms[swpix] = 0.2;
|
---|
219 | }
|
---|
220 |
|
---|
221 | }
|
---|
222 |
|
---|
223 | //
|
---|
224 | // Now we can run over the dc data to extract the psf.
|
---|
225 | //
|
---|
226 | const Int_t maxNumLines = 10000;
|
---|
227 |
|
---|
228 | Double_t ux[maxNumLines];
|
---|
229 | Double_t uy[maxNumLines];
|
---|
230 | Double_t sx[maxNumLines];
|
---|
231 | Double_t sy[maxNumLines];
|
---|
232 | Double_t chisqu[maxNumLines];
|
---|
233 | Double_t time[maxNumLines];
|
---|
234 |
|
---|
235 | MParList plist;
|
---|
236 |
|
---|
237 | MGeomCamMagic geomcam;
|
---|
238 | MCameraDC curr;
|
---|
239 | MTaskList tlist;
|
---|
240 |
|
---|
241 | plist.AddToList(&geomcam);
|
---|
242 | plist.AddToList(&curr);
|
---|
243 | plist.AddToList(&tlist);
|
---|
244 |
|
---|
245 | MReportFileRead read(fname);
|
---|
246 | read.SetHasNoHeader();
|
---|
247 | read.AddToList("MReportCurrents");
|
---|
248 |
|
---|
249 | tlist.AddToList(&read);
|
---|
250 |
|
---|
251 | MEvtLoop evtloop;
|
---|
252 | evtloop.SetParList(&plist);
|
---|
253 |
|
---|
254 | if (!evtloop.PreProcess())
|
---|
255 | return;
|
---|
256 |
|
---|
257 | MHCamera display(geomcam);
|
---|
258 | display.SetPrettyPalette();
|
---|
259 | display.Draw();
|
---|
260 | gPad->SetLogy();
|
---|
261 | gPad->cd(1);
|
---|
262 |
|
---|
263 | Double_t amin,edm,errdef;
|
---|
264 | Int_t nvpar,nparx,icstat;
|
---|
265 |
|
---|
266 | Double_t max,maxerror;
|
---|
267 | Double_t xmean,xsigma,ymean,ysigma;
|
---|
268 | Double_t xmeanerror,xsigmaerror,ymeanerror,ysigmaerror;
|
---|
269 |
|
---|
270 | TEllipse ellipse;
|
---|
271 | ellipse.SetFillStyle(4000);
|
---|
272 | ellipse.SetLineWidth(2);
|
---|
273 | ellipse.SetLineColor(2);
|
---|
274 |
|
---|
275 | ellipse.Draw();
|
---|
276 |
|
---|
277 | Int_t nbinsxy = 80;
|
---|
278 | Float_t minxy = -600*conv4mm2deg;
|
---|
279 | Float_t maxxy = 600*conv4mm2deg;
|
---|
280 | Float_t fromdegtobin = (maxxy-minxy)/nbinsxy;
|
---|
281 |
|
---|
282 | TH2D psfhist("psfhist","",nbinsxy,minxy,maxxy,nbinsxy,minxy,maxxy);
|
---|
283 | psfhist->GetXaxis()->SetTitle("[deg]");
|
---|
284 | psfhist->GetYaxis()->SetTitle("[deg]");
|
---|
285 | psfhist->GetZaxis()->SetTitle("DC [uA]");
|
---|
286 |
|
---|
287 | TCanvas *psfcanvas = new TCanvas("psfcanvas","Point Spread Funtion 2D",200,20,900,700);
|
---|
288 |
|
---|
289 |
|
---|
290 | //
|
---|
291 | // Using the first dc measurement we search the pixels which contains the star and define
|
---|
292 | // an area to be fitted by Minuit which is 3 rings of neightbords around the peak of the star.
|
---|
293 | //
|
---|
294 |
|
---|
295 | tlist.Process();
|
---|
296 |
|
---|
297 | Int_t numLines=0;
|
---|
298 | Float_t minDCStar = 6.0;
|
---|
299 |
|
---|
300 | Int_t numPixelsInStar = 0;
|
---|
301 | Float_t maxDC = 0;
|
---|
302 | Int_t swpixelmaxDC;
|
---|
303 |
|
---|
304 | Bool_t isPixelsFittedTmp[numPixels];
|
---|
305 |
|
---|
306 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
307 | isPixelsFittedTmp[swpixel] = kFALSE;
|
---|
308 |
|
---|
309 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
310 | {
|
---|
311 | if(curr[swpixel] > maxDC)
|
---|
312 | {
|
---|
313 | swpixelmaxDC = swpixel;
|
---|
314 | maxDC = curr[swpixelmaxDC];
|
---|
315 | }
|
---|
316 |
|
---|
317 | if(curr[swpixel]>minDCStar)
|
---|
318 | {
|
---|
319 | numPixelsInStar++;
|
---|
320 | isPixelsFitted[swpixel] = kTRUE;
|
---|
321 | }
|
---|
322 | else
|
---|
323 | isPixelsFitted[swpixel] = kFALSE;
|
---|
324 | }
|
---|
325 |
|
---|
326 | if (numPixelsInStar == 0)
|
---|
327 | {
|
---|
328 | cout << "PSF>> Warning: none pixel over minDCStar(" << minDCStar << ')' << endl;
|
---|
329 | return;
|
---|
330 | }
|
---|
331 |
|
---|
332 | //1st neighboor ring
|
---|
333 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
334 | if (isPixelsFitted[swpixel])
|
---|
335 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
---|
336 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
---|
337 |
|
---|
338 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
339 | if (isPixelsFittedTmp[swpixel])
|
---|
340 | isPixelsFitted[swpixel] = kTRUE;
|
---|
341 |
|
---|
342 | //2on neighboor ring
|
---|
343 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
344 | if (isPixelsFitted[swpixel])
|
---|
345 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
---|
346 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
---|
347 |
|
---|
348 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
349 | if (isPixelsFittedTmp[swpixel])
|
---|
350 | isPixelsFitted[swpixel] = kTRUE;
|
---|
351 |
|
---|
352 |
|
---|
353 | //3rt neighboor ring
|
---|
354 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
355 | if (isPixelsFitted[swpixel])
|
---|
356 | for(Int_t next=0; next<geomcam[swpixel].GetNumNeighbors(); next++)
|
---|
357 | isPixelsFittedTmp[geomcam[swpixel].GetNeighbor(next)] = kTRUE;
|
---|
358 |
|
---|
359 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
360 | if (isPixelsFittedTmp[swpixel])
|
---|
361 | isPixelsFitted[swpixel] = kTRUE;
|
---|
362 |
|
---|
363 |
|
---|
364 | for(Int_t swpixel=1; swpixel<numPixels; swpixel++)
|
---|
365 | curr[swpixel] = (Float_t)isPixelsFitted[swpixel];
|
---|
366 |
|
---|
367 | /* MHCamera display0(geomcam);
|
---|
368 | display0.SetPrettyPalette();
|
---|
369 | display0.Draw();
|
---|
370 |
|
---|
371 | display0.SetCamContent(curr);
|
---|
372 | cout << "PSF>> Fitted pixels drawn" << endl;
|
---|
373 | // Remove the comments if you want to go through the file
|
---|
374 | // event-by-event:
|
---|
375 | if (!HandleInput())
|
---|
376 | break;
|
---|
377 | */
|
---|
378 |
|
---|
379 | // Minuit initialization
|
---|
380 |
|
---|
381 | TMinuit *gMinuit = new TMinuit(7); //initialize TMinuit with a maximum of 5 params
|
---|
382 | gMinuit->SetFCN(fcn);
|
---|
383 |
|
---|
384 | Double_t arglist[10];
|
---|
385 | Int_t ierflg = 0;
|
---|
386 |
|
---|
387 | arglist[0] = 1;
|
---|
388 | gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);
|
---|
389 | // arglist[0] = -1;
|
---|
390 | arglist[0] = 0;
|
---|
391 | gMinuit->mnexcm("SET PRI", arglist ,1,ierflg);
|
---|
392 |
|
---|
393 | // Set starting values and step sizes for parameters
|
---|
394 | Double_t vstart[5];
|
---|
395 | Double_t step[5];
|
---|
396 | Double_t lowlimit[5] = {minDCStar, -2., 0.05, -2, 0.05};
|
---|
397 | Double_t uplimit[5] = {30., 2., 1., 2., 1.};
|
---|
398 |
|
---|
399 | vstart[0] = maxDC;
|
---|
400 | vstart[1] = geomcam[swpixelmaxDC].GetX()*conv4mm2deg;
|
---|
401 | vstart[2] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
---|
402 | vstart[3] = geomcam[swpixelmaxDC].GetY()*conv4mm2deg;
|
---|
403 | vstart[4] = 30*TMath::Sqrt(numPixelsInStar/2)*conv4mm2deg;
|
---|
404 |
|
---|
405 | for(Int_t i=0; i<5; i++)
|
---|
406 | {
|
---|
407 | if (vstart[i] != 0)
|
---|
408 | step[i] = TMath::Abs(vstart[i]/sqrt(2));
|
---|
409 | else
|
---|
410 | step[i] = uplimit[i]/2;
|
---|
411 | }
|
---|
412 |
|
---|
413 | gMinuit->mnparm(0, "max", vstart[0], step[0], lowlimit[0], uplimit[0],ierflg);
|
---|
414 | gMinuit->mnparm(1, "xmean", vstart[1], step[1], lowlimit[1], uplimit[1],ierflg);
|
---|
415 | gMinuit->mnparm(2, "xsigma", vstart[2], step[2], lowlimit[2], uplimit[2],ierflg);
|
---|
416 | gMinuit->mnparm(3, "ymean", vstart[3], step[3], lowlimit[3], uplimit[3],ierflg);
|
---|
417 | gMinuit->mnparm(4, "ysigma", vstart[4], step[4], lowlimit[4], uplimit[4],ierflg);
|
---|
418 |
|
---|
419 | while (tlist.Process() && numLines < maxNumLines)
|
---|
420 | {
|
---|
421 |
|
---|
422 | for (Int_t swpixel=1; swpixel<577; swpixel++)
|
---|
423 | {
|
---|
424 |
|
---|
425 | x[swpixel] = geomcam[swpixel].GetX()*conv4mm2deg;
|
---|
426 | y[swpixel] = geomcam[swpixel].GetY()*conv4mm2deg;
|
---|
427 | z[swpixel] = curr[swpixel]/currmean[swpixel];
|
---|
428 | errorz[swpixel] = TMath::Sqrt((curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))*(curr[swpixel]*currrms[swpixel]/(currmean[swpixel]*currmean[swpixel]))+(0.1)/currmean[swpixel]*(0.1)/currmean[swpixel]);
|
---|
429 |
|
---|
430 |
|
---|
431 | psfhist->SetBinContent((Int_t)((x[swpixel]+600*conv4mm2deg)/fromdegtobin),(Int_t)((y[swpixel]+600*conv4mm2deg)/fromdegtobin),z[swpixel]);
|
---|
432 | }
|
---|
433 |
|
---|
434 | psfcanvas->cd(1);
|
---|
435 | psfhist->Draw("lego2");
|
---|
436 |
|
---|
437 | // Now ready for minimization step
|
---|
438 | arglist[0] = 500;
|
---|
439 | arglist[1] = 1.;
|
---|
440 | gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
|
---|
441 |
|
---|
442 | // Print results
|
---|
443 | /* gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
|
---|
444 | gMinuit->mnprin(3,amin);
|
---|
445 | */
|
---|
446 | gMinuit->GetParameter(0,max,maxerror);
|
---|
447 | gMinuit->GetParameter(1,xmean,xmeanerror);
|
---|
448 | gMinuit->GetParameter(2,xsigma,xsigmaerror);
|
---|
449 | gMinuit->GetParameter(3,ymean,ymeanerror);
|
---|
450 | gMinuit->GetParameter(4,ysigma,ysigmaerror);
|
---|
451 |
|
---|
452 | /* cout << endl;
|
---|
453 | cout << "numLine " << numLines << endl;
|
---|
454 | cout << "max \t" << max << " +- " << maxerror << endl;
|
---|
455 | cout << "xmean \t" << xmean << " +- " << xmeanerror << endl;
|
---|
456 | cout << "xsigma \t" << TMath::Abs(xsigma) << " +- " << xsigmaerror << endl;
|
---|
457 | cout << "ymean \t" << ymean << " +- " << ymeanerror << endl;
|
---|
458 | cout << "ysigma \t" << TMath::Abs(ysigma) << " +- " << ysigmaerror << endl;
|
---|
459 | cout << "chisquare/ndof \t" << chisquare/(nPixelsFitted-5) << endl;
|
---|
460 | */
|
---|
461 |
|
---|
462 | chisqu[numLines]=chisquare/(nPixelsFitted-5);
|
---|
463 |
|
---|
464 | if(chisqu[numLines]<100.)
|
---|
465 | {
|
---|
466 | ux[numLines]=xmean;
|
---|
467 | uy[numLines]=ymean;
|
---|
468 | sx[numLines]=TMath::Abs(xsigma);
|
---|
469 | sy[numLines]=TMath::Abs(ysigma);
|
---|
470 | time[numLines]=numLines;
|
---|
471 |
|
---|
472 | display.SetCamContent(curr);
|
---|
473 | gPad->cd(1);
|
---|
474 | ellipse.SetX1(xmean/conv4mm2deg);
|
---|
475 | ellipse.SetY1(ymean/conv4mm2deg);
|
---|
476 | ellipse.SetR1(TMath::Abs(xsigma)/conv4mm2deg);
|
---|
477 | ellipse.SetR2(TMath::Abs(ysigma)/conv4mm2deg);
|
---|
478 |
|
---|
479 | gPad->Modified();
|
---|
480 | gPad->Update();
|
---|
481 |
|
---|
482 | // Remove the comments if you want to go through the file
|
---|
483 | //event-by-event:
|
---|
484 | if (userNumLines>0)
|
---|
485 | {
|
---|
486 | if (numLines>userNumLines)
|
---|
487 | break;
|
---|
488 | }
|
---|
489 | else
|
---|
490 | {
|
---|
491 | if (!HandleInput())
|
---|
492 | break;
|
---|
493 | }
|
---|
494 | numLines++;
|
---|
495 | }
|
---|
496 | }
|
---|
497 |
|
---|
498 |
|
---|
499 | evtloop.PostProcess();
|
---|
500 |
|
---|
501 | //
|
---|
502 | // Draw the ditributions of the sigmas the point spread function
|
---|
503 | //
|
---|
504 |
|
---|
505 | cout<<"PSF>> Number of lines "<<numLines<<endl;
|
---|
506 |
|
---|
507 | gROOT->Reset();
|
---|
508 | gStyle->SetCanvasColor(0);
|
---|
509 | gStyle->SetCanvasBorderMode(0);
|
---|
510 | gStyle->SetPadBorderMode(0);
|
---|
511 | gStyle->SetFrameBorderMode(0);
|
---|
512 | gStyle->SetOptStat(00000000);
|
---|
513 |
|
---|
514 | //
|
---|
515 | // Find in the file name the date, run and project name to put it in the title
|
---|
516 | //
|
---|
517 |
|
---|
518 | Size_t pos = fname.Last('/');
|
---|
519 | TString iRun = TString(fname(pos+24,5));
|
---|
520 | TString iYear = TString(fname(pos+4,4));
|
---|
521 | TString iMonth = TString(fname(pos+9,2));
|
---|
522 | TString iDay = TString(fname(pos+12,2));
|
---|
523 |
|
---|
524 | TString iHour = TString(fname(pos+15,2));
|
---|
525 | TString iMin = TString(fname(pos+18,2));
|
---|
526 | TString iSec = TString(fname(pos+21,2));
|
---|
527 |
|
---|
528 | Size_t poslast = fname.Last('.');
|
---|
529 | Size_t posfirst = poslast-1;
|
---|
530 | while (fname[posfirst] != '_')
|
---|
531 | posfirst--;
|
---|
532 |
|
---|
533 | TString iSource = TString(fname(posfirst+1,poslast-posfirst-1));
|
---|
534 |
|
---|
535 |
|
---|
536 | char str[100];
|
---|
537 |
|
---|
538 | // sprintf(str,"Date %s/%s/%s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iRun.Data(),iSource.Data());
|
---|
539 | sprintf(str,"Date %s/%s/%s StartTime %s:%s:%s Run %s Source %s",iYear.Data(),iMonth.Data(),iDay.Data(),iHour.Data(),iMin.Data(),iSec.Data(),iRun.Data(),iSource.Data());
|
---|
540 |
|
---|
541 | c1 = new TCanvas("c1",str,0,0,1200,850);
|
---|
542 | // c1 = new TCanvas("c1","Time evolution & distributions",0,0,1200,850);
|
---|
543 | c1->Divide(3,2);
|
---|
544 |
|
---|
545 | c1->cd(1);
|
---|
546 |
|
---|
547 | TMath math;
|
---|
548 |
|
---|
549 | Double_t minmeanx, maxmeanx;
|
---|
550 | minmeanx = ux[math.LocMin(numLines,ux)];
|
---|
551 | maxmeanx = ux[math.LocMax(numLines,ux)];
|
---|
552 |
|
---|
553 | Double_t minmeany, maxmeany;
|
---|
554 | minmeany = uy[math.LocMin(numLines,uy)];
|
---|
555 | maxmeany = uy[math.LocMax(numLines,uy)];
|
---|
556 |
|
---|
557 | Double_t minmean, maxmean;
|
---|
558 | minmean = math.Min(minmeanx,minmeany);
|
---|
559 | maxmean = math.Max(maxmeanx,maxmeany);
|
---|
560 |
|
---|
561 | Double_t diff;
|
---|
562 | diff = maxmean - minmean;
|
---|
563 | diff = 0.1*diff;
|
---|
564 | minmean = minmean - diff;
|
---|
565 | maxmean = maxmean + diff;
|
---|
566 |
|
---|
567 | Double_t mintime, maxtime;
|
---|
568 | mintime = time[math.LocMin(numLines,time)];
|
---|
569 | maxtime = time[math.LocMax(numLines,time)];
|
---|
570 |
|
---|
571 | TH2D *h1 = new TH2D("h1","",1,mintime-1,maxtime+1,1,minmean,maxmean);
|
---|
572 | h1->GetXaxis()->SetTitle("Event number");
|
---|
573 | h1->GetYaxis()->SetTitle("mean position (deg)");
|
---|
574 | h1->Draw();
|
---|
575 |
|
---|
576 | TGraph *grtimeevolmeanx = new TGraph(numLines,time,ux);
|
---|
577 | grtimeevolmeanx->SetMarkerColor(3);
|
---|
578 | grtimeevolmeanx->SetMarkerStyle(20);
|
---|
579 | grtimeevolmeanx->SetMarkerSize (0.4);
|
---|
580 | grtimeevolmeanx->Draw("P");
|
---|
581 |
|
---|
582 | TGraph *grtimeevolmeany = new TGraph(numLines,time,uy);
|
---|
583 | grtimeevolmeany->SetMarkerColor(6);
|
---|
584 | grtimeevolmeany->SetMarkerStyle(24);
|
---|
585 | grtimeevolmeany->SetMarkerSize (0.4);
|
---|
586 | grtimeevolmeany->Draw("P");
|
---|
587 |
|
---|
588 | legmeanxy = new TLegend(0.8,0.85,0.95,0.95);
|
---|
589 | legmeanxy.SetTextSize(0.03);
|
---|
590 | legmeanxy.AddEntry(grtimeevolmeanx,"mean x","P");
|
---|
591 | legmeanxy.AddEntry(grtimeevolmeany,"mean y","P");
|
---|
592 | legmeanxy.Draw();
|
---|
593 |
|
---|
594 | c1->cd(2);
|
---|
595 |
|
---|
596 | TMath math;
|
---|
597 |
|
---|
598 | Double_t minsigmax, maxsigmax;
|
---|
599 | minsigmax = sx[math.LocMin(numLines,sx)];
|
---|
600 | maxsigmax = sx[math.LocMax(numLines,sx)];
|
---|
601 |
|
---|
602 | Double_t minsigmay, maxsigmay;
|
---|
603 | minsigmay = sy[math.LocMin(numLines,sy)];
|
---|
604 | maxsigmay = sy[math.LocMax(numLines,sy)];
|
---|
605 |
|
---|
606 | Double_t minsigma, maxsigma;
|
---|
607 | minsigma = math.Min(minsigmax,minsigmay);
|
---|
608 | maxsigma = math.Max(maxsigmax,maxsigmay);
|
---|
609 |
|
---|
610 | diff = maxsigma - minsigma;
|
---|
611 | diff = 0.1*diff;
|
---|
612 | minsigma = minsigma - diff;
|
---|
613 | maxsigma = maxsigma + diff;
|
---|
614 |
|
---|
615 | TH2D *h2 = new TH2D("h2","",1,mintime-1,maxtime+1,1,minsigma,maxsigma);
|
---|
616 | h2->GetXaxis()->SetTitle("Event number");
|
---|
617 | h2->GetYaxis()->SetTitle("PSF Rms (deg)");
|
---|
618 | h2->Draw();
|
---|
619 |
|
---|
620 | TGraph* grtimeevolsigmax= new TGraph(numLines,time,sx);
|
---|
621 | grtimeevolsigmax->SetMarkerColor(3);
|
---|
622 | grtimeevolsigmax->SetMarkerStyle(20);
|
---|
623 | grtimeevolsigmax->SetMarkerSize (0.4);
|
---|
624 | grtimeevolsigmax->Draw("P");
|
---|
625 |
|
---|
626 | TGraph* grtimeevolsigmay= new TGraph(numLines,time,sy);
|
---|
627 | grtimeevolsigmay->SetMarkerColor(6);
|
---|
628 | grtimeevolsigmay->SetMarkerStyle(24);
|
---|
629 | grtimeevolsigmay->SetMarkerSize (0.4);
|
---|
630 | grtimeevolsigmay->Draw("P");
|
---|
631 |
|
---|
632 | legsigmaxy = new TLegend(0.8,0.85,0.95,0.95);
|
---|
633 | legsigmaxy.SetTextSize(0.03);
|
---|
634 | legsigmaxy.AddEntry(grtimeevolsigmax,"sigma x","P");
|
---|
635 | legsigmaxy.AddEntry(grtimeevolsigmay,"sigma y","P");
|
---|
636 | legsigmaxy.Draw();
|
---|
637 |
|
---|
638 | c1->cd(3);
|
---|
639 |
|
---|
640 | Double_t minchisqu, maxchisqu;
|
---|
641 |
|
---|
642 | minchisqu = chisqu[math.LocMin(numLines,chisqu)];
|
---|
643 | maxchisqu = chisqu[math.LocMax(numLines,chisqu)];
|
---|
644 |
|
---|
645 | diff = maxchisqu - minchisqu;
|
---|
646 | diff = 0.1*diff;
|
---|
647 | minchisqu = minchisqu - diff;
|
---|
648 | maxchisqu = maxchisqu + diff;
|
---|
649 |
|
---|
650 | TH2D *h3 = new TH2D("h3","",1,mintime-1,maxtime+1,1,minchisqu,maxchisqu);
|
---|
651 | h3->GetXaxis()->SetTitle("Event number");
|
---|
652 | h3->Draw();
|
---|
653 |
|
---|
654 | TGraph * grtimeevolchisqu = new TGraph(numLines,time,chisqu);
|
---|
655 | grtimeevolchisqu->SetMarkerColor(215);
|
---|
656 | grtimeevolchisqu->SetMarkerStyle(20);
|
---|
657 | grtimeevolchisqu->SetMarkerSize(0.4);
|
---|
658 | grtimeevolchisqu->Draw("P");
|
---|
659 |
|
---|
660 | legchisqu = new TLegend(0.55,0.90,0.95,0.95);
|
---|
661 | legchisqu.SetTextSize(0.03);
|
---|
662 | legchisqu.AddEntry(grtimeevolchisqu,"chi square / ndof","P");
|
---|
663 | legchisqu.Draw();
|
---|
664 |
|
---|
665 | //***************************************
|
---|
666 |
|
---|
667 | const Int_t nbins = 100;
|
---|
668 |
|
---|
669 | TH1D *xsigmahist = new TH1D("xsigmahist","",nbins,minsigma,maxsigma);
|
---|
670 | TH1D *ysigmahist = new TH1D("ysigmahist","",nbins,minsigma,maxsigma);
|
---|
671 | TH1D *xmeanhist = new TH1D("xmeanhist","",nbins,minmean,maxmean);
|
---|
672 | TH1D *ymeanhist = new TH1D("ymeanhist","",nbins,minmean,maxmean);
|
---|
673 | TH1D *chisquhist = new TH1D("chisquhist","",nbins,minchisqu,maxchisqu);
|
---|
674 |
|
---|
675 | for (Int_t i=0; i<numLines; i++)
|
---|
676 | {
|
---|
677 | xsigmahist->Fill(TMath::Abs(sx[i]));
|
---|
678 | ysigmahist->Fill(TMath::Abs(sy[i]));
|
---|
679 | xmeanhist->Fill(ux[i]);
|
---|
680 | ymeanhist->Fill(uy[i]);
|
---|
681 | chisquhist->Fill(chisqu[i]);
|
---|
682 |
|
---|
683 | }
|
---|
684 |
|
---|
685 | c1->cd(5);
|
---|
686 |
|
---|
687 | TMath math;
|
---|
688 | Double_t maxsigma;
|
---|
689 | Int_t binmaxx, binmaxy;
|
---|
690 | xsigmahist->SetLineColor(3);
|
---|
691 | xsigmahist->SetLineWidth(2);
|
---|
692 | xsigmahist->SetXTitle("RMS [deg]");
|
---|
693 | binmaxx = xsigmahist->GetMaximumBin();
|
---|
694 | binmaxx = xsigmahist->GetBinContent(binmaxx);
|
---|
695 |
|
---|
696 | ysigmahist->SetLineColor(6);
|
---|
697 | ysigmahist->SetLineWidth(2);
|
---|
698 | binmaxy = ysigmahist->GetMaximumBin();
|
---|
699 | binmaxy = ysigmahist->GetBinContent(binmaxy);
|
---|
700 |
|
---|
701 | maxsigma = math.Max(binmaxx,binmaxy);
|
---|
702 | maxsigma += 1;
|
---|
703 |
|
---|
704 | xsigmahist->SetMaximum(maxsigma);
|
---|
705 | ysigmahist->SetMaximum(maxsigma);
|
---|
706 | xsigmahist->DrawCopy();
|
---|
707 | ysigmahist->DrawCopy("Same");
|
---|
708 |
|
---|
709 | TLegend *legendsigma = new TLegend(.3,.8,.95,.95);
|
---|
710 | legendsigma->SetTextSize(0.03);
|
---|
711 | char xsigmatitle[100];
|
---|
712 | char ysigmatitle[100];
|
---|
713 | sprintf(xsigmatitle,"PSF Rms on X axis -- %5.2f +/- %5.2f deg",xsigmahist->GetMean(),xsigmahist->GetRMS());
|
---|
714 | sprintf(ysigmatitle,"PSF Rms on Y axis -- %5.2f +/- %5.2f deg",ysigmahist->GetMean(),ysigmahist->GetRMS());
|
---|
715 | legendsigma->AddEntry(xsigmahist,xsigmatitle,"F");
|
---|
716 | legendsigma->AddEntry(ysigmahist,ysigmatitle,"F");
|
---|
717 | legendsigma->Draw();
|
---|
718 |
|
---|
719 | c1->cd(4);
|
---|
720 |
|
---|
721 | Double_t maxmean;
|
---|
722 |
|
---|
723 | xmeanhist->SetLineColor(3);
|
---|
724 | xmeanhist->SetLineWidth(2);
|
---|
725 | xmeanhist->SetXTitle("mean [deg]");
|
---|
726 | binmaxx = xmeanhist->GetMaximumBin();
|
---|
727 | binmaxx = xmeanhist->GetBinContent(binmaxx);
|
---|
728 |
|
---|
729 | ymeanhist->SetLineColor(6);
|
---|
730 | ymeanhist->SetLineWidth(2);
|
---|
731 | binmaxy = ymeanhist->GetMaximumBin();
|
---|
732 | binmaxy = ymeanhist->GetBinContent(binmaxy);
|
---|
733 |
|
---|
734 | maxmean = math.Max(binmaxx,binmaxy);
|
---|
735 | maxmean += 1;
|
---|
736 |
|
---|
737 | xmeanhist->SetMaximum(maxmean);
|
---|
738 | ymeanhist->SetMaximum(maxmean);
|
---|
739 | xmeanhist->DrawCopy();
|
---|
740 | ymeanhist->DrawCopy("Same");
|
---|
741 |
|
---|
742 | TLegend *legendmean = new TLegend(.35,.8,.95,.95);
|
---|
743 | legendmean->SetTextSize(0.03);
|
---|
744 | char xmeantitle[100];
|
---|
745 | char ymeantitle[100];
|
---|
746 | sprintf(xmeantitle,"mean on X axis -- %5.2f +/- %5.2f deg",xmeanhist->GetMean(),xmeanhist->GetRMS());
|
---|
747 | sprintf(ymeantitle,"mean on Y axis -- %5.2f +/- %5.2f deg",ymeanhist->GetMean(),ymeanhist->GetRMS());
|
---|
748 | legendmean->AddEntry(xmeanhist,xmeantitle,"F");
|
---|
749 | legendmean->AddEntry(ymeanhist,ymeantitle,"F");
|
---|
750 | legendmean->Draw();
|
---|
751 |
|
---|
752 | //meancanvas->Modified();
|
---|
753 | //meancanvas->Update();
|
---|
754 |
|
---|
755 | c1->cd(6);
|
---|
756 |
|
---|
757 | chisquhist->SetLineColor(215);
|
---|
758 | chisquhist->SetLineWidth(2);
|
---|
759 | chisquhist->SetXTitle("chi square / ndof");
|
---|
760 | TAxis * axis = chisquhist->GetXaxis();
|
---|
761 | axis->SetLabelSize(0.025);
|
---|
762 | chisquhist->DrawCopy();
|
---|
763 |
|
---|
764 | TLegend *legendchisqu = new TLegend(.4,.85,.95,.95);
|
---|
765 | legendchisqu->SetTextSize(0.03);
|
---|
766 | char chisqutitle[100];
|
---|
767 | sprintf(chisqutitle,"chi square / ndof -- %5.2f +/- %5.2f ",chisquhist->GetMean(),chisquhist->GetRMS());
|
---|
768 | legendchisqu->AddEntry(chisquhist,chisqutitle,"F");
|
---|
769 | legendchisqu->Draw();
|
---|
770 |
|
---|
771 |
|
---|
772 | return;
|
---|
773 |
|
---|
774 | }
|
---|
775 |
|
---|
776 |
|
---|
777 |
|
---|
778 |
|
---|
779 |
|
---|
780 |
|
---|
781 |
|
---|
782 |
|
---|
783 |
|
---|
784 |
|
---|
785 |
|
---|
786 |
|
---|
787 |
|
---|
788 |
|
---|
789 |
|
---|
790 |
|
---|
791 |
|
---|
792 |
|
---|
793 |
|
---|
794 |
|
---|
795 |
|
---|
796 |
|
---|
797 |
|
---|
798 |
|
---|
799 |
|
---|
800 |
|
---|
801 |
|
---|
802 |
|
---|
803 |
|
---|
804 |
|
---|
805 |
|
---|
806 |
|
---|
807 |
|
---|
808 |
|
---|
809 |
|
---|
810 |
|
---|
811 |
|
---|
812 |
|
---|
813 |
|
---|
814 |
|
---|
815 |
|
---|
816 |
|
---|
817 |
|
---|
818 |
|
---|
819 |
|
---|
820 |
|
---|
821 |
|
---|
822 |
|
---|
823 |
|
---|
824 |
|
---|
825 |
|
---|
826 |
|
---|
827 |
|
---|
828 |
|
---|
829 |
|
---|
830 |
|
---|
831 |
|
---|
832 |
|
---|
833 |
|
---|
834 |
|
---|
835 |
|
---|
836 |
|
---|
837 |
|
---|
838 |
|
---|
839 |
|
---|
840 |
|
---|
841 |
|
---|
842 |
|
---|
843 |
|
---|
844 |
|
---|
845 |
|
---|
846 |
|
---|
847 |
|
---|
848 |
|
---|
849 |
|
---|
850 |
|
---|
851 |
|
---|
852 |
|
---|
853 |
|
---|