Index: trunk/MagicSoft/Cosy/main/MBending.cc =================================================================== --- trunk/MagicSoft/Cosy/main/MBending.cc (revision 1803) +++ 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 1803) +++ 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;