| 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): Thomas Bretz, 2003 <mailto:tbretz@astro.uni-wuerzburg.de>
|
|---|
| 19 | !
|
|---|
| 20 | ! Copyright: MAGIC Software Development, 2000-2003
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 | /////////////////////////////////////////////////////////////////////////////
|
|---|
| 26 | //
|
|---|
| 27 | // MBending
|
|---|
| 28 | //
|
|---|
| 29 | // Double_t fIe ; // [rad] Index Error in Elevation
|
|---|
| 30 | // Double_t fIa ; // [rad] Index Error in Azimuth
|
|---|
| 31 | //
|
|---|
| 32 | // Double_t fFlop ; // [rad] Vertical Sag
|
|---|
| 33 | // * do not use if not data: Zd<0
|
|---|
| 34 | //
|
|---|
| 35 | // Double_t fNpae ; // [rad] Az-El Nonperpendicularity
|
|---|
| 36 | //
|
|---|
| 37 | // Double_t fCa ; // [rad] Left-Right Collimation Error
|
|---|
| 38 | //
|
|---|
| 39 | // Double_t fAn ; // [rad] Azimuth Axis Misalignment (N-S)
|
|---|
| 40 | // Double_t fAw ; // [rad] Azimuth Axis Misalignment (E-W)
|
|---|
| 41 | //
|
|---|
| 42 | // Double_t fTf ; // [rad] Tube fluxture (sin)
|
|---|
| 43 | // * same as ecec if no data: Zd<0
|
|---|
| 44 | // Double_t fTx ; // [rad] Tube fluxture (tan)
|
|---|
| 45 | // * do not use with NPAE if no data: Zd<0
|
|---|
| 46 | //
|
|---|
| 47 | // Double_t fNrx ; // [rad] Nasmyth rotator displacement, horizontan
|
|---|
| 48 | // Double_t fNry ; // [rad] Nasmyth rotator displacement, vertical
|
|---|
| 49 | //
|
|---|
| 50 | // Double_t fCrx ; // [rad] Alt/Az Coude Displacement (N-S)
|
|---|
| 51 | // Double_t fCry ; // [rad] Alt/Az Coude Displacement (E-W)
|
|---|
| 52 | //
|
|---|
| 53 | // Double_t fEces ; // [rad] Elevation Centering Error (sin)
|
|---|
| 54 | // Double_t fAces ; // [rad] Azimuth Centering Error (sin)
|
|---|
| 55 | // Double_t fEcec ; // [rad] Elevation Centering Error (cos)
|
|---|
| 56 | // Double_t fAcec ; // [rad] Azimuth Centering Error (cos)
|
|---|
| 57 | //
|
|---|
| 58 | ////////////////////////////////////////////////////////////////////////////
|
|---|
| 59 | #include "MBending.h"
|
|---|
| 60 |
|
|---|
| 61 | #include <fstream.h>
|
|---|
| 62 | #include <iomanip.h>
|
|---|
| 63 | #include <TVector3.h>
|
|---|
| 64 |
|
|---|
| 65 | #include <TMinuit.h>
|
|---|
| 66 |
|
|---|
| 67 | #include "MTime.h"
|
|---|
| 68 |
|
|---|
| 69 | ClassImp(MBending);
|
|---|
| 70 |
|
|---|
| 71 | #undef DEBUG
|
|---|
| 72 | //#define DEBUG(txt) txt
|
|---|
| 73 | #define DEBUG(txt)
|
|---|
| 74 |
|
|---|
| 75 | const Int_t MBending::fNumPar=19;
|
|---|
| 76 |
|
|---|
| 77 | void MBending::Init()
|
|---|
| 78 | {
|
|---|
| 79 | fCoeff = new Double_t*[fNumPar];
|
|---|
| 80 | fName = new TString[fNumPar];
|
|---|
| 81 | fDescr = new TString[fNumPar];
|
|---|
| 82 |
|
|---|
| 83 | fCoeff[ 0] = &fIa; fName[ 0] = "IA";
|
|---|
| 84 | fCoeff[ 1] = &fIe; fName[ 1] = "IE";
|
|---|
| 85 | fCoeff[ 2] = &fFlop; fName[ 2] = "FLOP";
|
|---|
| 86 | fCoeff[ 3] = &fAn; fName[ 3] = "AN";
|
|---|
| 87 | fCoeff[ 4] = &fAw; fName[ 4] = "AW";
|
|---|
| 88 | fCoeff[ 5] = &fNpae; fName[ 5] = "NPAE";
|
|---|
| 89 | fCoeff[ 6] = &fCa; fName[ 6] = "CA";
|
|---|
| 90 | fCoeff[ 7] = &fTf; fName[ 7] = "TF";
|
|---|
| 91 | fCoeff[ 8] = &fTx; fName[ 8] = "TX";
|
|---|
| 92 | fCoeff[ 9] = &fEces; fName[ 9] = "ECES";
|
|---|
| 93 | fCoeff[10] = &fAces; fName[10] = "ACES";
|
|---|
| 94 | fCoeff[11] = &fEcec; fName[11] = "ECEC";
|
|---|
| 95 | fCoeff[12] = &fAcec; fName[12] = "ACEC";
|
|---|
| 96 | fCoeff[13] = &fNrx; fName[13] = "NRX";
|
|---|
| 97 | fCoeff[14] = &fNry; fName[14] = "NRY";
|
|---|
| 98 | fCoeff[15] = &fCrx; fName[15] = "CRX";
|
|---|
| 99 | fCoeff[16] = &fCry; fName[16] = "CRY";
|
|---|
| 100 | fCoeff[17] = &fMagic1; fName[17] = "MAGIC1";
|
|---|
| 101 | fCoeff[18] = &fMagic2; fName[18] = "MAGIC2";
|
|---|
| 102 |
|
|---|
| 103 | fDescr[ 0] = "Index Error Azimuth";
|
|---|
| 104 | fDescr[ 1] = "Index Error Zenith Distance";
|
|---|
| 105 | fDescr[ 2] = "Vertical Sag";
|
|---|
| 106 | fDescr[ 3] = "Azimuth Axis Misalignment (N-S)";
|
|---|
| 107 | fDescr[ 4] = "Azimuth Axis Misalignment (E-W)";
|
|---|
| 108 | fDescr[ 5] = "Az-El Nonperpendicularity";
|
|---|
| 109 | fDescr[ 6] = "Left-Right Collimation Error";
|
|---|
| 110 | fDescr[ 7] = "Tube fluxture (sin)";
|
|---|
| 111 | fDescr[ 8] = "Tube fluxture (tan)";
|
|---|
| 112 | fDescr[ 9] = "Elevation Centering Error (sin)";
|
|---|
| 113 | fDescr[10] = "Azimuth Centering Error (sin)";
|
|---|
| 114 | fDescr[11] = "Elevation Centering Error (cos)";
|
|---|
| 115 | fDescr[12] = "Azimuth Centering Error (cos)";
|
|---|
| 116 | fDescr[13] = "Nasmyth rotator displacement (horizontal)";
|
|---|
| 117 | fDescr[14] = "Nasmyth rotator displacement (vertical)";
|
|---|
| 118 | fDescr[15] = "Alt/Az Coude Displacement (N-S)";
|
|---|
| 119 | fDescr[16] = "Alt/Az Coude Displacement (E-W)";
|
|---|
| 120 | fDescr[17] = "n/a <ZA Hysteresis>";
|
|---|
| 121 | fDescr[18] = "n/a";
|
|---|
| 122 | }
|
|---|
| 123 | void MBending::Reset()
|
|---|
| 124 | {
|
|---|
| 125 | Clear();
|
|---|
| 126 | }
|
|---|
| 127 |
|
|---|
| 128 | void MBending::Load(const char *name)
|
|---|
| 129 | {
|
|---|
| 130 | /*
|
|---|
| 131 | ! MMT 1987 July 8
|
|---|
| 132 | ! T 36 7.3622 41.448 -0.0481
|
|---|
| 133 | ! IA -37.5465 20.80602
|
|---|
| 134 | ! IE -13.9180 1.25217
|
|---|
| 135 | ! NPAE +7.0751 26.44763
|
|---|
| 136 | ! CA -6.9149 32.05358
|
|---|
| 137 | ! AN +0.5053 1.40956
|
|---|
| 138 | ! AW -2.2016 1.37480
|
|---|
| 139 | ! END
|
|---|
| 140 | */
|
|---|
| 141 |
|
|---|
| 142 | ifstream fin(name);
|
|---|
| 143 | if (!fin)
|
|---|
| 144 | {
|
|---|
| 145 | cout << "Error: Cannot open file '" << name << "'" << endl;
|
|---|
| 146 | return;
|
|---|
| 147 | }
|
|---|
| 148 |
|
|---|
| 149 | char c;
|
|---|
| 150 | while (fin && fin.get()!='\n');
|
|---|
| 151 | fin >> c;
|
|---|
| 152 |
|
|---|
| 153 | if (c!='S' && c!='s')
|
|---|
| 154 | {
|
|---|
| 155 | cout << "Error: This in not a model correcting the star position (" << c << ")" << endl;
|
|---|
| 156 | return;
|
|---|
| 157 | }
|
|---|
| 158 |
|
|---|
| 159 | Clear();
|
|---|
| 160 |
|
|---|
| 161 | cout << endl;
|
|---|
| 162 |
|
|---|
| 163 | Double_t val;
|
|---|
| 164 | fin >> val;
|
|---|
| 165 | cout << "Number of observed stars: " << val << endl;
|
|---|
| 166 | fin >> val;
|
|---|
| 167 | cout << "Sky RMS: " << val << "\"" << endl;
|
|---|
| 168 | fin >> val;
|
|---|
| 169 | cout << "Refraction Constant A: " << val << "\"" << endl;
|
|---|
| 170 | fin >> val;
|
|---|
| 171 | cout << "Refraction Constant B: " << val << "\"" << endl;
|
|---|
| 172 |
|
|---|
| 173 | cout << endl;
|
|---|
| 174 |
|
|---|
| 175 | cout << " & = Name Value Sigma" << endl;
|
|---|
| 176 | cout << "--------------------------------------------------" << endl;
|
|---|
| 177 |
|
|---|
| 178 | while (fin)
|
|---|
| 179 | {
|
|---|
| 180 | TString str;
|
|---|
| 181 | fin >> str;
|
|---|
| 182 |
|
|---|
| 183 | if (str=="END")
|
|---|
| 184 | break;
|
|---|
| 185 |
|
|---|
| 186 | if (str[0]=='&')
|
|---|
| 187 | {
|
|---|
| 188 | cout << " & ";
|
|---|
| 189 | str.Remove(0);
|
|---|
| 190 | }
|
|---|
| 191 | else
|
|---|
| 192 | cout << " ";
|
|---|
| 193 |
|
|---|
| 194 | if (str[1]=='=')
|
|---|
| 195 | {
|
|---|
| 196 | cout << "= ";
|
|---|
| 197 | str.Remove(0);
|
|---|
| 198 | }
|
|---|
| 199 | else
|
|---|
| 200 | cout << " ";
|
|---|
| 201 |
|
|---|
| 202 | fin >> val;
|
|---|
| 203 | cout << str << "\t" << setw(11) << val << "° \t";
|
|---|
| 204 | val *= kDeg2Rad;
|
|---|
| 205 |
|
|---|
| 206 | // Find parameter
|
|---|
| 207 | Int_t n = -1;
|
|---|
| 208 | for (int i=0; i<fNumPar; i++)
|
|---|
| 209 | if (str==fName[i])
|
|---|
| 210 | {
|
|---|
| 211 | n = i;
|
|---|
| 212 | *fCoeff[i] = val;
|
|---|
| 213 | break;
|
|---|
| 214 | }
|
|---|
| 215 |
|
|---|
| 216 | fin >> val;
|
|---|
| 217 | cout << setw(9) << val << "°" << endl;
|
|---|
| 218 |
|
|---|
| 219 | // corresponding error
|
|---|
| 220 | fError[n] = val*kDeg2Rad;
|
|---|
| 221 | }
|
|---|
| 222 | cout << endl;
|
|---|
| 223 | }
|
|---|
| 224 |
|
|---|
| 225 | void MBending::Save(const char *name)
|
|---|
| 226 | {
|
|---|
| 227 | /*
|
|---|
| 228 | ! MMT 1987 July 8
|
|---|
| 229 | ! T 36 7.3622 41.448 -0.0481
|
|---|
| 230 | ! IA -37.5465 20.80602
|
|---|
| 231 | ! IE -13.9180 1.25217
|
|---|
| 232 | ! NPAE +7.0751 26.44763
|
|---|
| 233 | ! CA -6.9149 32.05358
|
|---|
| 234 | ! AN +0.5053 1.40956
|
|---|
| 235 | ! AW -2.2016 1.37480
|
|---|
| 236 | ! END
|
|---|
| 237 | */
|
|---|
| 238 |
|
|---|
| 239 | ofstream fout(name);
|
|---|
| 240 | if (!fout)
|
|---|
| 241 | {
|
|---|
| 242 | cout << "Error: Cannot open file '" << name << "'" << endl;
|
|---|
| 243 | return;
|
|---|
| 244 | }
|
|---|
| 245 |
|
|---|
| 246 | MTime t;
|
|---|
| 247 | t.Now();
|
|---|
| 248 |
|
|---|
| 249 | fout << "MAGIC1 " << t << endl;
|
|---|
| 250 | fout << "S 00 000000 000000 0000000" << endl;
|
|---|
| 251 | fout << setprecision(8);
|
|---|
| 252 | for (int i=0; i<fNumPar; i++)
|
|---|
| 253 | {
|
|---|
| 254 | fout << " " << setw(6) << GetName(i) << " ";
|
|---|
| 255 | fout << setw(13) << *fCoeff[i]*kRad2Deg << " ";
|
|---|
| 256 | fout << setw(11) << fError[i]*kRad2Deg << endl;
|
|---|
| 257 | }
|
|---|
| 258 | fout << "END" << endl;
|
|---|
| 259 | }
|
|---|
| 260 |
|
|---|
| 261 | Double_t MBending::Sign(Double_t val, Double_t alt)
|
|---|
| 262 | {
|
|---|
| 263 | // Some pointing corrections are defined as Delta ZA, which
|
|---|
| 264 | // is (P. Wallace) defined [0,90]deg while Alt is defined
|
|---|
| 265 | // [0,180]deg
|
|---|
| 266 | return (TMath::Pi()/2-alt < 0 ? -val : val);
|
|---|
| 267 | }
|
|---|
| 268 |
|
|---|
| 269 | AltAz MBending::AddOffsets(const AltAz &aa) const
|
|---|
| 270 | {
|
|---|
| 271 | // Correct [rad]
|
|---|
| 272 | // zdaz [rad]
|
|---|
| 273 | AltAz p = aa;
|
|---|
| 274 |
|
|---|
| 275 | const AltAz I(fIe, fIa);
|
|---|
| 276 | p += I;
|
|---|
| 277 |
|
|---|
| 278 | return p;
|
|---|
| 279 | }
|
|---|
| 280 |
|
|---|
| 281 | AltAz MBending::SubtractOffsets(const AltAz &aa) const
|
|---|
| 282 | {
|
|---|
| 283 | // Correct [rad]
|
|---|
| 284 | // zdaz [rad]
|
|---|
| 285 | AltAz p = aa;
|
|---|
| 286 |
|
|---|
| 287 | const AltAz I(fIe, fIa);
|
|---|
| 288 | p -= I;
|
|---|
| 289 |
|
|---|
| 290 | return p;
|
|---|
| 291 | }
|
|---|
| 292 |
|
|---|
| 293 | AltAz MBending::CalcAnAw(const AltAz &p, Int_t sign) const
|
|---|
| 294 | {
|
|---|
| 295 | // Corrections for AN and AW without approximations
|
|---|
| 296 | // as done by Patrick Wallace. The approximation cannot
|
|---|
| 297 | // be used for MAGIC because the correctioon angle
|
|---|
| 298 | // AW (~1.5deg) is not small enough.
|
|---|
| 299 |
|
|---|
| 300 | // Vector in cartesian coordinates
|
|---|
| 301 | TVector3 v1;
|
|---|
| 302 |
|
|---|
| 303 | // Set Azimuth and Elevation
|
|---|
| 304 | v1.SetMagThetaPhi(1, TMath::Pi()/2-p.Alt(), p.Az());
|
|---|
| 305 |
|
|---|
| 306 |
|
|---|
| 307 | TVector3 v2(v1);
|
|---|
| 308 | // cout << sign << endl;
|
|---|
| 309 |
|
|---|
| 310 | // cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
|
|---|
| 311 |
|
|---|
| 312 | // Rotate around the x- and y-axis
|
|---|
| 313 | v1.RotateY(sign*fAn);
|
|---|
| 314 | v1.RotateX(sign*fAw);
|
|---|
| 315 |
|
|---|
| 316 | // cout << "v1: " << v1.Theta()*TMath::RadToDeg() << " " << v1.Phi()*TMath::RadToDeg() << endl;
|
|---|
| 317 | // cout << "v2: " << v2.Theta()*TMath::RadToDeg() << " " << v2.Theta()*TMath::RadToDeg() << endl;
|
|---|
| 318 |
|
|---|
| 319 | // cout << "dv: " << (v2.Theta()-v1.Theta())*TMath::RadToDeg() << " " << (v2.Phi()-v1.Phi())*TMath::RadToDeg() << endl;
|
|---|
| 320 |
|
|---|
| 321 | Double_t dalt = v1.Theta()-v2.Theta();
|
|---|
| 322 | Double_t daz = v1.Phi() -v2.Phi();
|
|---|
| 323 |
|
|---|
| 324 | //cout << dalt*TMath::RadToDeg() << " " << daz*TMath::RadToDeg() << endl;
|
|---|
| 325 |
|
|---|
| 326 | if (daz>TMath::Pi())
|
|---|
| 327 | daz -= TMath::TwoPi();
|
|---|
| 328 | if (daz<-TMath::Pi())
|
|---|
| 329 | daz += TMath::TwoPi();
|
|---|
| 330 |
|
|---|
| 331 | // if (daz>TMath::Pi()/2)
|
|---|
| 332 | // {
|
|---|
| 333 | // }
|
|---|
| 334 |
|
|---|
| 335 | AltAz d(dalt, daz);
|
|---|
| 336 | return d;
|
|---|
| 337 |
|
|---|
| 338 | // Calculate Delta Azimuth and Delta Elevation
|
|---|
| 339 | /*
|
|---|
| 340 | AltAz d(TMath::Pi()/2-v1.Theta(), v1.Phi());
|
|---|
| 341 |
|
|---|
| 342 | cout << "p : " << p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
|
|---|
| 343 | cout << "d : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
|
|---|
| 344 | d -= p;
|
|---|
| 345 | cout << "d-p: " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
|
|---|
| 346 | d *= sign;
|
|---|
| 347 | cout << "d* : " << d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
|
|---|
| 348 |
|
|---|
| 349 |
|
|---|
| 350 | cout << "p2: " << 90-p.Alt()*TMath::RadToDeg() << " " << p.Az()*TMath::RadToDeg() << endl;
|
|---|
| 351 | cout << "d2: " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
|
|---|
| 352 |
|
|---|
| 353 | Int_t s1 = 90-d.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
|
|---|
| 354 | Int_t s2 = 90-p.Alt()*TMath::RadToDeg() < 0 ? -1 : 1;
|
|---|
| 355 |
|
|---|
| 356 |
|
|---|
| 357 | if (s1 != s2)
|
|---|
| 358 | {
|
|---|
| 359 | //90-d.Alt() <-- -90+d.Alt()
|
|---|
| 360 |
|
|---|
| 361 | d.Alt(d.Alt()-TMath::Pi());
|
|---|
| 362 | cout << "Alt-" << endl;
|
|---|
| 363 | }
|
|---|
| 364 | cout << "d': " << 90-d.Alt()*TMath::RadToDeg() << " " << d.Az()*TMath::RadToDeg() << endl;
|
|---|
| 365 | /*
|
|---|
| 366 | // Fix 'direction' of output depending on input vector
|
|---|
| 367 | if (TMath::Pi()/2-sign*p.Alt()<0)
|
|---|
| 368 | {
|
|---|
| 369 | d.Alt(d.Alt()-TMath::Pi());
|
|---|
| 370 | cout << "Alt-" << endl;
|
|---|
| 371 | }
|
|---|
| 372 | //if (TMath::Pi()/2-sign*p.Alt()>TMath::Pi())
|
|---|
| 373 | //{
|
|---|
| 374 | // d.Alt(TMath::Pi()-d.Alt());
|
|---|
| 375 | // cout << "Alt+" << endl;
|
|---|
| 376 | //}
|
|---|
| 377 |
|
|---|
| 378 | // Align correction into [-180,180]
|
|---|
| 379 | while (d.Az()>TMath::Pi())
|
|---|
| 380 | {
|
|---|
| 381 | d.Az(d.Az()-TMath::Pi()*2);
|
|---|
| 382 | cout << "Az-" << endl;
|
|---|
| 383 | }
|
|---|
| 384 | while (d.Az()<-TMath::Pi())
|
|---|
| 385 | {
|
|---|
| 386 | d.Az(d.Az()+TMath::Pi()*2);
|
|---|
| 387 | cout << "Az+" << endl;
|
|---|
| 388 | }
|
|---|
| 389 | */
|
|---|
| 390 | return d;
|
|---|
| 391 | }
|
|---|
| 392 |
|
|---|
| 393 |
|
|---|
| 394 | AltAz MBending::Correct(const AltAz &aa) const
|
|---|
| 395 | {
|
|---|
| 396 | // Correct [rad]
|
|---|
| 397 | // zdaz [rad]
|
|---|
| 398 | AltAz p = aa;
|
|---|
| 399 |
|
|---|
| 400 | DEBUG(cout << setprecision(16));
|
|---|
| 401 | DEBUG(cout << "Bend7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 402 |
|
|---|
| 403 | const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
|
|---|
| 404 | const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
|
|---|
| 405 | p += CRX;
|
|---|
| 406 | p += CRY;
|
|---|
| 407 |
|
|---|
| 408 | DEBUG(cout << "Bend6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 409 |
|
|---|
| 410 | const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
|
|---|
| 411 | const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
|
|---|
| 412 | p += NRX;
|
|---|
| 413 | p += NRY;
|
|---|
| 414 |
|
|---|
| 415 | DEBUG(cout << "Bend5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 416 |
|
|---|
| 417 | const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
|
|---|
| 418 | const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
|
|---|
| 419 | p += CES;
|
|---|
| 420 | p += CEC;
|
|---|
| 421 |
|
|---|
| 422 | DEBUG(cout << "Bend4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 423 |
|
|---|
| 424 | const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
|
|---|
| 425 | const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
|
|---|
| 426 | //p += TX;
|
|---|
| 427 | p += TF;
|
|---|
| 428 |
|
|---|
| 429 |
|
|---|
| 430 | DEBUG(cout << "Bend3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 431 |
|
|---|
| 432 | /*
|
|---|
| 433 | //New Corrections for NPAE and CA:
|
|---|
| 434 | TVector3 v(1.,1.,1.); // Vector in cartesian coordinates
|
|---|
| 435 |
|
|---|
| 436 | //Set Azimuth and Elevation
|
|---|
| 437 | v.SetPhi(p.Az());
|
|---|
| 438 | v.SetTheta(TMath::Pi()/2-p.Alt());
|
|---|
| 439 | //Rotation Vectors:
|
|---|
| 440 | TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
|
|---|
| 441 | TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
|
|---|
| 442 | //Rotate around the vectors vNpae and vCa
|
|---|
| 443 | v.Rotate(fNpae, vNpae);
|
|---|
| 444 | v.Rotate(fCa, vCa);
|
|---|
| 445 |
|
|---|
| 446 | p.Az(v.Phi());
|
|---|
| 447 | p.Alt(TMath::Pi()/2-v.Theta());
|
|---|
| 448 | */
|
|---|
| 449 |
|
|---|
| 450 | //Old correction terms for Npae and Ca:
|
|---|
| 451 | const AltAz CA(0, -fCa/cos(p.Alt()));
|
|---|
| 452 | p += CA;
|
|---|
| 453 |
|
|---|
| 454 | const AltAz NPAE(0, -fNpae*tan(p.Alt()));
|
|---|
| 455 | p += NPAE;
|
|---|
| 456 |
|
|---|
| 457 | DEBUG(cout << "Bend2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 458 |
|
|---|
| 459 | const AltAz ANAW(CalcAnAw(p, -1));
|
|---|
| 460 | p += ANAW;
|
|---|
| 461 |
|
|---|
| 462 | /* Old correction terms for An and Aw:
|
|---|
| 463 | const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
|
|---|
| 464 | const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
|
|---|
| 465 | p += AW;
|
|---|
| 466 | p += AN;
|
|---|
| 467 | */
|
|---|
| 468 |
|
|---|
| 469 | DEBUG(cout << "Bend1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 470 |
|
|---|
| 471 | const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
|
|---|
| 472 | p += FLOP;
|
|---|
| 473 |
|
|---|
| 474 | const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
|
|---|
| 475 | p += MAGIC1;
|
|---|
| 476 |
|
|---|
| 477 | const AltAz I(fIe, fIa);
|
|---|
| 478 | p += I;
|
|---|
| 479 |
|
|---|
| 480 | DEBUG(cout << "Bend0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 481 |
|
|---|
| 482 | return p;
|
|---|
| 483 | }
|
|---|
| 484 |
|
|---|
| 485 | AltAz MBending::CorrectBack(const AltAz &aa) const
|
|---|
| 486 | {
|
|---|
| 487 | // Correct [rad]
|
|---|
| 488 | // zdaz [rad]
|
|---|
| 489 | AltAz p = aa;
|
|---|
| 490 |
|
|---|
| 491 | DEBUG(cout << setprecision(16));
|
|---|
| 492 | DEBUG(cout << "Back0: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 493 |
|
|---|
| 494 | const AltAz I(fIe, fIa);
|
|---|
| 495 | p -= I;
|
|---|
| 496 |
|
|---|
| 497 | const AltAz MAGIC1(fMagic1*TMath::Sign(1., sin(p.Az())), 0);
|
|---|
| 498 | p -= MAGIC1;
|
|---|
| 499 |
|
|---|
| 500 | //const AltAz MAGIC1(fMagic1*sin(p.Az()), 0);
|
|---|
| 501 | //p -= MAGIC1;
|
|---|
| 502 |
|
|---|
| 503 | const AltAz FLOP(Sign(fFlop, p.Alt()), 0);
|
|---|
| 504 | p -= FLOP;
|
|---|
| 505 |
|
|---|
| 506 | DEBUG(cout << "Back1: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 507 |
|
|---|
| 508 | /* Old correction terms for An and Aw:
|
|---|
| 509 | const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt()));
|
|---|
| 510 | const AltAz AW( fAw*sin(p.Az()), -fAw*cos(p.Az())*tan(p.Alt()));
|
|---|
| 511 | p -= AN;
|
|---|
| 512 | p -= AW;
|
|---|
| 513 | */
|
|---|
| 514 |
|
|---|
| 515 | const AltAz ANAW(CalcAnAw(p, -1));
|
|---|
| 516 | p -= ANAW;
|
|---|
| 517 |
|
|---|
| 518 | DEBUG(cout << "Back2: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 519 |
|
|---|
| 520 | //Old Correction terms for Npae and Ca:
|
|---|
| 521 | const AltAz NPAE(0, -fNpae*tan(p.Alt()));
|
|---|
| 522 | p -= NPAE;
|
|---|
| 523 |
|
|---|
| 524 | const AltAz CA(0, -fCa/cos(p.Alt()));
|
|---|
| 525 | p -= CA;
|
|---|
| 526 |
|
|---|
| 527 | /*
|
|---|
| 528 | //New Correction term for Npae and Ca:
|
|---|
| 529 | TVector3 v2(1.,1.,1.); // Vector in cartesian coordinates
|
|---|
| 530 | //Set Azimuth and Elevation
|
|---|
| 531 | v2.SetPhi(p.Az());
|
|---|
| 532 | v2.SetTheta(TMath::Pi()/2-p.Alt());
|
|---|
| 533 | //Rotation Vectors:
|
|---|
| 534 | TVector3 vNpae( cos(p.Az()), sin(p.Az()), 0);
|
|---|
| 535 | TVector3 vCa( -cos(p.Az())*cos(p.Alt()), -sin(p.Az())*cos(p.Alt()), sin(p.Alt()));
|
|---|
| 536 | //Rotate around the vectors vCa and vNpae
|
|---|
| 537 | v2.Rotate(-fCa, vCa);
|
|---|
| 538 | v2.Rotate(-fNpae, vNpae);
|
|---|
| 539 |
|
|---|
| 540 | p.Az(v2.Phi());
|
|---|
| 541 | p.Alt(TMath::Pi()/2-v2.Theta());
|
|---|
| 542 | */
|
|---|
| 543 |
|
|---|
| 544 | DEBUG(cout << "Back3: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 545 |
|
|---|
| 546 | const AltAz TF(Sign(fTf*cos(p.Alt()), p.Alt()), 0);
|
|---|
| 547 | const AltAz TX(Sign(fTx/tan(p.Alt()), p.Alt()), 0);
|
|---|
| 548 | p -= TF;
|
|---|
| 549 | //p -= TX;
|
|---|
| 550 |
|
|---|
| 551 | DEBUG(cout << "Back4: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 552 |
|
|---|
| 553 | const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az()));
|
|---|
| 554 | const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az()));
|
|---|
| 555 | p -= CEC;
|
|---|
| 556 | p -= CES;
|
|---|
| 557 |
|
|---|
| 558 | DEBUG(cout << "Back5: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 559 |
|
|---|
| 560 | const AltAz NRY(fNry*cos(p.Alt()), -fNry*tan(p.Alt()));
|
|---|
| 561 | const AltAz NRX(fNrx*sin(p.Alt()), -fNrx);
|
|---|
| 562 | p -= NRY;
|
|---|
| 563 | p -= NRX;
|
|---|
| 564 |
|
|---|
| 565 | DEBUG(cout << "Back6: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 566 |
|
|---|
| 567 | const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt()));
|
|---|
| 568 | const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt()));
|
|---|
| 569 | p -= CRY;
|
|---|
| 570 | p -= CRX;
|
|---|
| 571 |
|
|---|
| 572 | DEBUG(cout << "Back7: " << 90-p.Alt()*180/TMath::Pi() << " " << p.Az()*180/TMath::Pi() << endl);
|
|---|
| 573 |
|
|---|
| 574 | return p;
|
|---|
| 575 | }
|
|---|
| 576 |
|
|---|
| 577 | ZdAz MBending::Correct(const ZdAz &zdaz) const
|
|---|
| 578 | {
|
|---|
| 579 | // Correct [rad]
|
|---|
| 580 | // zdaz [rad]
|
|---|
| 581 | AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
|
|---|
| 582 | AltAz c = Correct(p);
|
|---|
| 583 | return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
|
|---|
| 584 | }
|
|---|
| 585 |
|
|---|
| 586 | ZdAz MBending::CorrectBack(const ZdAz &zdaz) const
|
|---|
| 587 | {
|
|---|
| 588 | // Correct [rad]
|
|---|
| 589 | // zdaz [rad]
|
|---|
| 590 | AltAz p(TMath::Pi()/2-zdaz.Zd(), zdaz.Az());
|
|---|
| 591 | AltAz c = CorrectBack(p);
|
|---|
| 592 | return ZdAz(TMath::Pi()/2-c.Alt(), c.Az());
|
|---|
| 593 | }
|
|---|
| 594 |
|
|---|
| 595 | void MBending::SetParameters(const Double_t *par, Int_t n)
|
|---|
| 596 | {
|
|---|
| 597 | Clear();
|
|---|
| 598 |
|
|---|
| 599 | while (n--)
|
|---|
| 600 | *fCoeff[n] = par[n]/kRad2Deg;
|
|---|
| 601 | }
|
|---|
| 602 |
|
|---|
| 603 | void MBending::GetParameters(Double_t *par, Int_t n) const
|
|---|
| 604 | {
|
|---|
| 605 | while (n--)
|
|---|
| 606 | par[n] = *fCoeff[n]*kRad2Deg;
|
|---|
| 607 | }
|
|---|
| 608 |
|
|---|
| 609 | void MBending::SetMinuitParameters(TMinuit &m, Int_t n) const
|
|---|
| 610 | {
|
|---|
| 611 | if (n<0)
|
|---|
| 612 | n = fNumPar;
|
|---|
| 613 |
|
|---|
| 614 | Int_t ierflg = 0;
|
|---|
| 615 |
|
|---|
| 616 | while (n--)
|
|---|
| 617 | m.mnparm(n, fName[n], *fCoeff[n]*kRad2Deg, 1, -360, 360, ierflg);
|
|---|
| 618 | }
|
|---|
| 619 |
|
|---|
| 620 | void MBending::GetMinuitParameters(TMinuit &m, Int_t n)
|
|---|
| 621 | {
|
|---|
| 622 | if (n<0 || n>m.GetNumPars())
|
|---|
| 623 | n = m.GetNumPars();
|
|---|
| 624 |
|
|---|
| 625 | while (n--)
|
|---|
| 626 | {
|
|---|
| 627 | m.GetParameter(n, *fCoeff[n], fError[n]);
|
|---|
| 628 | *fCoeff[n] /= kRad2Deg;
|
|---|
| 629 | fError[n] /= kRad2Deg;
|
|---|
| 630 | }
|
|---|
| 631 | }
|
|---|
| 632 | /*
|
|---|
| 633 | void FormatPar(TMinuit &m, Int_t n)
|
|---|
| 634 | {
|
|---|
| 635 | Double_t par, err;
|
|---|
| 636 | m.GetParameter(n, par, err);
|
|---|
| 637 |
|
|---|
| 638 | int expp = (int)log10(par);
|
|---|
| 639 | int expe = (int)log10(err);
|
|---|
| 640 |
|
|---|
| 641 | if (err<2*pow(10, expe))
|
|---|
| 642 | expe--;
|
|---|
| 643 |
|
|---|
| 644 | Int_t exp = expe>expp ? expp : expe;
|
|---|
| 645 |
|
|---|
| 646 | par = (int)(par/pow(10, exp)) * pow(10, exp);
|
|---|
| 647 | err = (int)(err/pow(10, exp)) * pow(10, exp);
|
|---|
| 648 |
|
|---|
| 649 | cout << par << " +- " << err << flush;
|
|---|
| 650 | }
|
|---|
| 651 | */
|
|---|
| 652 | void MBending::PrintMinuitParameters(TMinuit &m, Int_t n) const
|
|---|
| 653 | {
|
|---|
| 654 | if (n<0)
|
|---|
| 655 | n = m.GetNumPars();
|
|---|
| 656 |
|
|---|
| 657 | cout << setprecision(3);
|
|---|
| 658 |
|
|---|
| 659 | Double_t par, err;
|
|---|
| 660 |
|
|---|
| 661 | while (n--)
|
|---|
| 662 | {
|
|---|
| 663 | m.GetParameter(n, par, err);
|
|---|
| 664 | cout << Form(" %2d %6s: ", n, (const char*)fName[n]);
|
|---|
| 665 | cout << setw(8) << par << " \xb1 " << setw(6) << err << endl;
|
|---|
| 666 | }
|
|---|
| 667 | }
|
|---|