Changeset 4410 for trunk/MagicSoft
- Timestamp:
- 07/20/04 12:15:06 (20 years ago)
- Location:
- trunk/MagicSoft/Mars/mtemp/mifae/library
- Files:
-
- 5 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Mars/mtemp/mifae/library/MLiveTime.cc
r4405 r4410 69 69 totalLiveTime += fLiveTimeBin[bin]; 70 70 } 71 *fLog << GetName() << ": Total live time " << totalLiveTime << endl;71 *fLog << GetName() << ": Total live time " << totalLiveTime << " sec" << endl; 72 72 } 73 73 else -
trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPlace.h
r4117 r4410 25 25 26 26 TH2F* fHistPos; // histogram of the used source positions 27 OnOffMode_t fMode; // On/Off data mode (write/read to/from the histogram)27 OnOffMode_t fMode; // On/Off data mode (write/read to/from the histogram) 28 28 Bool_t fCreateHisto; // flag to decide whether internal histo is created or not 29 29 … … 48 48 void SetInternalHistoName(TString name) {fHistoName=name;} 49 49 void SetInternalHistoBinSize(Float_t size){fHistoBinPrec=size;} 50 void SetPositionHisto(TH2F* histo) {fHistPos=histo;} 50 51 void SetCreateHisto(Bool_t inp=kTRUE) {fCreateHisto=inp;} 51 52 -
trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.cc
r4295 r4410 40 40 41 41 #include <TList.h> 42 #include <TH2F.h> 43 #include <TF2.h> 44 #include <TTimer.h> 45 #include <TString.h> 46 #include <TCanvas.h> 42 47 43 48 #include "MSrcPosCam.h" … … 61 66 // Default constructor. 62 67 // 63 MSrcPosFromStars::MSrcPosFromStars( Float_t first, Float_t second,const char *name, const char *title)64 : fStars(NULL) , fNumStars(2), fDistanceFirstStar(first), fDistanceSecondStar(second)68 MSrcPosFromStars::MSrcPosFromStars(const char *name, const char *title) 69 : fStars(NULL) 65 70 { 66 71 fName = name ? name : gsDefName.Data(); 67 72 fTitle = title ? title : gsDefTitle.Data(); 73 74 fDistances.Set(0); 68 75 69 76 } … … 84 91 } 85 92 93 fNumStars = fDistances.GetSize(); 94 86 95 return kTRUE; 87 96 } … … 90 99 // 91 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 92 125 Int_t MSrcPosFromStars::ComputeNewSrcPosition() 93 126 { 94 127 95 if (f DistanceFirstStar == 0. || fDistanceSecondStar == 0.)128 if (fNumStars == 0) 96 129 { 97 130 if (fStars->GetNumStars() > 0) … … 135 168 } 136 169 } 137 else if (fStars->GetNumStars() == fNumStars)170 else if (fStars->GetNumStars() >= fNumStars) 138 171 { 139 172 140 MStarLocalPos& firstStar = (*fStars)[0]; 141 MStarLocalPos& secondStar = (*fStars)[1]; 142 143 if (firstStar.GetChiSquareNdof() > 0. && firstStar.GetChiSquareNdof() < 10. && 144 secondStar.GetChiSquareNdof() > 0. && secondStar.GetChiSquareNdof() < 10.) 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++) 145 196 { 146 147 Float_t firstStarPosX = firstStar.GetMeanX(); 148 Float_t firstStarPosY = firstStar.GetMeanY(); 149 150 Float_t secondStarPosX = secondStar.GetMeanX(); 151 Float_t secondStarPosY = secondStar.GetMeanY(); 152 153 Float_t distanceStars = sqrt(pow(firstStarPosX - secondStarPosX,2) + pow(firstStarPosY - secondStarPosY,2)); 154 Float_t sin_alpha = (secondStarPosY - firstStarPosY) / distanceStars; 155 Float_t cos_alpha = (secondStarPosX - firstStarPosX) / distanceStars; 156 157 Float_t x = (pow(fDistanceFirstStar,2) - pow(fDistanceSecondStar,2) + pow(distanceStars,2)) / (2 * distanceStars); 158 159 Float_t arg = 4 * pow(distanceStars,2) * pow(fDistanceFirstStar,2) - pow(pow(fDistanceFirstStar,2) - pow(fDistanceSecondStar,2) + pow(distanceStars,2),2); 160 161 if (arg >= 0.) 162 { 163 164 Float_t y = sqrt(arg) / (2 * distanceStars); 165 166 Float_t xc1 = firstStarPosX + x * cos_alpha - y * sin_alpha; 167 Float_t yc1 = firstStarPosY + x * sin_alpha + y * cos_alpha; 168 169 Float_t xc2 = firstStarPosX + x * cos_alpha + y * sin_alpha; 170 Float_t yc2 = firstStarPosY + x * sin_alpha - y * cos_alpha; 171 172 if (sqrt(xc1*xc1 + yc1*yc1) < sqrt(xc2*xc2 + yc2*yc2)) 173 GetOutputSrcPosCam()->SetXY(xc1,yc1); 174 else 175 GetOutputSrcPosCam()->SetXY(xc2,yc2); 176 } 177 else 178 *fLog << warn << GetName() << " negative argument [" << arg << "] for sqrt()" << endl; 179 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 } 180 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; 181 330 } 182 331 332 183 333 return kTRUE; 184 334 } -
trunk/MagicSoft/Mars/mtemp/mifae/library/MSrcPosFromStars.h
r4295 r4410 1 1 #ifndef MARS_MSrcPosFromStars 2 2 #define MARS_MSrcPosFromStars 3 4 #ifndef ROOT_TArrayF 5 #include "TArrayF.h" 6 #endif 7 8 #ifndef ROOT_TMath 9 #include "TMath.h" 10 #endif 3 11 4 12 #ifndef MARS_MSrcPlace … … 15 23 MStarLocalCam *fStars; 16 24 17 const UInt_t fNumStars; 18 Float_t fDistanceFirstStar; 19 Float_t fDistanceSecondStar; 25 UInt_t fNumStars; 26 TArrayF fDistances; 20 27 21 28 virtual Int_t ComputeNewSrcPosition(); … … 24 31 public: 25 32 26 MSrcPosFromStars(Float_t first=0., Float_t second=0., const char *name=NULL, const char *title=NULL); 33 MSrcPosFromStars(const char *name=NULL, const char *title=NULL); 34 35 36 // This fuction is to compute the probability using a ring distribution 37 Double_t Prob(Double_t d, Double_t dist, Double_t sigma) 38 { 39 return TMath::Exp(-(d-dist)*(d-dist)/(2*sigma*sigma));///(sigma*TMath::Sqrt(2*3.141592654)); 40 } 41 42 Double_t Dist(Double_t x1, Double_t x2, Double_t y1, Double_t y2) 43 { 44 return TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)); 45 } 27 46 28 void SetDistanceFirstStar(Float_t dist) { fDistanceFirstStar = dist; } 29 void SetDistanceSecondStar(Float_t dist) { fDistanceSecondStar = dist; } 47 void SetDistances(TArrayF dist) { fDistances = dist; } 30 48 31 49 ClassDef(MSrcPosFromStars, 0) // task to calculate the position of the source using the position of stars -
trunk/MagicSoft/Mars/mtemp/mifae/library/Makefile
r4409 r4410 66 66 MDispCalc.cc \ 67 67 MControlPlots.cc \ 68 MSrcPosFromStars.c \68 MSrcPosFromStars.cc \ 69 69 MLiveTime.cc \ 70 70 MLiveTimeCalc.cc
Note:
See TracChangeset
for help on using the changeset viewer.