source: branches/Corsika7500Compatibility/mtemp/mifae/library/MHPSFFromStars.cc@ 18752

Last change on this file since 18752 was 4295, checked in by jlopez, 21 years ago
*** empty log message ***
File size: 14.5 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): Javier López 05/2004 <mailto:jlopezs@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24/////////////////////////////////////////////////////////////////////////////
25//
26// MHPSFFromStars
27//
28/////////////////////////////////////////////////////////////////////////////
29#include "MHPSFFromStars.h"
30
31#include <TVirtualPad.h>
32#include <TCanvas.h>
33#include <TPad.h>
34#include <TText.h>
35#include <TPaveText.h>
36#include <TF1.h>
37
38#include "MLog.h"
39#include "MLogManip.h"
40
41#include "MGeomCam.h"
42#include "MGeomPix.h"
43#include "MCameraDC.h"
44
45#include "MStarLocalCam.h"
46#include "MStarLocalPos.h"
47
48#include "MParList.h"
49
50ClassImp(MHPSFFromStars);
51
52using namespace std;
53
54// --------------------------------------------------------------------------
55//
56// Default Constructor.
57//
58MHPSFFromStars::MHPSFFromStars(Int_t starpos, const char *name, const char *title)
59 : fGeom(NULL), fCurr(NULL), fSelectedStarPos(starpos)
60{
61 fName = name ? name : "MHPSFFromStars";
62 fTitle = title ? title : "Class to fill the Point Spread Function histograms";
63
64 fNumEvents=0;
65
66 fHistMeanX.SetName("MeanX");
67 fHistMeanX.SetTitle("Mean X of the Fit");
68 fHistMeanX.SetDirectory(NULL);
69 fHistMeanX.UseCurrentStyle();
70 fHistMeanX.SetXTitle("Mean X [\\circ]");
71 fHistMeanX.SetYTitle("Counts [#]");
72 fHistMeanX.SetBins(800,-2.0,2.0);
73
74 fHistMeanY.SetName("MeanY");
75 fHistMeanY.SetTitle("Mean Y of the Fit");
76 fHistMeanY.SetDirectory(NULL);
77 fHistMeanY.UseCurrentStyle();
78 fHistMeanY.SetXTitle("Mean Y [\\circ]");
79 fHistMeanY.SetYTitle("Counts [#]");
80 fHistMeanY.SetBins(800,-2.0,2.0);
81
82 fHistSigmaMinor.SetName("SigmaMinor");
83 fHistSigmaMinor.SetTitle("Sigma X of the Fit");
84 fHistSigmaMinor.SetDirectory(NULL);
85 fHistSigmaMinor.UseCurrentStyle();
86 fHistSigmaMinor.SetXTitle("Sigma X [\\circ]");
87 fHistSigmaMinor.SetYTitle("Counts [#]");
88 fHistSigmaMinor.SetBins(160,0,0.8);
89
90 fHistSigmaMajor.SetName("SigmaMajor");
91 fHistSigmaMajor.SetTitle("Sigma Y of the Fit");
92 fHistSigmaMajor.SetDirectory(NULL);
93 fHistSigmaMajor.UseCurrentStyle();
94 fHistSigmaMajor.SetXTitle("Sigma Y [\\circ]");
95 fHistSigmaMajor.SetYTitle("Counts [#]");
96 fHistSigmaMajor.SetBins(160,0,0.8);
97
98 fProjectionX.SetName("ProjetionX");
99 fProjectionX.SetTitle("Star projection in X axis");
100 fProjectionX.SetDirectory(NULL);
101 fProjectionX.UseCurrentStyle();
102 fProjectionX.SetXTitle("X [\\circ]");
103 fProjectionX.SetYTitle("DC [uA]");
104 fProjectionX.SetBins(16,-0.8,0.8);
105
106 fProjectionY.SetName("ProjetionY");
107 fProjectionY.SetTitle("Star projection in Y axis");
108 fProjectionY.SetDirectory(NULL);
109 fProjectionY.UseCurrentStyle();
110 fProjectionY.SetXTitle("Y [\\circ]");
111 fProjectionY.SetYTitle("DC [uA]");
112 fProjectionY.SetBins(16,-0.8,0.8);
113
114 fProfile.SetName("Profile");
115 fProfile.SetTitle("Star profile");
116 fProfile.SetDirectory(NULL);
117 fProfile.UseCurrentStyle();
118 fProfile.SetXTitle("X [\\circ]");
119 fProfile.SetYTitle("X [\\circ]");
120 fProfile.SetZTitle("DC [uA]");
121 fProfile.SetBins(16,-0.8,0.8,16,-0.8,0.8);
122
123 fvsTimeMeanX.Set(0);
124 fvsTimeMeanY.Set(0);
125 fvsTimeSigmaMinor.Set(0);
126 fvsTimeSigmaMajor.Set(0);
127
128 // fTime.Set(0);
129 fEvent.Set(0);
130
131}
132
133MHPSFFromStars::~MHPSFFromStars()
134{
135 delete fGraphMeanX;
136 delete fGraphMeanY;
137 delete fGraphSigmaMinor;
138 delete fGraphSigmaMajor;
139}
140
141
142// --------------------------------------------------------------------------
143//
144// TObject *MHPSFFromStars::lone(const char *) const
145// {
146// MHPSFFromStars *hpsf = new MHPSFFromStars();
147
148// return hpsf;
149// }
150
151// --------------------------------------------------------------------------
152//
153// Gets the pointers to:
154//
155Bool_t MHPSFFromStars::SetupFill(const MParList *pList)
156{
157
158 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
159 if (!fGeom)
160 {
161 *fLog << err << GetDescriptor()
162 << ": MGeomCam not found... aborting." << endl;
163 return kFALSE;
164 }
165
166 TArrayC PixelsUsed;
167 PixelsUsed.Set(577);
168 PixelsUsed.Reset((Char_t)kTRUE);
169 fCamera.SetGeometry(*fGeom,"StarsDisplay","StarsDisplay");
170 fCamera.SetUsed(PixelsUsed);
171
172 fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));
173
174 if (!fCurr)
175 {
176 *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
177 return kFALSE;
178 }
179
180 return kTRUE;
181}
182
183//--------------------------------------------------------------------------------
184//
185//
186Bool_t MHPSFFromStars::Fill(const MParContainer *par, const Stat_t w)
187{
188 const UInt_t numPixels = fGeom->GetNumPixels();
189
190 MStarLocalCam* stars = (MStarLocalCam*)par;
191 if (!stars)
192 {
193 *fLog << err << "No argument in " << GetName() << "::Fill... abort." << endl;
194 return kFALSE;
195 }
196
197 //
198 MStarLocalPos* star = SelectStar(stars);
199 if (star == NULL)
200 return kCONTINUE;
201
202
203 //Check good conditions of the fitting procedure
204 Float_t expectedStarXPos = star->GetMeanX();
205 Float_t expectedStarYPos = star->GetMeanY();
206
207 if ( star->GetChiSquareNdof() == 0
208 || star->GetChiSquareNdof()>5.)
209// || TMath::Sqrt(expectedStarXPos*expectedStarXPos+expectedStarYPos*expectedStarYPos)>150.)
210 {
211 *fLog << warn << "Star selected were bad reconstructed" << endl;
212 return kCONTINUE;
213 }
214
215
216
217 //Fill Histograms
218 fHistMeanX.Fill(star->GetMeanX()*fGeom->GetConvMm2Deg());
219 fHistMeanY.Fill(star->GetMeanY()*fGeom->GetConvMm2Deg());
220 if (star->GetSigmaMinorAxis() < 0)
221 *fLog << err << GetName() << " Sigma of the fit " << star->GetSigmaMinorAxis() << " negative" << endl;
222 if (star->GetSigmaMajorAxis() < 0)
223 *fLog << err << GetName() << " Sigma of the fit " << star->GetSigmaMajorAxis() << " negative" << endl;
224 fHistSigmaMinor.Fill(star->GetSigmaMinorAxis()*fGeom->GetConvMm2Deg());
225 fHistSigmaMajor.Fill(star->GetSigmaMajorAxis()*fGeom->GetConvMm2Deg());
226
227 fCamera.AddCamContent(*fCurr);
228
229 //Selected star X(andY) projections
230 Float_t ringofinterest = 4*(star->GetSigmaMajorAxis()>star->GetSigmaMinorAxis()?star->GetSigmaMajorAxis():star->GetSigmaMinorAxis());
231
232 // First define a area of interest around the expected position of the star
233 for (UInt_t pix=1; pix<numPixels; pix++)
234 {
235 Float_t pixelXCenter = (*fGeom)[pix].GetX();
236 Float_t pixelYCenter = (*fGeom)[pix].GetY();
237 Float_t dist = TMath::Sqrt((pixelXCenter-expectedStarXPos)*(pixelXCenter-expectedStarXPos)+
238 (pixelYCenter-expectedStarYPos)*(pixelYCenter-expectedStarYPos));
239
240 if ((dist < ringofinterest))
241 {
242 fProjectionX.Fill((pixelXCenter-expectedStarXPos)*fGeom->GetConvMm2Deg(), (*fCurr)[pix]);
243 fProjectionY.Fill((pixelYCenter-expectedStarYPos)*fGeom->GetConvMm2Deg(), (*fCurr)[pix]);
244 //FIXME fProfile.Fill((pixelXCenter-expectedStarXPos)*fGeom->GetConvMm2Deg(),(pixelYCenter-expectedStarYPos)*fGeom->GetConvMm2Deg(), (*fCurr)[pix]);
245 }
246 }
247
248
249 //
250 fvsTimeMeanX.Set(fNumEvents+1);
251 fvsTimeMeanY.Set(fNumEvents+1);
252 fvsTimeSigmaMinor.Set(fNumEvents+1);
253 fvsTimeSigmaMajor.Set(fNumEvents+1);
254
255 // fTime.Set(fNumEvents+1);
256 fEvent.Set(fNumEvents+1);
257
258 fvsTimeMeanX[fNumEvents] = star->GetMeanX()*fGeom->GetConvMm2Deg();
259 fvsTimeMeanY[fNumEvents] = star->GetMeanY()*fGeom->GetConvMm2Deg();
260 fvsTimeSigmaMinor[fNumEvents] = star->GetSigmaMinorAxis()*fGeom->GetConvMm2Deg();
261 fvsTimeSigmaMajor[fNumEvents] = star->GetSigmaMajorAxis()*fGeom->GetConvMm2Deg();
262
263 // fTime[fNumEvents] = ?;
264 fEvent[fNumEvents] = fNumEvents;
265
266 fNumEvents++;
267
268// fGraphMeanX.SetPoint(fNumEvents,fNumEvents,star->GetMeanX()*fGeom->GetConvMm2Deg());
269// fGraphMeanY.SetPoint(fNumEvents,fNumEvents,star->GetMeanY()*fGeom->GetConvMm2Deg());
270// fGraphSigmaMinor.SetPoint(fNumEvents,fNumEvents,star->GetSigmaMinorAxis()*fGeom->GetConvMm2Deg());
271// fGraphSigmaMajor.SetPoint(fNumEvents,fNumEvents,star->GetSigmaMajorAxis()*fGeom->GetConvMm2Deg());
272
273 return kTRUE;
274}
275
276// --------------------------------------------------------------------------
277//
278//
279Bool_t MHPSFFromStars::Finalize()
280{
281 *fLog << inf << GetName() << " Statistics: " << fNumEvents << " events" << endl;
282 //Fit the profile plot with a gaussian(+baseline)
283
284 return kTRUE;
285}
286
287MStarLocalPos* MHPSFFromStars::SelectStar(MStarLocalCam* stars)
288{
289 TIter Next(stars->GetList());
290 MStarLocalPos* star = NULL;
291 MStarLocalPos* selectedStar = NULL;
292
293 Float_t dist;
294
295 if (fSelectedStarPos >= 0)
296 for (Int_t i=0; i<=fSelectedStarPos; i++)
297 star = (MStarLocalPos*)Next();
298 else //Look for the closer star to the center
299 {
300 Float_t mindist = 600; //mm
301 while ((star=(MStarLocalPos*)Next()))
302 {
303 dist = TMath::Sqrt(star->GetMeanX()*star->GetMeanX() +
304 star->GetMeanY()*star->GetMeanY());
305 if (dist < mindist)
306 {
307 selectedStar = star;
308 mindist = dist;
309 }
310 }
311
312 star = selectedStar;
313 if (star == NULL)
314 *fLog << warn << "Not found star closer to center" << endl;
315
316 }
317
318 return star;
319}
320
321static Double_t fitfunc(Double_t *x, Double_t *par)
322{
323 Double_t fitval =
324 par[0] + par[1]*TMath::Exp(-0.5*(x[0]-par[2])*(x[0]-par[2])/(par[3]*par[3]));
325
326 return fitval;
327}
328
329// -----------------------------------------------------------------------------
330//
331// Default draw:
332//
333//
334void MHPSFFromStars::Draw(const Option_t *opt)
335{
336
337 TString option(opt);
338 option.ToLower();
339
340 TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);
341 pad->SetBorderMode(0);
342
343 if (fNumEvents==0)
344 {
345 *fLog << err << GetName() << " no entries to be drawn" << endl;
346 return;
347 }
348
349
350
351 if (option.Contains("mean"))
352 {
353 fGraphMeanX = new TGraph((Int_t)fNumEvents,fEvent.GetArray(),fvsTimeMeanX.GetArray());
354 fGraphMeanX->SetName("GraphMeanX");
355 fGraphMeanX->SetTitle("Star X position v.s. event ");
356 fGraphMeanX->UseCurrentStyle();
357 fGraphMeanX->GetXaxis()->SetTitle("Event [#]");
358 fGraphMeanX->GetYaxis()->SetTitle("X [\\circ]");
359
360 fGraphMeanY = new TGraph((Int_t)fNumEvents,fEvent.GetArray(),fvsTimeMeanY.GetArray());
361 fGraphMeanY->SetName("GraphMeanY");
362 fGraphMeanY->SetTitle("Star Y position v.s. event ");
363 fGraphMeanY->UseCurrentStyle();
364 fGraphMeanY->GetXaxis()->SetTitle("Event [#]");
365 fGraphMeanY->GetYaxis()->SetTitle("Y [\\circ]");
366
367 pad->Divide(2,2);
368
369 pad->cd(1);
370 fHistMeanX.DrawClone();
371 pad->cd(2);
372 fHistMeanY.DrawClone();
373 pad->cd(3);
374 fGraphMeanX->SetMarkerStyle(20);
375 fGraphMeanX->SetMarkerSize(0.4);
376 fGraphMeanX->Draw("AP");
377 pad->cd(4);
378 fGraphMeanY->SetMarkerStyle(20);
379 fGraphMeanY->SetMarkerSize(0.4);
380 fGraphMeanY->Draw("AP");
381 }
382 else if (option.Contains("sig"))
383 {
384
385 fGraphSigmaMinor = new TGraph((Int_t)fNumEvents,fEvent.GetArray(),fvsTimeSigmaMinor.GetArray());
386 fGraphSigmaMinor->SetName("GraphSigmaMinor");
387 fGraphSigmaMinor->SetTitle("Star X sigma v.s. event ");
388 fGraphSigmaMinor->UseCurrentStyle();
389 fGraphSigmaMinor->GetXaxis()->SetTitle("Event [#]");
390 fGraphSigmaMinor->GetYaxis()->SetTitle("Sigma X [\\circ]");
391
392 fGraphSigmaMajor = new TGraph((Int_t)fNumEvents,fEvent.GetArray(),fvsTimeSigmaMajor.GetArray());
393 fGraphSigmaMajor->SetName("GraphSigmaMajor");
394 fGraphSigmaMajor->SetTitle("Star Y sigma v.s. event ");
395 fGraphSigmaMajor->UseCurrentStyle();
396 fGraphSigmaMajor->GetXaxis()->SetTitle("Event [#]");
397 fGraphSigmaMajor->GetYaxis()->SetTitle("Sigma Y [\\circ]");
398
399 pad->Divide(2,2);
400
401 pad->cd(1);
402 fHistSigmaMinor.DrawClone();
403 pad->cd(2);
404 fHistSigmaMajor.DrawClone();
405 pad->cd(3);
406 fGraphSigmaMinor->SetMarkerStyle(21);
407 fGraphSigmaMinor->SetMarkerSize(0.4);
408// fGraphSigmaMinor->SetMarkerColor(kBlue);
409 fGraphSigmaMinor->Draw("AP");
410 fGraphSigmaMajor->SetMarkerStyle(21);
411 fGraphSigmaMajor->SetMarkerSize(0.4);
412 pad->cd(4);
413 fGraphSigmaMajor->Draw("AP");
414 }
415 else if (option.Contains("proj"))
416 {
417 pad->Divide(2,1);
418
419 pad->cd(1);
420// Creates a Root function based on function fitf above
421 TF1 *funcX = new TF1("fitfuncX",fitfunc,-0.4,0.4,4);
422
423 Float_t sigguess = fHistSigmaMinor.GetMean();
424 Float_t max = fProjectionX.GetMaximum();
425// Sets initial values and parameter names
426 funcX->SetParNames("base","Max","Mean","Sigma");
427 funcX->SetParameters(1.,max,0.,sigguess);
428
429// Fit histogram in range defined by function
430 fProjectionX.Fit("fitfuncX","QR");
431 // fProjectionX.DrawClone();
432
433 pad->cd(2);
434 TF1 *funcY = new TF1("fitfuncY",fitfunc,-0.4,0.4,4);
435
436 sigguess = fHistSigmaMajor.GetMean();
437 max = fProjectionY.GetMaximum();
438// Sets initial values and parameter names
439 funcY->SetParNames("base","Max","Mean","Sigma");
440 funcY->SetParameters(1.,max,0.,sigguess);
441
442// Fit histogram in range defined by function
443 fProjectionY.Fit("fitfuncX","QR");
444 fProjectionY.DrawClone();
445 }
446// else if (option.Contains("prof"))
447// {
448// fProfile.DrawClone();
449// }
450 else if (option.Contains("cam"))
451 {
452 pad->cd(1);
453
454 TArrayF x(fNumEvents);
455 TArrayF y(fNumEvents);
456 for (UInt_t i=0; i<fNumEvents; i++)
457 {
458 x[i] = fvsTimeMeanX[i]/fGeom->GetConvMm2Deg();
459 y[i] = fvsTimeMeanY[i]/fGeom->GetConvMm2Deg();
460 }
461
462 fGraphPath = new TGraph((Int_t)fNumEvents,x.GetArray(),y.GetArray());
463 fGraphPath->SetName("GraphPath");
464 fGraphPath->UseCurrentStyle();
465
466 fCamera.DrawClone();
467 fCamera.SetPrettyPalette();
468 fGraphPath->Draw("P");
469 fGraphPath->SetMarkerStyle(21);
470 fGraphPath->SetMarkerSize(0.4);
471 }
472 else
473 *fLog << err << GetName() << "::Draw Uknown option " << option << endl;
474
475}
476
Note: See TracBrowser for help on using the repository browser.