| 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): Roger Firpo 04/2004 <mailto:jlopez@ifae.es>
|
|---|
| 19 | ! Author(s): Javier López 05/2004 <mailto:jlopez@ifae.es>
|
|---|
| 20 | !
|
|---|
| 21 | ! Copyright: MAGIC Software Development, 2000-2004
|
|---|
| 22 | !
|
|---|
| 23 | !
|
|---|
| 24 | \* ======================================================================== */
|
|---|
| 25 |
|
|---|
| 26 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 27 | //
|
|---|
| 28 | // MSrcPosFromStars
|
|---|
| 29 | //
|
|---|
| 30 | // Task to set the position of the source using the positions of the stars
|
|---|
| 31 | // in the field of view of the source
|
|---|
| 32 | //
|
|---|
| 33 | // Output Containers:
|
|---|
| 34 | // MSrcPosCam
|
|---|
| 35 | //
|
|---|
| 36 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 37 |
|
|---|
| 38 |
|
|---|
| 39 | #include "MSrcPosFromStars.h"
|
|---|
| 40 |
|
|---|
| 41 | #include <TList.h>
|
|---|
| 42 | #include <TH2F.h>
|
|---|
| 43 | #include <TF2.h>
|
|---|
| 44 | #include <TTimer.h>
|
|---|
| 45 | #include <TString.h>
|
|---|
| 46 | #include <TCanvas.h>
|
|---|
| 47 |
|
|---|
| 48 | #include "MSrcPosCam.h"
|
|---|
| 49 | #include "MStarLocalCam.h"
|
|---|
| 50 | #include "MStarLocalPos.h"
|
|---|
| 51 |
|
|---|
| 52 | #include "MLog.h"
|
|---|
| 53 | #include "MLogManip.h"
|
|---|
| 54 |
|
|---|
| 55 | #include "MParList.h"
|
|---|
| 56 |
|
|---|
| 57 | ClassImp(MSrcPosFromStars);
|
|---|
| 58 |
|
|---|
| 59 | using namespace std;
|
|---|
| 60 |
|
|---|
| 61 | static const TString gsDefName = "MSrcPosFromStars";
|
|---|
| 62 | static const TString gsDefTitle = "task to calculate the position of the source using the position of stars";
|
|---|
| 63 |
|
|---|
| 64 | // -------------------------------------------------------------------------
|
|---|
| 65 | //
|
|---|
| 66 | // Default constructor.
|
|---|
| 67 | //
|
|---|
| 68 | MSrcPosFromStars::MSrcPosFromStars(const char *name, const char *title)
|
|---|
| 69 | : fStars(NULL)
|
|---|
| 70 | {
|
|---|
| 71 | fName = name ? name : gsDefName.Data();
|
|---|
| 72 | fTitle = title ? title : gsDefTitle.Data();
|
|---|
| 73 |
|
|---|
| 74 | fDistances.Set(0);
|
|---|
| 75 |
|
|---|
| 76 | }
|
|---|
| 77 |
|
|---|
| 78 | // -------------------------------------------------------------------------
|
|---|
| 79 | //
|
|---|
| 80 | Int_t MSrcPosFromStars::PreProcess(MParList *pList)
|
|---|
| 81 | {
|
|---|
| 82 |
|
|---|
| 83 | if(!MSrcPlace::PreProcess(pList))
|
|---|
| 84 | return kFALSE;
|
|---|
| 85 |
|
|---|
| 86 | fStars = (MStarLocalCam*)pList->FindObject(AddSerialNumber("MStarLocalCam"));
|
|---|
| 87 | if (!fStars)
|
|---|
| 88 | {
|
|---|
| 89 | *fLog << err << AddSerialNumber("MStarLocalCam") << " not found ... aborting" << endl;
|
|---|
| 90 | return kFALSE;
|
|---|
| 91 | }
|
|---|
| 92 |
|
|---|
| 93 | fNumStars = fDistances.GetSize();
|
|---|
| 94 |
|
|---|
| 95 | return kTRUE;
|
|---|
| 96 | }
|
|---|
| 97 |
|
|---|
| 98 | // --------------------------------------------------------------------------
|
|---|
| 99 | //
|
|---|
| 100 | //
|
|---|
| 101 | static Bool_t HandleInput()
|
|---|
| 102 | {
|
|---|
| 103 | TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
|
|---|
| 104 | while (1)
|
|---|
| 105 | {
|
|---|
| 106 | //
|
|---|
| 107 | // While reading the input process gui events asynchronously
|
|---|
| 108 | //
|
|---|
| 109 | timer.TurnOn();
|
|---|
| 110 | cout << "Type 'q' to exit, <return> to go on: " << endl;
|
|---|
| 111 | TString input;
|
|---|
| 112 | input = getchar();
|
|---|
| 113 | timer.TurnOff();
|
|---|
| 114 |
|
|---|
| 115 | if (input=="q\n")
|
|---|
| 116 | return kFALSE;
|
|---|
| 117 |
|
|---|
| 118 | if (input=="\n")
|
|---|
| 119 | return kTRUE;
|
|---|
| 120 | };
|
|---|
| 121 |
|
|---|
| 122 | return kFALSE;
|
|---|
| 123 | }
|
|---|
| 124 |
|
|---|
| 125 | Int_t MSrcPosFromStars::ComputeNewSrcPosition()
|
|---|
| 126 | {
|
|---|
| 127 |
|
|---|
| 128 | if (fNumStars == 0)
|
|---|
| 129 | {
|
|---|
| 130 | if (fStars->GetNumStars() > 0)
|
|---|
| 131 | {
|
|---|
| 132 |
|
|---|
| 133 | //Look for the star closer to the center of the camera
|
|---|
| 134 | TIter Next(fStars->GetList());
|
|---|
| 135 | MStarLocalPos* star;
|
|---|
| 136 | Float_t mindist = 600; //mm
|
|---|
| 137 | UInt_t starnum = 0;
|
|---|
| 138 | Int_t select = -1;
|
|---|
| 139 | while ((star=(MStarLocalPos*)Next()))
|
|---|
| 140 | {
|
|---|
| 141 | Float_t dist = TMath::Sqrt(star->GetMeanX()*star->GetMeanX() +
|
|---|
| 142 | star->GetMeanY()*star->GetMeanY());
|
|---|
| 143 | if (dist < mindist)
|
|---|
| 144 | {
|
|---|
| 145 | mindist = dist;
|
|---|
| 146 | select = starnum;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | starnum++;
|
|---|
| 150 | }
|
|---|
| 151 |
|
|---|
| 152 | if (select < 0)
|
|---|
| 153 | {
|
|---|
| 154 | *fLog << err << "Not found star closer to center" << endl;
|
|---|
| 155 | return kFALSE;
|
|---|
| 156 | }
|
|---|
| 157 |
|
|---|
| 158 | MStarLocalPos& selecStar = (*fStars)[select];
|
|---|
| 159 |
|
|---|
| 160 | if (selecStar.GetChiSquareNdof() > 0. && selecStar.GetChiSquareNdof() < 10.)
|
|---|
| 161 | {
|
|---|
| 162 |
|
|---|
| 163 | Float_t selecStarPosX = selecStar.GetMeanX();
|
|---|
| 164 | Float_t selecStarPosY = selecStar.GetMeanY();
|
|---|
| 165 |
|
|---|
| 166 | GetOutputSrcPosCam()->SetXY(selecStarPosX,selecStarPosY);
|
|---|
| 167 | }
|
|---|
| 168 | }
|
|---|
| 169 | }
|
|---|
| 170 | else if (fStars->GetNumStars() >= fNumStars)
|
|---|
| 171 | {
|
|---|
| 172 |
|
|---|
| 173 | Float_t diameterInnerPixel = 30; //[mm]
|
|---|
| 174 | Float_t diameterOuterPixel = 60; //[mm]
|
|---|
| 175 |
|
|---|
| 176 | Double_t probability = 1.;
|
|---|
| 177 |
|
|---|
| 178 | // Resolution and computing time are proportional to bins^2
|
|---|
| 179 | const Int_t bins = 200;
|
|---|
| 180 | Double_t max_x_mm = 600;
|
|---|
| 181 | Double_t min_x_mm = -max_x_mm;
|
|---|
| 182 | Double_t max_y_mm = max_x_mm;
|
|---|
| 183 | Double_t min_y_mm = -max_x_mm;
|
|---|
| 184 |
|
|---|
| 185 | // bins to mmrees
|
|---|
| 186 | const Double_t bin2mm = 2 * max_x_mm / bins;
|
|---|
| 187 |
|
|---|
| 188 | // First run, to find the maximum peak and define the area
|
|---|
| 189 | TH2F *hgrid = new TH2F("hgrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm);
|
|---|
| 190 |
|
|---|
| 191 | for (Int_t ny = 1; ny <= bins; ny++)
|
|---|
| 192 | for (Int_t nx = 1; nx <= bins; nx++)
|
|---|
| 193 | hgrid->SetBinContent(nx, ny, 1.);
|
|---|
| 194 |
|
|---|
| 195 | for (UInt_t numstar = 0; numstar < fNumStars; numstar++)
|
|---|
| 196 | {
|
|---|
| 197 | probability = 1;
|
|---|
| 198 |
|
|---|
| 199 | MStarLocalPos& star = (*fStars)[numstar];
|
|---|
| 200 |
|
|---|
| 201 | if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.)
|
|---|
| 202 | {
|
|---|
| 203 |
|
|---|
| 204 | Float_t starPosX = star.GetMeanX();
|
|---|
| 205 | Float_t starPosY = star.GetMeanY();
|
|---|
| 206 | Float_t starDist = Dist(0.,0.,starPosX,starPosY);
|
|---|
| 207 | Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel);
|
|---|
| 208 |
|
|---|
| 209 | // *fLog << dbg << "Star[" << numstar << "] pos (" << starPosX << "," << starPosY << ") dist " << starDist << " size " << starSigma << endl;
|
|---|
| 210 |
|
|---|
| 211 | if (starSigma != 0)
|
|---|
| 212 | {
|
|---|
| 213 | for (Int_t ny = 1; ny < bins + 1; ny++)
|
|---|
| 214 | {
|
|---|
| 215 | for (Int_t nx = 0; nx < bins + 1; nx++)
|
|---|
| 216 | {
|
|---|
| 217 | Float_t dist = Dist(min_x_mm + nx * bin2mm, starPosX, min_y_mm + ny * bin2mm, starPosY);
|
|---|
| 218 | Float_t prob = Prob(dist, fDistances[numstar], starSigma);
|
|---|
| 219 |
|
|---|
| 220 | // if ( prob > 0.8)
|
|---|
| 221 | // *fLog << dbg << " x " << min_x_mm + nx * bin2mm << " y " << min_y_mm + ny * bin2mm << " dist " << dist << " stardist " << fDistances[numstar] << " prob " << prob << endl;
|
|---|
| 222 |
|
|---|
| 223 | probability = hgrid->GetBinContent(nx, ny)*prob;
|
|---|
| 224 | hgrid->SetBinContent(nx, ny, probability);
|
|---|
| 225 |
|
|---|
| 226 | }
|
|---|
| 227 | }
|
|---|
| 228 | }
|
|---|
| 229 | else probability *= 1;
|
|---|
| 230 |
|
|---|
| 231 | }
|
|---|
| 232 | }
|
|---|
| 233 |
|
|---|
| 234 | // Finding the peak
|
|---|
| 235 | Double_t peak[2] = {0,0};
|
|---|
| 236 | Double_t peak_value = 0;
|
|---|
| 237 |
|
|---|
| 238 | for (Int_t ny = 0; ny < bins + 1; ny++)
|
|---|
| 239 | {
|
|---|
| 240 | for (Int_t nx = 0; nx < bins + 1; nx++)
|
|---|
| 241 | {
|
|---|
| 242 | if (hgrid->GetBinContent(nx, ny) > peak_value)
|
|---|
| 243 | {
|
|---|
| 244 | peak_value = hgrid->GetBinContent(nx, ny);
|
|---|
| 245 | peak[0] = min_x_mm + nx * bin2mm;
|
|---|
| 246 | peak[1] = min_y_mm + ny * bin2mm;
|
|---|
| 247 | }
|
|---|
| 248 | }
|
|---|
| 249 | }
|
|---|
| 250 |
|
|---|
| 251 | *fLog << dbg << "The peak is at (x, y) = (" << peak[0] << ", " << peak[1] << ") mm" << endl;
|
|---|
| 252 |
|
|---|
| 253 |
|
|---|
| 254 | // TCanvas c1;
|
|---|
| 255 | // hgrid->Draw("lego");
|
|---|
| 256 | // if(!HandleInput())
|
|---|
| 257 | // exit(1);
|
|---|
| 258 |
|
|---|
| 259 |
|
|---|
| 260 | // Here we define a small area surrounding the peak to avoid wasting time and resolution anywhere else
|
|---|
| 261 |
|
|---|
| 262 | min_x_mm = peak[0] - diameterInnerPixel;
|
|---|
| 263 | max_x_mm = peak[0] + diameterInnerPixel;
|
|---|
| 264 | min_y_mm = peak[1] - diameterInnerPixel;
|
|---|
| 265 | max_y_mm = peak[1] + diameterInnerPixel;
|
|---|
| 266 |
|
|---|
| 267 | const Double_t xbin2mm = (max_x_mm - min_x_mm) / bins;
|
|---|
| 268 | const Double_t ybin2mm = (max_y_mm - min_y_mm) / bins;
|
|---|
| 269 |
|
|---|
| 270 | // Other run centered in the peak for more precision
|
|---|
| 271 | TH2F *hagrid = new TH2F("hagrid", "", bins, min_x_mm, max_x_mm, bins, min_y_mm, max_y_mm);
|
|---|
| 272 |
|
|---|
| 273 | for (Int_t ny = 1; ny <= bins; ny++)
|
|---|
| 274 | for (Int_t nx = 1; nx <= bins; nx++)
|
|---|
| 275 | hagrid->SetBinContent(nx, ny, 1.);
|
|---|
| 276 |
|
|---|
| 277 |
|
|---|
| 278 | for (UInt_t numstar = 0; numstar < fNumStars; numstar++)
|
|---|
| 279 | {
|
|---|
| 280 | probability = 1;
|
|---|
| 281 |
|
|---|
| 282 | MStarLocalPos& star = (*fStars)[numstar];
|
|---|
| 283 |
|
|---|
| 284 | if (star.GetChiSquareNdof() > 0. && star.GetChiSquareNdof() < 10.)
|
|---|
| 285 | {
|
|---|
| 286 |
|
|---|
| 287 | Float_t starPosX = star.GetMeanX();
|
|---|
| 288 | Float_t starPosY = star.GetMeanY();
|
|---|
| 289 | Float_t starDist = Dist(0.,0.,starPosX,starPosY);
|
|---|
| 290 | Float_t starSigma = (starDist<350.?diameterInnerPixel:diameterOuterPixel);
|
|---|
| 291 |
|
|---|
| 292 | if (starSigma != 0)
|
|---|
| 293 | {
|
|---|
| 294 | for (Int_t ny = 0; ny < bins; ny++)
|
|---|
| 295 | {
|
|---|
| 296 | for (Int_t nx = 0; nx < bins; nx++)
|
|---|
| 297 | {
|
|---|
| 298 | Float_t prob = Prob(Dist(min_x_mm + nx * xbin2mm, starPosX, min_y_mm + ny * ybin2mm, starPosY), fDistances[numstar], starSigma);
|
|---|
| 299 |
|
|---|
| 300 | probability = hagrid->GetBinContent(nx, ny)*prob;
|
|---|
| 301 | hagrid->SetBinContent(nx, ny, probability);
|
|---|
| 302 | }
|
|---|
| 303 | }
|
|---|
| 304 | }
|
|---|
| 305 | else probability *= 1;
|
|---|
| 306 |
|
|---|
| 307 | }
|
|---|
| 308 | }
|
|---|
| 309 |
|
|---|
| 310 | // This time we fit the histogram with a 2D gaussian
|
|---|
| 311 | // Although we don't expect our function to be gaussian...
|
|---|
| 312 | TF2 *gauss2d = new TF2("gauss2d","[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))*exp(-(y-[3])*(y-[3])/(2*[4]*[4]))", min_x_mm, max_x_mm, min_y_mm, max_y_mm);
|
|---|
| 313 |
|
|---|
| 314 | gauss2d->SetParName(0,"f0");
|
|---|
| 315 | gauss2d->SetParName(1,"meanx");
|
|---|
| 316 | gauss2d->SetParName(2,"sigmax");
|
|---|
| 317 | gauss2d->SetParName(3,"meany");
|
|---|
| 318 | gauss2d->SetParName(4,"sigmay");
|
|---|
| 319 |
|
|---|
| 320 | gauss2d->SetParameter(0,10);
|
|---|
| 321 | gauss2d->SetParameter(1,peak[0]);
|
|---|
| 322 | gauss2d->SetParameter(2,0.002);
|
|---|
| 323 | gauss2d->SetParameter(3,peak[1]);
|
|---|
| 324 | gauss2d->SetParameter(4,0.002);
|
|---|
| 325 |
|
|---|
| 326 | GetOutputSrcPosCam()->SetXY(gauss2d->GetParameter(1), gauss2d->GetParameter(3));
|
|---|
| 327 |
|
|---|
| 328 | delete hgrid;
|
|---|
| 329 | delete hagrid;
|
|---|
| 330 | }
|
|---|
| 331 |
|
|---|
| 332 |
|
|---|
| 333 | return kTRUE;
|
|---|
| 334 | }
|
|---|
| 335 |
|
|---|
