Index: /trunk/MagicSoft/Cosy/main/MBending.cc =================================================================== --- /trunk/MagicSoft/Cosy/main/MBending.cc (revision 1804) +++ /trunk/MagicSoft/Cosy/main/MBending.cc (revision 1805) @@ -2,7 +2,9 @@ #include +#include #include +#include "timer.h" ClassImp(MBending); @@ -10,18 +12,5 @@ void MBending::Reset() { - fIa = 0; - fIe = 0; - fCa = 0; - fAn = 0; - fAw = 0; - fNrx = 0; - fNry = 0; - fCrx = 0; - fCry = 0; - fNpae = 0; - fEces = 0; - fAces = 0; - fEcec = 0; - fAcec = 0; + Clear(); } @@ -104,18 +93,20 @@ val *= kDeg2Rad; - if (str=="IA") fIa = val; - if (str=="IE") fIe = val; - if (str=="CA") fCa = val; - if (str=="AN") fAn = val; - if (str=="AW") fAw = val; - if (str=="NRX") fNrx = val; - if (str=="NRY") fNry = val; - if (str=="CRX") fCrx = val; - if (str=="CRY") fCry = val; - if (str=="NPAE") fNpae = val; - if (str=="ECES") fEces = val; - if (str=="ACES") fAces = val; - if (str=="ECEC") fEcec = val; - if (str=="ACEC") fAcec = val; + if (str=="IA") fIa = val; + if (str=="IE") fIe = val; + if (str=="CA") fCa = val; + if (str=="AN") fAn = val; + if (str=="AW") fAw = val; + if (str=="NRX") fNrx = val; + if (str=="NRY") fNry = val; + if (str=="CRX") fCrx = val; + if (str=="CRY") fCry = val; + if (str=="NPAE") fNpae = val; + if (str=="ECES") fEces = val; + if (str=="ACES") fAces = val; + if (str=="ECEC") fEcec = val; + if (str=="ACEC") fAcec = val; + if (str=="MAGIC1") fMagic1 = val; + if (str=="MAGIC2") fMagic2 = val; fin >> val; @@ -125,5 +116,49 @@ } -void MBending::Save(const char *name) {} +void MBending::Save(const char *name) +{ + /* + ! MMT 1987 July 8 + ! T 36 7.3622 41.448 -0.0481 + ! IA -37.5465 20.80602 + ! IE -13.9180 1.25217 + ! NPAE +7.0751 26.44763 + ! CA -6.9149 32.05358 + ! AN +0.5053 1.40956 + ! AW -2.2016 1.37480 + ! END + */ + + ofstream fout(name); + if (!fout) + { + cout << "Error: Cannot open file '" << name << "'" << endl; + return; + } + + Timer t; + t.Now(); + + fout << "MAGIC1 " << t.GetTimeStr() << endl; + fout << "S 00 000000 000000 0000000" << endl; + fout << setprecision(8); + fout << " IA " << kRad2Deg*fIa << " -1" << endl; + fout << " IE " << kRad2Deg*fIe << " -1" << endl; + fout << " CA " << kRad2Deg*fNpae << " -1" << endl; + fout << " NPAE " << kRad2Deg*fCa << " -1" << endl; + fout << " AN " << kRad2Deg*fAn << " -1" << endl; + fout << " AW " << kRad2Deg*fAw << " -1" << endl; + fout << " NRX " << kRad2Deg*fNrx << " -1" << endl; + fout << " NRY " << kRad2Deg*fNry << " -1" << endl; + fout << " CRX " << kRad2Deg*fCrx << " -1" << endl; + fout << " CRY " << kRad2Deg*fCry << " -1" << endl; + fout << " ECES " << kRad2Deg*fEces << " -1" << endl; + fout << " ACES " << kRad2Deg*fAces << " -1" << endl; + fout << " ECEC " << kRad2Deg*fEcec << " -1" << endl; + fout << " ACEC " << kRad2Deg*fAcec << " -1" << endl; + fout << " MAGIC1 " << kRad2Deg*fMagic1 << " -1" << endl; + fout << " MAGIC2 " << kRad2Deg*fMagic2 << " -1" << endl; + fout << "END" << endl; +} AltAz MBending::Correct(const AltAz &aa) const @@ -133,11 +168,11 @@ AltAz p = aa; - const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az())); - const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az())); + const AltAz CES(-fEces*sin(p.Alt()), -fAces*sin(p.Az())); + const AltAz CEC(-fEcec*cos(p.Alt()), -fAcec*cos(p.Az())); p += CES; p += CEC; - const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()), fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt())); - const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt())); + const AltAz CRX(-fCrx*sin(p.Az()-p.Alt()), fCrx*cos(p.Az()-p.Alt())/cos(p.Alt())); + const AltAz CRY(-fCry*cos(p.Az()-p.Alt()), -fCry*sin(p.Az()-p.Alt())/cos(p.Alt())); p += CRX; p += CRY; @@ -148,16 +183,22 @@ p += NRY; - const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt())); - const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt())); + const AltAz MAGIC(-fMagic1*sin(p.Az()-fMagic2), 0); + p += MAGIC; + + const AltAz AW( fAw*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt())); + const AltAz AN(-fAn*cos(p.Az()), -fAn*sin(p.Az())*tan(p.Alt())); p += AW; p += AN; - const AltAz CA(0, -fCa*kDeg2Rad/cos(p.Alt())); + //zd1 += 11.22; + //zd1 += 1.15*sin((az0+0.05)*3.1415/180); + + const AltAz CA(0, -fCa/cos(p.Alt())); p += CA; - const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt())); + const AltAz NPAE(0, -fNpae*tan(p.Alt())); p += NPAE; - const AltAz I(fIe/**kDeg2Rad*/, fIa/**kDeg2Rad*/); + const AltAz I(fIe, fIa); p += I; @@ -171,32 +212,32 @@ AltAz p = aa; - const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az())); - const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az())); - p -= CES; - p -= CEC; - - const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()), fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt())); - const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt())); - p -= CRX; - p -= CRY; - - const AltAz NRX(fNrx*kDeg2Rad*sin(p.Alt()), -fNrx*kDeg2Rad); - const AltAz NRY(fNry*kDeg2Rad*cos(p.Alt()), -fNry*kDeg2Rad*tan(p.Alt())); - p -= NRX; - p -= NRY; - - const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), -fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt())); - const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt())); - p -= AW; - p -= AN; + const AltAz I(fIe/**kDeg2Rad*/, fIa/**kDeg2Rad*/); + p -= I; + + const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt())); + p -= NPAE; const AltAz CA(0, -fCa*kDeg2Rad/cos(p.Alt())); p -= CA; - const AltAz NPAE(0, -fNpae*kDeg2Rad*tan(p.Alt())); - p -= NPAE; - - const AltAz I(fIe/**kDeg2Rad*/, fIa/**kDeg2Rad*/); - p -= I; + const AltAz AN(-fAn*kDeg2Rad*cos(p.Az()), -fAn*kDeg2Rad*sin(p.Az())*tan(p.Alt())); + const AltAz AW( fAw*kDeg2Rad*sin(p.Az()), 0/*-fAw*kDeg2Rad*cos(p.Az())*tan(p.Alt())*/); + p -= AN; + p -= AW; + + const AltAz NRY(fNry*kDeg2Rad*cos(p.Alt()), -fNry*kDeg2Rad*tan(p.Alt())); + const AltAz NRX(fNrx*kDeg2Rad*sin(p.Alt()), -fNrx*kDeg2Rad); + p -= NRY; + p -= NRX; + + const AltAz CRY(-fCry*kDeg2Rad*cos(p.Az()-p.Alt()), -fCry*kDeg2Rad*sin(p.Az()-p.Alt())/cos(p.Alt())); + const AltAz CRX(-fCrx*kDeg2Rad*sin(p.Az()-p.Alt()), fCrx*kDeg2Rad*cos(p.Az()-p.Alt())/cos(p.Alt())); + p -= CRY; + p -= CRX; + + const AltAz CEC(-fEcec*kDeg2Rad*cos(p.Alt()), -fAcec*kDeg2Rad*cos(p.Az())); + const AltAz CES(-fEces*kDeg2Rad*sin(p.Alt()), -fAces*kDeg2Rad*sin(p.Az())); + p -= CEC; + p -= CES; return p; @@ -221,5 +262,5 @@ } -void MBending::SetParameters(const Double_t *par, Int_t n=14) +void MBending::SetParameters(const Double_t *par, Int_t n) { Clear(); @@ -227,67 +268,75 @@ switch (n) { + case 16: + fMagic2=par[15]/kRad2Deg; // Magic User Defined + case 15: + fMagic1=par[14]/kRad2Deg; // Magic User Defined case 14: - fAcec=par[13] ; // Azimuth Centering Error (cos) + fAcec =par[13]/kRad2Deg; // Azimuth Centering Error (cos) case 13: - fEcec=par[12] ; // Elevation Centering Error (cos) + fEcec =par[12]/kRad2Deg; // Elevation Centering Error (cos) case 12: - fAces=par[11] ; // Azimuth Centering Error (sin) + fAces =par[11]/kRad2Deg; // Azimuth Centering Error (sin) case 11: - fEces=par[10] ; // Elevation Centering Error (sin) + fEces =par[10]/kRad2Deg; // Elevation Centering Error (sin) case 10: - fCry =par[9] ; // Alt/Az Coude Displacement (E-W) + fCry =par[9]/kRad2Deg; // Alt/Az Coude Displacement (E-W) case 9: - fCrx =par[8] ; // Alt/Az Coude Displacement (N-S) + fCrx =par[8]/kRad2Deg; // Alt/Az Coude Displacement (N-S) case 8: - fNry =par[7] ; // Nasmyth rotator displacement, vertical + fNry =par[7]/kRad2Deg; // Nasmyth rotator displacement, vertical case 7: - fNrx =par[6] ; // Nasmyth rotator displacement, horizontan + fNrx =par[6]/kRad2Deg; // Nasmyth rotator displacement, horizontan case 6: - fAw =par[5] ; // Azimuth Axis Misalignment (E-W) + fAw =par[5]/kRad2Deg; // Azimuth Axis Misalignment (E-W) case 5: - fAn =par[4] ; // Azimuth Axis Misalignment (N-S) + fAn =par[4]/kRad2Deg; // Azimuth Axis Misalignment (N-S) case 4: - fCa =par[3] ; // Left-Right Collimation Error + fCa =par[3]/kRad2Deg; // Left-Right Collimation Error case 3: - fNpae=par[2] ; // Az-El Nonperpendicularity + fNpae =par[2]/kRad2Deg; // Az-El Nonperpendicularity case 2: - fIe =par[1] ; // Index Error in Elevation + fIe =par[1]/kRad2Deg; // Index Error in Elevation case 1: - fIa =par[0] ; // Index Error in Azimuth - } -} - -void MBending::GetParameters(Double_t *par, Int_t n=14) const + fIa =par[0]/kRad2Deg; // Index Error in Azimuth + } +} + +void MBending::GetParameters(Double_t *par, Int_t n) const { switch (n) { + case 16: + par[15]=fMagic2*kRad2Deg; // + case 15: + par[14]=fMagic1*kRad2Deg; // case 14: - par[13]=fAcec; // Azimuth Centering Error (cos) + par[13]=fAcec*kRad2Deg; // Azimuth Centering Error (cos) case 13: - par[12]=fEcec; // Elevation Centering Error (cos) + par[12]=fEcec*kRad2Deg; // Elevation Centering Error (cos) case 12: - par[11]=fAces ; // Azimuth Centering Error (sin) + par[11]=fAces*kRad2Deg; // Azimuth Centering Error (sin) case 11: - par[10]=fEces ; // Elevation Centering Error (sin) + par[10]=fEces*kRad2Deg; // Elevation Centering Error (sin) case 10: - par[9]=fCry ; // Alt/Az Coude Displacement (E-W) + par[9]=fCry*kRad2Deg; // Alt/Az Coude Displacement (E-W) case 9: - par[8]=fCrx ; // Alt/Az Coude Displacement (N-S) + par[8]=fCrx*kRad2Deg; // Alt/Az Coude Displacement (N-S) case 8: - par[7]=fNry ; // Nasmyth rotator displacement, vertical + par[7]=fNry*kRad2Deg; // Nasmyth rotator displacement, vertical case 7: - par[6]=fNrx ; // Nasmyth rotator displacement, horizontan + par[6]=fNrx*kRad2Deg; // Nasmyth rotator displacement, horizontan case 6: - par[5]=fAw ; // Azimuth Axis Misalignment (E-W) + par[5]=fAw*kRad2Deg; // Azimuth Axis Misalignment (E-W) case 5: - par[4]=fAn ; // Azimuth Axis Misalignment (N-S) + par[4]=fAn*kRad2Deg; // Azimuth Axis Misalignment (N-S) case 4: - par[3]=fCa ; // Left-Right Collimation Error + par[3]=fCa*kRad2Deg; // Left-Right Collimation Error case 3: - par[2]=fNpae ; // Az-El Nonperpendicularity + par[2]=fNpae*kRad2Deg; // Az-El Nonperpendicularity case 2: - par[1]=fIe ; // Index Error in Elevation + par[1]=fIe*kRad2Deg; // Index Error in Elevation case 1: - par[0]=fIa ; // Index Error in Azimuth + par[0]=fIa*kRad2Deg; // Index Error in Azimuth } } @@ -302,24 +351,51 @@ switch (n) { + case 16: + m.mnparm(15,"MAGIC2", fMagic2*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; + case 15: + m.mnparm(14,"MAGIC1", fMagic1*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 14: + m.mnparm(13,"ACEC", fAcec*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 13: + m.mnparm(12,"ECEC", fEcec*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 12: + m.mnparm(11,"ACES", fAcec*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 11: + m.mnparm(10,"ECES", fEcec*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 10: + m.mnparm(9, "CRY", fCry*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 9: + m.mnparm(8, "CRX", fCrx*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 8: + m.mnparm(7, "NRY", fNry*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 7: + m.mnparm(6, "NRX", fNrx*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; + case 6: + m.mnparm(5, "AW", fAw*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 5: + m.mnparm(4, "AN", fAn*kRad2Deg, 1, -360, 360, ierflg); + // cout << "Init 3 CA: " << fCa << endl; case 4: - m.mnparm(3, "CA", fCa, 1, -360, 360, ierflg); + m.mnparm(3, "CA", fCa*kRad2Deg, 1, -360, 360, ierflg); // cout << "Init 3 CA: " << fCa << endl; case 3: - m.mnparm(2, "NPAE", fNpae, 1, -360, 360, ierflg); + m.mnparm(2, "NPAE", fNpae*kRad2Deg, 1, -360, 360, ierflg); // cout << "Init 2 NPAE: " << fNpae << endl; case 2: - m.mnparm(1, "IE", fIe, 1, -360, 360, ierflg); + m.mnparm(1, "IE", fIe*kRad2Deg, 1, -360, 360, ierflg); // cout << "Init 1 IE: " << fIe << endl; case 1: - m.mnparm(0, "IA", fIa, 1, -360, 360, ierflg); + m.mnparm(0, "IA", fIa*kRad2Deg, 1, -360, 360, ierflg); // cout << "Init 0 IA: " << fIa << endl; } @@ -335,22 +411,52 @@ switch (n) { + case 16: + m.GetParameter(15, fMagic2, err); + fMagic2 /= kRad2Deg; + case 15: + m.GetParameter(14, fMagic1, err); + fMagic1 /= kRad2Deg; case 14: + m.GetParameter(13, fAcec, err); + fAcec /= kRad2Deg; case 13: + m.GetParameter(12, fEcec, err); + fEcec /= kRad2Deg; case 12: + m.GetParameter(11, fAces, err); + fAces /= kRad2Deg; case 11: + m.GetParameter(10, fEces, err); + fEces /= kRad2Deg; case 10: + m.GetParameter(9, fCry, err); + fCry /= kRad2Deg; case 9: + m.GetParameter(8, fCrx, err); + fCrx /= kRad2Deg; case 8: + m.GetParameter(7, fNry, err); + fNry /= kRad2Deg; case 7: + m.GetParameter(6, fNrx, err); + fNrx /= kRad2Deg; case 6: + m.GetParameter(5, fAw, err); + fAw /= kRad2Deg; case 5: + m.GetParameter(4, fAn, err); + fAn /= kRad2Deg; case 4: m.GetParameter(3, fCa, err); + fCa /= kRad2Deg; case 3: m.GetParameter(2, fNpae, err); + fNpae /= kRad2Deg; case 2: m.GetParameter(1, fIe, err); + fIe /= kRad2Deg; case 1: m.GetParameter(0, fIa, err); + fIa /= kRad2Deg; } } @@ -365,26 +471,52 @@ switch (n) { + case 16: + m.GetParameter(15, par, err); + cout << " 15 MAGIC2: " << par << " +- " << err << endl; + case 15: + m.GetParameter(14, par, err); + cout << " 14 MAGIC1: " << par << " +- " << err << endl; case 14: + m.GetParameter(13, par, err); + cout << " 13 ACEC: " << par << " +- " << err << endl; case 13: + m.GetParameter(12, par, err); + cout << " 12 ECEC: " << par << " +- " << err << endl; case 12: + m.GetParameter(11, par, err); + cout << " 11 ACES: " << par << " +- " << err << endl; case 11: + m.GetParameter(10, par, err); + cout << " 10 ECES: " << par << " +- " << err << endl; case 10: + m.GetParameter(9, par, err); + cout << " 9 CRY: " << par << " +- " << err << endl; case 9: + m.GetParameter(8, par, err); + cout << " 8 CRX: " << par << " +- " << err << endl; case 8: + m.GetParameter(7, par, err); + cout << " 7 NRY: " << par << " +- " << err << endl; case 7: + m.GetParameter(6, par, err); + cout << " 6 NRX: " << par << " +- " << err << endl; case 6: + m.GetParameter(5, par, err); + cout << " 5 AW: " << par << " +- " << err << endl; case 5: + m.GetParameter(4, par, err); + cout << " 4 AN: " << par << " +- " << err << endl; case 4: m.GetParameter(3, par, err); - cout << " 3 CA: " << par << " +- " << err << endl; + cout << " 3 CA: " << par << " +- " << err << endl; case 3: m.GetParameter(2, par, err); - cout << " 2 NPAE: " << par << " +- " << err << endl; + cout << " 2 NPAE: " << par << " +- " << err << endl; case 2: m.GetParameter(1, par, err); - cout << " 1 IE: " << par << " +- " << err << endl; + cout << " 1 IE: " << par << " +- " << err << endl; case 1: m.GetParameter(0, par, err); - cout << " 0 IA: " << par << " +- " << err << endl; - } -} + cout << " 0 IA: " << par << " +- " << err << endl; + } +} Index: /trunk/MagicSoft/Cosy/main/MBending.h =================================================================== --- /trunk/MagicSoft/Cosy/main/MBending.h (revision 1804) +++ /trunk/MagicSoft/Cosy/main/MBending.h (revision 1805) @@ -27,4 +27,7 @@ Double_t fEcec ; // [rad] Elevation Centering Error (cos) Double_t fAcec ; // [rad] Azimuth Centering Error (cos) + Double_t fMagic1; // [rad] Magic Term (what is it?) + Double_t fMagic2; // [rad] Magic Term (what is it?) + void Clear() @@ -44,4 +47,6 @@ fEcec=0 ; // Elevation Centering Error (cos) fAcec=0 ; // Azimuth Centering Error (cos) + fMagic1=0; // Azimuth Centering Error (cos) + fMagic2=0; // Azimuth Centering Error (cos) } @@ -65,5 +70,5 @@ ZdAz operator()(const ZdAz &zdaz, void (*fcn)(ZdAz &zdaz, Double_t *par)) const { - Double_t par[14]; + Double_t par[16]; GetParameters(par); ZdAz za = zdaz; @@ -74,5 +79,5 @@ AltAz operator()(const AltAz &aaz, void (*fcn)(AltAz &aaz, Double_t *par)) const { - Double_t par[14]; + Double_t par[16]; GetParameters(par); AltAz aa = aaz; @@ -80,8 +85,7 @@ return aa; } - - void SetParameters(const Double_t *par, Int_t n=14); - void GetParameters(Double_t *par, Int_t n=14) const; + void SetParameters(const Double_t *par, Int_t n=16); + void GetParameters(Double_t *par, Int_t n=16) const; void SetMinuitParameters(TMinuit &m, Int_t n=-1) const; Index: /trunk/MagicSoft/Cosy/tpoint/tpointfit.C =================================================================== --- /trunk/MagicSoft/Cosy/tpoint/tpointfit.C (revision 1805) +++ /trunk/MagicSoft/Cosy/tpoint/tpointfit.C (revision 1805) @@ -0,0 +1,453 @@ +#include +#include + +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include "coord.h" +#include "MBending.h" + +// +// Sekans = 1/cos +// + +void MyAdjust(AltAz &p, Double_t *par) +{ + // p[rad] + MBending bend; + bend.SetParameters(par); // [deg] + p=bend(p); +} + +class Set : public TObject +{ + friend ifstream &operator>>(ifstream &fin, Set &set); +private: + Double_t fOrigStarAz; + Double_t fOrigStarEl; + + Double_t fStarAz; + Double_t fStarEl; + + Double_t fRawAz; + Double_t fRawEl; + +public: + Set(Double_t sel=0, Double_t saz=0, Double_t rel=0, Double_t raz=0) : + fStarAz(saz*TMath::Pi()/180), + fStarEl(sel*TMath::Pi()/180), + fRawAz(raz*TMath::Pi()/180), + fRawEl(rel*TMath::Pi()/180) + { + } + + Double_t GetResidual() const + { + Double_t daz = fRawEl-fStarEl; + Double_t del = fRawAz-fStarAz; + + //Double_t rzd = TMath::Pi()/2-fRawEl; + Double_t rzd = TMath::Pi()/2-fOrigStarEl; + + Double_t dphi2 = daz/2.; + Double_t cos2 = cos(dphi2)*cos(dphi2); + Double_t sin2 = sin(dphi2)*sin(dphi2); + Double_t d = cos(del)*cos2 - cos(2*rzd+del)*sin2; + + Double_t dist = acos(d); + + return dist * 180 / TMath::Pi(); + } + + void operator=(Set &set) + { + fStarAz = set.fStarAz; + fStarEl = set.fStarEl; + fRawAz = set.fRawAz; + fRawEl = set.fRawEl; + } + + Double_t GetDEl() const { return (fRawEl-fStarEl)*180/TMath::Pi(); } + Double_t GetDAz() const { return (fRawAz-fStarAz)*180/TMath::Pi(); } + Double_t GetStarEl() const { return fStarEl*180/TMath::Pi(); } + Double_t GetStarZd() const { return 90.-fStarEl*180/TMath::Pi(); } + Double_t GetStarAz() const { return fStarAz*180/TMath::Pi(); } + Double_t GetRawEl() const { return fRawEl*180/TMath::Pi(); } + Double_t GetRawAz() const { return fRawAz*180/TMath::Pi(); } + Double_t GetRawZd() const { return 90.-fRawEl*180/TMath::Pi(); } + + ZdAz GetStarZdAz() const { return ZdAz(TMath::Pi()/2-fStarEl, fStarAz); } + AltAz GetStarAltAz() const { return AltAz(fStarEl, fStarAz); } + + ZdAz GetRawZdAz() const { return ZdAz(TMath::Pi()/2-fRawEl, fRawAz); } + AltAz GetRawAltAz() const { return AltAz(fRawEl, fRawAz); } + + void AdjustEl(Double_t del) { fStarEl += del*TMath::Pi()/180; } + void AdjustAz(Double_t daz) { fStarAz += daz*TMath::Pi()/180; } + + void Adjust(const MBending &bend) + { + AltAz p = bend(GetStarAltAz()); + fStarEl = p.Alt(); + fStarAz = p.Az(); + } + + void Adjust(const MBending &bend, void (*fcn)(AltAz &zdaz, Double_t *par)) + { + AltAz star = GetStarAltAz(); + + AltAz p = bend(star, fcn); + fStarEl = p.Alt(); + fStarAz = p.Az(); + } +}; + +ifstream &operator>>(ifstream &fin, Set &set) +{ + Double_t v[4]; + fin >> v[2]; + fin >> v[3]; + fin >> v[0]; + fin >> v[1]; + + set.fStarAz = v[0]*TMath::Pi()/180; + set.fStarEl = v[1]*TMath::Pi()/180; + + set.fOrigStarAz = v[0]*TMath::Pi()/180; + set.fOrigStarEl = v[1]*TMath::Pi()/180; + + set.fRawAz = v[2]*TMath::Pi()/180; + set.fRawEl = v[3]*TMath::Pi()/180; + + return fin; +} + +TList gCoordinates; + +void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) +{ + f = 0; + + MBending bend; + bend.SetParameters(par/*, npar*/); // Set Parameters [deg] to MBending + + for (int i=0; iGetView(); + + if (!view) + { + cout << "No View!" << endl; + return; + } + + TMarker mark0; + TMarker mark1; + mark0.SetMarkerStyle(kStar); + mark1.SetMarkerStyle(kStar); + mark1.SetMarkerColor(kRed); + + r0 /= 90; + r1 /= 90; + phi0 *= TMath::Pi()/180; + phi1 *= TMath::Pi()/180; + + Double_t x0[3] = { r0*cos(phi0), r0*sin(phi0), 0}; + Double_t x1[3] = { r1*cos(phi1), r1*sin(phi1), 0}; + + Double_t y0[3], y1[3]; + + view->WCtoNDC(x0, y0); + view->WCtoNDC(x1, y1); + + mark0.DrawMarker(y0[0], y0[1]); + mark1.DrawMarker(y1[0], y1[1]); +} + +void DrawPolLine(TVirtualPad *pad, Double_t r0, Double_t phi0, Double_t r1, Double_t phi1) +{ + TView *view = pad->GetView(); + + if (!view) + { + cout << "No View!" << endl; + return; + } + + if (r0<0) + { + r0 = -r0; + phi0 += 180; + } + if (r1<0) + { + r1 = -r1; + phi1 += 180; + } + + /* + phi0 = fmod(phi0+360, 360); + phi1 = fmod(phi1+360, 360); + */ + + TLine line; + line.SetLineWidth(2); + line.SetLineColor(kBlue); + + Double_t p0 = phi0phi1) + { + Double_t d = r1; + r1 = r0; + r0 = d; + } + + r0 /= 90; + r1 /= 90; + + Double_t dr = r1-r0; + Double_t dp = p1-p0; + + Double_t x0[3] = { r0*cos(p0*TMath::Pi()/180), r0*sin(p0*TMath::Pi()/180), 0}; + + for (double i=p0+10; ip1) + i=p1; + + Double_t r = dr/dp*(i-p0)+r0; + Double_t p = TMath::Pi()/180*i; + + Double_t x1[3] = { r*cos(p), r*sin(p), 0}; + + Double_t y0[3], y1[3]; + + view->WCtoNDC(x0, y0); + view->WCtoNDC(x1, y1); + + line.DrawLine(y0[0], y0[1], y1[0], y1[1]); + + x0[0] = x1[0]; + x0[1] = x1[1]; + } +} + +void DrawSet(TVirtualPad *pad, Set &set, Float_t scale=-1) +{ + Double_t r0 = set.GetRawZd(); + Double_t phi0 = set.GetRawAz(); + Double_t r1 = set.GetStarZd(); + Double_t phi1 = set.GetStarAz(); + + if (scale>0) + { + Double_t d = r1-r0; + if (d<0) d=-d; + r0 -= scale*d; + r1 += scale*d; + d = phi1-phi0; + if (d<0) d=-d; + phi0 -= scale*d; + phi1 += scale*d; + } + + DrawPolLine(pad, r0, phi0, r1, phi1); + DrawMarker(pad, r0, phi0, r1, phi1); +} + +void tpointfit() +{ + cout << "Starting..." << endl; + + gCoordinates.SetOwner(); + + TH2F hdaz("dAz", "\\Delta Az vs. Alt", 32, 0, 90, 32, -1, 1); + TH2F hdzd("dZd", "\\Delta Alt vs. Az", 32, 0, 360, 32, -1, 1); + + TH1F hres1("Res1", " Residuals before correction ", 10, 0, 180); + TH1F hres2("Res2", " Residuals after correction ", 10, 0, 1); + + TH2F h2res1("Res2D1", " Arb. Residuals before correction ", 32, 0, 360, 10, 0, 90); + TH2F h2res2("Res2D2", " Arb. Residuals after correction ", 32, 0, 360, 10, 0, 90); + + hdaz.SetXTitle("Zd [\\circ]"); + hdaz.SetYTitle("\\Delta Az [\\circ]"); + + hdzd.SetXTitle("Az [\\circ]"); + hdzd.SetYTitle("\\Delta Zd [\\circ]"); + + hres1.SetXTitle("\\Delta [\\circ]"); + hres1.SetYTitle("Counts"); + + hres2.SetXTitle("\\Delta [\\circ]"); + hres2.SetYTitle("Counts"); + + TGraph gdaz; + TGraph gdzd; + + ifstream fin("tpoint/tpoint2.txt"); + + while (fin && fin.get()!='\n'); + while (fin && fin.get()!='\n'); + while (fin && fin.get()!='\n'); + if (!fin) + { + cout << "File not found!" << endl; + return; + } + + while (1) + { + Set set; + + fin >> set; // Read data from file [deg], it is stored in [rad] + + if (!fin) + break; + + gCoordinates.Add(new Set(set)); + } + + cout << "Found " << gCoordinates.GetSize() << " sets of coordinates." << endl; + + TMinuit minuit(16); //initialize TMinuit with a maximum of 5 params + minuit.SetPrintLevel(-1); + minuit.SetFCN(fcn); + + MBending bending; + bending.SetMinuitParameters(minuit, 16); // Init Parameters [deg] + + Int_t ierflg = 0; + + // minuit.FixParameter(2); // NPAE + // minuit.FixParameter(3); // CA + // minuit.FixParameter(4); // AN + // minuit.FixParameter(5); // AW + minuit.FixParameter(6); // NRX + minuit.FixParameter(7); // NRY + minuit.FixParameter(8); // CRX + minuit.FixParameter(9); // CRY + // minuit.FixParameter(10); // ECES + // minuit.FixParameter(11); // ACES + minuit.FixParameter(12); // ECEC + minuit.FixParameter(13); // ACEC + //minuit.FixParameter(14); // MAGIC1 + minuit.FixParameter(15); // MAGIC2 + + //minuit.Command("SHOW PARAMETERS"); + ierflg = minuit.Migrad(); + cout << endl << "Migrad returns " << ierflg << endl << endl; + minuit.Command("SHOW LIMITS"); + + cout << endl; + + /* + minuit.FixParameter(0); + minuit.FixParameter(1); + minuit.Release(2); + + ierflg = minuit.Migrad(); + cout << endl << "Migrad returns " << ierflg << endl << endl; + */ + + bending.GetMinuitParameters(minuit); + bending.PrintMinuitParameters(minuit); + bending.Save("bending_magic.txt"); + + TList list; + list.SetOwner(); + + for (int i=0; iDivide(3,2); + + c1->cd(1); + hdaz.DrawCopy("cont"); + gdaz.DrawClone("*"); + + c1->cd(4); + hdzd.DrawCopy("cont"); + gdzd.DrawClone("*"); + + c1->cd(2); + hres1.SetLineColor(kRed); + hres1.DrawCopy(); + + c1->cd(5); + hres2.SetLineColor(kBlue); + hres2.DrawCopy(); + +// c1->cd(3); +// gPad->SetTheta(90); +// gPad->SetPhi(90); +// h2res1.DrawCopy("surf1pol"); +// gPad->Modified(); +// gPad->Update(); +// for (int i=0; icd(6); +// gPad->SetTheta(90); +// gPad->SetPhi(90); +// h2res1.DrawCopy("surf1pol"); +// gPad->Modified(); +// gPad->Update(); +// for (int i=0; i