Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3567) +++ trunk/MagicSoft/Mars/Changelog (revision 3568) @@ -18,12 +18,67 @@ -*-*- END OF LINE -*-*- + + 2004/03/22: Thomas Bretz + + * mpointing/MSrcPosCalc.[h,cc], MSrcPosCam.[h,cc]: + - added + + * mastro/MAstro.[h,cc]: + - added code to calculate rotationangle previously in MObservatory + - changed definition of rotation angle such, that it is now + 180deg if Ra and Az grid is parallel + + * mastro/MAstroCatalog.[h,cc]: + - fixes and enhancements to the display (such as misscalculated + number of grid lines, title display, etc) + - enhancements to the output + - generalized creation of grid - for further usage + + * mastro/MAstroSky2Local.[h,cc]: + - replaced calculation of rotation angle by the function in + MAstro + + * mastro/MObservatory.[h,cc]: + - small changes to Print output + - moved code for calculation of rotation angle to MAstro + + * mbase/MEvtLoop.cc: + - do not output number of events per second if no events processed + + * mbase/MParList.cc: + - updated some comments + + * mfileio/MCT1ReadAscii.cc, mfileio/MCT1ReadPreProc.cc, + mfileio/MReadRflFile.cc, mraw/MRawFileRead.cc, + mreport/MReportFileRead.cc: + - output error string if file cannot be opened + + * mfileio/MReadTree.cc: + - output name of chain which is scanned + + * mimage/MConcentration.cc: + - replaced loop by iterator + - removed obsolete (unused) variables + + * mimage/MHNewImagePar.[h,cc]: + - fixed display colors + + * mpointing/MPointingPos.[h,cc]: + - added member function to calculate rotation angle + - added comments + + * mpointing/Makefile: + - added include MAstro + + + 2004/03/19: Markus Gaug * mcalib/MHCalibrationChargePix.cc - - added some style sto the default Draw in order to see the + - added some style to the default Draw in order to see the label and axis titles better * mcalib/MHCalibrationChargeCam.[h,cc] - - store and display more informaiton on the average pxiels + - store and display more information on the average pxiels * mcalib/MCalibrationCam.cc @@ -42,6 +97,6 @@ - included MCalibrationQECam to be initialized - * mcalib/MCalibrationChargePix.[h,cc] - * mcalib/MCalibrationQEPix.[h,cc] + * mcalib/MCalibrationChargePix.[h,cc], + mcalib/MCalibrationQEPix.[h,cc]: - replace DefinePixId by SetPixId @@ -59,26 +114,25 @@ - remove uncommented piece of code - * msignal/MExtractSignal.cc - * msignal/MExtractSignal2.cc + * msignal/MExtractSignal.cc, msignal/MExtractSignal2.cc: - remove warning about pixel with low gain saturation, now in MBadPixelsPix - * mbadpixels/MBadPixelsPix.[h,cc] - * mcalib/MCalibrationChargeCam.cc + * mbadpixels/MBadPixelsPix.[h,cc], mcalib/MCalibrationChargeCam.cc: - added new flag: kDeviatingNumPhes * mcalib/MCalibrationChargePix.cc - - check for mean arr. time in last bin replaced by check in last two bins - - * mcalib/MCalibrationChargePix.[h,cc] - * mcalib/MCalibrationChargeCam.cc - * mcalib/MHCalibrationChargeCam.cc - - removed flag kHiGainFitted, kLoGainFitted, since they are available - from MBadPixelsPix - - * macros/calibration.C - * macros/calibrate_data.C + - check for mean arr. time in last bin replaced by check in last + two bins + + * mcalib/MCalibrationChargePix.[h,cc], + mcalib/MCalibrationChargeCam.cc, + mcalib/MHCalibrationChargeCam.cc: + - removed flag kHiGainFitted, kLoGainFitted, since they are + available from MBadPixelsPix + + * macros/calibration.C, macros/calibrate_data.C - a few flags from MCalibrationChargeCam::GetPixelContent were wrong, corrected them + Index: trunk/MagicSoft/Mars/mastro/MAstro.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstro.cc (revision 3568) @@ -310,2 +310,97 @@ return MTime(ut1).GetGmst(); } + +// -------------------------------------------------------------------------- +// +// RotationAngle +// +// calculates the angle for the rotation of the sky coordinate system +// with respect to the local coordinate system. This is identical +// to the rotation angle of the sky image in the camera. +// +// sinl [rad]: sine of observers latitude +// cosl [rad]: cosine of observers latitude +// theta [rad]: polar angle/zenith distance +// phi [rad]: rotation angle/azimuth +// +// Return sin/cos component of angle +// +// The convention is such, that the rotation angle is -pi/pi if +// right ascension and local rotation angle are counted in the +// same direction, 0 if counted in the opposite direction. +// +// (In other words: The rotation angle is 0 when the source culminates) +// +// Using vectors it can be done like: +// TVector3 v, p; +// v.SetMagThetaPhi(1, theta, phi); +// p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0); +// v = v.Cross(l)); +// v.RotateZ(-phi); +// v.Rotate(-theta) +// rho = TMath::ATan2(v(2), v(1)); +// +// For more information see TDAS 00-11, eqs. (18) and (20) +// +void MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) +{ + const Double_t sint = TMath::Sin(theta); + const Double_t cost = TMath::Cos(theta); + + const Double_t snlt = sinl*sint; + const Double_t cslt = cosl*cost; + + const Double_t sinp = TMath::Sin(phi); + const Double_t cosp = TMath::Cos(phi); + + const Double_t v1 = sint*sinp; + const Double_t v2 = cosl - snlt*cosp; + + const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2); + + cos = -cosl*sinp / denom; + sin = (snlt-cslt*cosp) / denom; +} + +// -------------------------------------------------------------------------- +// +// RotationAngle +// +// calculates the angle for the rotation of the sky coordinate system +// with respect to the local coordinate system. This is identical +// to the rotation angle of the sky image in the camera. +// +// sinl [rad]: sine of observers latitude +// cosl [rad]: cosine of observers latitude +// theta [rad]: polar angle/zenith distance +// phi [rad]: rotation angle/azimuth +// +// Return angle [rad] in the range -pi, pi +// +// The convention is such, that the rotation angle is -pi/pi if +// right ascension and local rotation angle are counted in the +// same direction, 0 if counted in the opposite direction. +// +// (In other words: The rotation angle is 0 when the source culminates) +// +// Using vectors it can be done like: +// TVector3 v, p; +// v.SetMagThetaPhi(1, theta, phi); +// p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0); +// v = v.Cross(l)); +// v.RotateZ(-phi); +// v.Rotate(-theta) +// rho = TMath::ATan2(v(2), v(1)); +// +// For more information see TDAS 00-11, eqs. (18) and (20) +// +Double_t MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi) +{ + const Double_t snlt = sinl*TMath::Sin(theta); + const Double_t cslt = cosl*TMath::Cos(theta); + + const Double_t sinp = TMath::Sin(phi); + const Double_t cosp = TMath::Cos(phi); + + return TMath::ATan2(-cosl*sinp, snlt-cslt*cosp); +} Index: trunk/MagicSoft/Mars/mastro/MAstro.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstro.h (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstro.h (revision 3568) @@ -50,4 +50,8 @@ static Double_t UT2GMST(Double_t ut1); + // Rotation angle between local and sky coordinate system + static void RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos); + static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi); + ClassDef(MAstro, 0) }; Index: trunk/MagicSoft/Mars/mastro/MAstroCamera.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroCamera.cc (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstroCamera.cc (revision 3568) @@ -57,4 +57,5 @@ using namespace std; +// -------------------------------------------------------------------------- MAstroCamera::MAstroCamera() : fGeom(0), fMirrors(0) { @@ -63,4 +64,5 @@ } +// -------------------------------------------------------------------------- MAstroCamera::~MAstroCamera() { @@ -73,10 +75,11 @@ } -void MAstroCamera::SetMirrors(TClonesArray *arr) -{ - if (!arr || arr->GetClass()!=MGeomMirror::Class()) +// -------------------------------------------------------------------------- +void MAstroCamera::SetMirrors(TClonesArray &arr) +{ + if (arr.GetClass()!=MGeomMirror::Class()) return; - const Int_t n = arr->GetSize(); + const Int_t n = arr.GetSize(); if (!fMirrors) @@ -86,8 +89,9 @@ for (int i=0; iGetReflection(w, fGeom->GetCameraDist())*1000; - +// -------------------------------------------------------------------------- +Int_t MAstroCamera::ConvertToPad(const TVector3 &w, TVector2 &v) const +{ /* --- Use this to plot the 'mean grid' instead of the grid of a @@ -113,5 +116,5 @@ */ - + const TVector3 spot = fMirror0->GetReflection(w, fGeom->GetCameraDist())*1000; v.Set(spot(0), spot(1)); @@ -120,16 +123,5 @@ } -void MAstroCamera::DrawNet(const TRotation &trans) -{ - const TRotation rot(MAstroSky2Local(*fTime, *fObservatory)); - - TVector2 radec(fRaDec.Phi(), fRaDec.Theta()); - MAstroCatalog::DrawNet(radec, trans*rot, 2); - - const TVector3 zdaz0 = MAstroSky2Local(*fTime, *fObservatory)*fRaDec; - TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta()); - MAstroCatalog::DrawNet(zdaz, trans, 1); -} - +// -------------------------------------------------------------------------- TObject *FindObjectInPad(const char *name, TVirtualPad *pad) { @@ -155,4 +147,5 @@ } +// -------------------------------------------------------------------------- void MAstroCamera::AddPrimitives(Option_t *o) { @@ -171,13 +164,4 @@ const Bool_t hasdot = opt.Contains(".", TString::kIgnoreCase); const Bool_t usecam = opt.Contains("c", TString::kIgnoreCase); - - MAstroSky2Local rot(*fTime, *fObservatory); - - const Float_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta()); - - TString str = fTime->GetSqlDateTime(); - str += Form(" (\\alpha=%.1fh \\delta=%.1f\\circ) \\rho=%.1f\\circ", - fRaDec.Phi()/TMath::Pi()*12, 90-fRaDec.Theta()*TMath::RadToDeg(), - rho *TMath::RadToDeg()); // Get camera @@ -198,5 +182,5 @@ } - camera->SetTitle(str); + camera->SetTitle(GetPadTitle()); gPad->cd(1); @@ -221,20 +205,5 @@ } - TVector3 zdaz0 = fRaDec; - zdaz0 *= rot; - - cout << zdaz0.Phi()*TMath::RadToDeg() << " " << zdaz0.Theta()*TMath::RadToDeg() << endl; - - TVector3 test = zdaz0; - test *= rot.Inverse(); - - cout << test.Phi()*TMath::RadToDeg()/15 << " " << test.Theta()*TMath::RadToDeg() << endl; - - TRotation rot2; - rot2.RotateZ(-zdaz0.Phi()); - rot2.RotateY(-zdaz0.Theta()); - rot2.RotateZ(-TMath::Pi()/2); // align coordinate system - - DrawNet(rot2); + const TRotation rot(GetGrid(kTRUE)); MVector3 *radec; @@ -247,8 +216,6 @@ TVector3 star(*radec); - // Calculate local coordinates + // Rotate Star into telescope system star *= rot; - // Rotate Star into telescope system - star *= rot2; TVector3 mean; @@ -270,6 +237,4 @@ { TMarker *m=new TMarker(spot(0), spot(1), 1); - m->SetBit(kCannotPick); - m->SetBit(kCanDelete); m->SetMarkerColor(kMagenta); m->SetMarkerStyle(kDot); @@ -303,14 +268,10 @@ case kKey_Plus: fTime->SetMjd(fTime->GetMjd()+0.25/24); - SetBit(kHasChanged); - gPad->Modified(); - gPad->Update(); + Update(kTRUE); return; case kKey_Minus: fTime->SetMjd(fTime->GetMjd()-0.25/24); - SetBit(kHasChanged); - gPad->Modified(); - gPad->Update(); + Update(kTRUE); return; } Index: trunk/MagicSoft/Mars/mastro/MAstroCamera.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroCamera.h (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstroCamera.h (revision 3568) @@ -21,8 +21,7 @@ MGeomMirror *fMirror0; //! - Int_t ConvertToPad(const TVector3 &w, TVector2 &v); + Int_t ConvertToPad(const TVector3 &w, TVector2 &v) const; void AddPrimitives(Option_t *o); void SetRangePad() { } - void DrawNet(const TRotation &rot); void ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2); @@ -31,5 +30,5 @@ ~MAstroCamera(); - void SetMirrors(TClonesArray *arr); + void SetMirrors(TClonesArray &arr); void SetGeom(const MGeomCam &cam); Index: trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc (revision 3568) @@ -33,6 +33,8 @@ #include "MAstroCatalog.h" +#include #include #include +#include // INT_MAX (Suse 7.3/gcc 2.95) #include @@ -45,4 +47,5 @@ #include #include +#include #include @@ -59,5 +62,5 @@ using namespace std; - +/* class MRotation : public TRotation { @@ -76,5 +79,5 @@ } }; - +*/ /* MVector3 MVector3::GetZdAz(const MObservatory &obs, Double_t gmst) const @@ -134,5 +137,6 @@ } */ -MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(99), fToolTip(0), fObservatory(0), fTime(0) +// -------------------------------------------------------------------------- +MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(90), fToolTip(0), fObservatory(0), fTime(0) { fList.SetOwner(); @@ -140,4 +144,5 @@ } +// -------------------------------------------------------------------------- MAstroCatalog::~MAstroCatalog() { @@ -161,4 +166,5 @@ } +// -------------------------------------------------------------------------- TString MAstroCatalog::FindToken(TString &line, Char_t tok) { @@ -176,4 +182,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::atoi(const TSubString &sub) { @@ -181,4 +188,5 @@ } +// -------------------------------------------------------------------------- Float_t MAstroCatalog::atof(const TSubString &sub) { @@ -186,4 +194,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::atoi(const TString &s) { @@ -191,4 +200,5 @@ } +// -------------------------------------------------------------------------- Float_t MAstroCatalog::atof(const TString &s) { @@ -196,4 +206,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadXephem(TString catalog) { @@ -203,4 +214,10 @@ ifstream fin(catalog); + if (!fin) + { + gLog << err << "Cannot open file " << catalog << ": "; + gLog << strerror(errno) << endl; + return 0; + } Int_t add =0; @@ -276,4 +293,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadNGC2000(TString catalog) { @@ -283,4 +301,10 @@ ifstream fin(catalog); + if (!fin) + { + gLog << err << "Cannot open file " << catalog << ": "; + gLog << strerror(errno) << endl; + return 0; + } Int_t add=0; @@ -339,4 +363,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::ReadBSC(TString catalog) { @@ -346,4 +371,10 @@ ifstream fin(catalog); + if (!fin) + { + gLog << err << "Cannot open file " << catalog << ": "; + gLog << strerror(errno) << endl; + return 0; + } Int_t add=0; @@ -404,8 +435,8 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::Paint(Option_t *o) { -// if (!gPad->IsBatch()) -// gVirtualX->ClearWindow(); + SetRangePad(); if (TestBit(kHasChanged)) @@ -413,4 +444,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t transparent, const char *txt) { @@ -428,10 +460,9 @@ TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotLarge);; tip->SetMarkerColor(kBlack); - tip->SetBit(kCanDelete); - tip->SetBit(kCannotPick); AddMap(tip, new TString(str)); } -void MAstroCatalog::Update() +// -------------------------------------------------------------------------- +void MAstroCatalog::Update(Bool_t upd) { if (gPad && TestBit(kMustCleanup)) @@ -439,7 +470,14 @@ SetBit(kHasChanged); gPad->Modified(); - } -} - + if (upd) + gPad->Update(); + } +} + +// -------------------------------------------------------------------------- +// +// Set the observation time. Necessary to use local coordinate +// system. +// void MAstroCatalog::SetTime(const MTime &time) { @@ -449,4 +487,9 @@ } +// -------------------------------------------------------------------------- +// +// Set the observatory location. Necessary to use local coordinate +// system. +// void MAstroCatalog::SetObservatory(const MObservatory &obs) { @@ -456,26 +499,34 @@ } -Int_t MAstroCatalog::ConvertToPad(const TVector3 &w0, TVector2 &v) +// -------------------------------------------------------------------------- +// +// Convert the vector to pad coordinates. After conversion +// the x- coordinate of the vector must be the x coordinate +// of the pad - the same for y. If the coordinate is inside +// the current draw area return kTRUE, otherwise kFALSE. +// If it is an invalid coordinate return kERROR +// +Int_t MAstroCatalog::ConvertToPad(const TVector3 &w0, TVector2 &v) const { TVector3 w(w0); - // Stretch such, that the X-component is alwas the same. Now - // Y and Z contains the crossing point between the star-light - // and the plain of a virtual screen (ccd...) + // Stretch such, that the Z-component is alwas the same. Now + // X and Y contains the intersection point between the star-light + // and the plain of a virtual plain screen (ccd...) if (TestBit(kPlainScreen)) - w *= TMath::RadToDeg()/w.X(); - else - w.SetMag(TMath::RadToDeg()); - - v.Set(w(1), w(2)); - - if (w(0)<0) + w *= 1./w(2); + + w *= TMath::RadToDeg(); + v.Set(w(0), w(1)); + + if (w(2)<0) return kERROR; - return w(1)>gPad->GetX1() && w(2)>gPad->GetY1() && - w(1)GetX2() && w(2)GetY2(); -} - -Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v) + return w(0)>gPad->GetX1() && w(1)>gPad->GetY1() && + w(0)GetX2() && w(1)GetY2(); +} + +// -------------------------------------------------------------------------- +Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v) const { MVector3 w; @@ -486,4 +537,5 @@ } +// -------------------------------------------------------------------------- Bool_t MAstroCatalog::DrawLine(const TVector2 &v, Double_t dx, Double_t dy, const TRotation &rot, Int_t type) { @@ -501,18 +553,14 @@ TLine *line = new TLine(v0.X(), v0.Y(), v1.X(), v1.Y()); - line->SetBit(kCanDelete); line->SetLineStyle(kDashDotted); //kDashed, kDotted, kDashDotted line->SetLineColor(kWhite+type*2); - line->SetBit(kCannotPick); AddMap(line); const TVector2 deg = v*TMath::RadToDeg(); TString txt = type==1 ? - Form("Ra=%.1fh Dec=%.1fd", fmod(deg.X()/15+48, 24), fmod(90-deg.Y()+270,180)-90) : + Form("Ra=%.2fh Dec=%.1fd", fmod(deg.X()/15+48, 24), fmod(90-deg.Y()+270,180)-90) : Form("Zd=%.1fd Az=%.1fd", fmod(deg.Y()+270,180)-90, fmod(deg.X()+720, 360)); TMarker *tip=new TMarker(v0.X(), v0.Y(), kDot); - tip->SetBit(kCanDelete); - tip->SetBit(kCannotPick); tip->SetMarkerColor(kWhite+type*2); AddMap(tip, new TString(txt)); @@ -521,5 +569,9 @@ } - +// -------------------------------------------------------------------------- +// +// Use "local" draw option to align the display to the local +// coordinate system instead of the sky coordinate system. +// void MAstroCatalog::Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type) { @@ -570,5 +622,6 @@ } -void MAstroCatalog::DrawNet(const TVector2 &v0, const TRotation &rot, Int_t type) +// -------------------------------------------------------------------------- +void MAstroCatalog::DrawGrid(const TVector3 &v0, const TRotation &rot, Int_t type) { TArrayI dx(1); @@ -576,23 +629,26 @@ // align to 1deg boundary - TVector2 v = v0*TMath::RadToDeg(); + TVector2 v(v0.Phi()*TMath::RadToDeg(), v0.Theta()*TMath::RadToDeg()); v.Set((Float_t)TMath::Nint(v.X()), (Float_t)TMath::Nint(v.Y())); // calculate stepsizes based on visible FOV - Int_t stepx=1; - - if (fabs(90-v.Y())>90-fRadiusFOV || fabs(90-v.Y())180-fRadiusFOV) + stepx=36; else { // This is a rough estimate how many degrees are visible - const Float_t m = log(fRadiusFOV/180.)/log(90./fRadiusFOV-1); + const Float_t m = log(fRadiusFOV/180.)/log(90./(fRadiusFOV+1)+1); const Float_t t = log(180.)-m*log(fRadiusFOV); - const Int_t n = (Int_t)(exp(m*log(90-fabs(90-v.Y()))+t)+0.5); - stepx = n<6 ? 1 : n/6; - } - - const Int_t n = (Int_t)(fRadiusFOV+1); - Int_t stepy = n<4 ? 1 : n/4; + const Float_t f = m*log(90-fabs(90-v.Y()))+t; + const Int_t nx = (Int_t)(exp(f)+0.5); + stepx = nx<4 ? 1 : nx/4; + if (stepx>36) + stepx=36; + } + + const Int_t ny = (Int_t)(fRadiusFOV+1); + Int_t stepy = ny<4 ? 1 : ny/4; // align stepsizes to be devisor or 180 and 90 @@ -602,10 +658,18 @@ stepy++; - // align to step-size boundary + // align to step-size boundary (search for the nearest one) + Int_t dv = 1; while ((int)(v.X())%stepx) - v.Set(v.X()+1, v.Y()); - + { + v.Set(v.X()+dv, v.Y()); + dv = -TMath::Sign(TMath::Abs(dv)+1, dv); + } + + dv = 1; while ((int)(v.Y())%stepy) - v.Set(v.X(), v.Y()+1); + { + v.Set(v.X(), v.Y()+dv); + dv = -TMath::Sign(TMath::Abs(dv)+1, dv); + } // draw... @@ -615,21 +679,99 @@ } +// -------------------------------------------------------------------------- +// +// Get a rotation matrix which aligns the pointing position +// to the center of the x,y plain +// +TRotation MAstroCatalog::AlignCoordinates(const TVector3 &v) const +{ + TRotation trans; + trans.RotateZ(-v.Phi()); + trans.RotateY(-v.Theta()); + trans.RotateZ(-TMath::Pi()/2); + return trans; +} + +// -------------------------------------------------------------------------- +// +// Return the rotation matrix which converts either sky or +// local coordinates to coordinates which pole is the current +// pointing direction. +// +TRotation MAstroCatalog::GetGrid(Bool_t local) +{ + const Bool_t enable = fTime && fObservatory; + + if (!local) + { + const TRotation trans(AlignCoordinates(fRaDec)); + + DrawGrid(fRaDec, trans, 1); + + if (enable) + { + const MAstroSky2Local rot(*fTime, *fObservatory); + DrawGrid(rot*fRaDec, trans*rot.Inverse(), 2); + } + + return trans; + } + + if (local && enable) + { + const MAstroSky2Local rot(*fTime, *fObservatory); + + const TRotation trans(AlignCoordinates(rot*fRaDec)); + + DrawGrid(fRaDec, trans*rot, 1); + DrawGrid(rot*fRaDec, trans, 2); + + return trans*rot; + } + + return TRotation(); +} + +// -------------------------------------------------------------------------- +TString MAstroCatalog::GetPadTitle() const +{ + const Double_t ra = fRaDec.Phi()*TMath::RadToDeg(); + const Double_t dec = (TMath::Pi()/2-fRaDec.Theta())*TMath::RadToDeg(); + + TString txt; + txt += Form("\\alpha=%.2fh ", fmod(ra/15+48, 24)); + txt += Form("\\delta=%.1f\\circ ", fmod(dec+270,180)-90); + txt += Form("/ FOV=%.1f\\circ", fRadiusFOV); + + if (!fTime || !fObservatory) + return txt; + + const MAstroSky2Local rot(*fTime, *fObservatory); + const TVector3 loc = rot*fRaDec; + + const Double_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta()); + + const Double_t zd = TMath::RadToDeg()*loc.Theta(); + const Double_t az = TMath::RadToDeg()*loc.Phi(); + + txt.Prepend("#splitline{"); + + txt += Form(" \\theta=%.1fh ", fmod(zd+270,180)-90); + txt += Form("\\phi=%.1f\\circ ", fmod(az+720, 360)); + txt += Form(" / \\rho=%.1f\\circ", rho*TMath::RadToDeg()); + txt += "}{<"; + txt += fTime->GetSqlDateTime(); + txt += ">}"; + return txt; +} + +// -------------------------------------------------------------------------- void MAstroCatalog::AddPrimitives(Option_t *o) { - // Precalc Sin/Cos... - TRotation trans; - trans.RotateZ(-fRaDec.Phi()); - trans.Rotate(TMath::Pi()/2-fRaDec.Theta(), TVector3(0, 1, 0)); - - if (fTime && fObservatory) - { - const TRotation rot(MAstroSky2Local(*fTime, *fObservatory)); - const TVector3 zdaz0 = rot*fRaDec; - const TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta()); - DrawNet(zdaz, trans*rot.Inverse(), 2); - } - - const TVector2 radec(fRaDec.Phi(), fRaDec.Theta()); - DrawNet(radec, trans, 1); + const Bool_t local = TString(o).Contains("local", TString::kIgnoreCase); + + cout << "Opt: " << o << endl; + + const TRotation rot(GetGrid(local)); TIter Next(&fList); @@ -639,9 +781,19 @@ // FIXME: Check Magnitude! TVector2 s(v->Phi(), v->Theta()); - if (Convert(trans, s)==kTRUE) + if (Convert(rot, s)==kTRUE) DrawStar(s.X(), TMath::Pi()/2-s.Y(), *v, kFALSE); } -} - + + TPaveText *pv = new TPaveText(0.01, 0.90, 0.63, 0.99, "brNDC"); + pv->AddText(GetPadTitle()); + AddMap(pv); + + TMarker *mk=new TMarker(0, 0, kMultiply); + mk->SetMarkerColor(kBlack); + mk->SetMarkerSize(1.5); + AddMap(mk); +} + +// -------------------------------------------------------------------------- void MAstroCatalog::SetRangePad() { @@ -658,4 +810,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::DrawPrimitives(Option_t *o) { @@ -680,4 +833,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::DrawMap() { @@ -688,4 +842,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::Draw(Option_t *o) { @@ -770,4 +925,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::ExecuteEventKbd(Int_t keycode, Int_t keysym) { @@ -821,4 +977,5 @@ } +// -------------------------------------------------------------------------- Int_t MAstroCatalog::DistancetoPrimitive(Int_t px, Int_t py) { @@ -843,4 +1000,5 @@ } +// -------------------------------------------------------------------------- void MAstroCatalog::ShowToolTip(Int_t px, Int_t py, const char *txt) { @@ -858,3 +1016,2 @@ fToolTip->Show(x+4, y+4); } - Index: trunk/MagicSoft/Mars/mastro/MAstroCatalog.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroCatalog.h (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstroCatalog.h (revision 3568) @@ -74,8 +74,9 @@ { private: - Double_t fLimMag; // [1] Limiting Magnitude - Double_t fRadiusFOV; // [deg] Radius of Field of View - - TGToolTip *fToolTip; //! + Double_t fLimMag; // [1] Limiting Magnitude + Double_t fRadiusFOV; // [deg] Radius of Field of View + + TExMap fMapG; //! A map with all gui primitives and tooltips + TGToolTip *fToolTip; //! The tooltip currently displayed void ShowToolTip(Int_t px, Int_t py, const char *txt); @@ -92,37 +93,18 @@ //#endif -protected: - enum { - kHasChanged = BIT(15), - kGuiActive = BIT(16), - kPlainScreen = BIT(17) - }; - - TExMap fMapG; - TList fList; // List of stars loaded - MVector3 fRaDec; // pointing position - - MObservatory *fObservatory; // Possible obervatora location - MTime *fTime; // Possible observation time - - virtual Int_t ConvertToPad(const TVector3 &w, TVector2 &v); - virtual Int_t Convert(const TRotation &rot, TVector2 &v); - virtual Bool_t DrawLine(const TVector2 &v0, Double_t dx, Double_t dy, const TRotation &rot, Int_t type); - virtual void AddPrimitives(Option_t *o); - virtual void DrawPrimitives(Option_t *o); - virtual void SetRangePad(); - void Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type); - void DrawNet(const TVector2 &v0, const TRotation &rot, Int_t type); - void DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t t, const char *txt=0); - void Paint(Option_t *o=""); - Int_t DistancetoPrimitive(Int_t px, Int_t py); - void Update(); - - void ExecuteEventKbd(Int_t keycode, Int_t keysym); - void ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2); - - void DrawMap(); - void AddMap(void *k, void *v=0) { fMapG.Add(fMapG.GetSize(), (Long_t)k, (Long_t)v); } - void DeleteMap() + virtual Int_t ConvertToPad(const TVector3 &w, TVector2 &v) const; + virtual void AddPrimitives(Option_t *o); + virtual void SetRangePad(); + + Int_t Convert(const TRotation &rot, TVector2 &v) const; + void Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type); + void DrawPrimitives(Option_t *o); + Bool_t DrawLine(const TVector2 &v0, Double_t dx, Double_t dy, const TRotation &rot, Int_t type); + void DrawGrid(const TVector3 &v0, const TRotation &rot, Int_t type); + TRotation AlignCoordinates(const TVector3 &v) const; + void Paint(Option_t *o=""); + Int_t DistancetoPrimitive(Int_t px, Int_t py); + void DrawMap(); + void DeleteMap() { Long_t key, val; @@ -148,4 +130,32 @@ } +protected: + enum { + kHasChanged = BIT(15), + kGuiActive = BIT(16), + kPlainScreen = BIT(17) + }; + + TList fList; // List of stars loaded + MVector3 fRaDec; // pointing position + + MObservatory *fObservatory; // Possible obervatora location + MTime *fTime; // Possible observation time + + virtual TString GetPadTitle() const; + TRotation GetGrid(Bool_t local); + void DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t t, const char *txt=0); + void Update(Bool_t upd=kFALSE); + + void ExecuteEventKbd(Int_t keycode, Int_t keysym); + void ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2); + + void AddMap(TObject *k, void *v=0) + { + k->SetBit(kCanDelete); + k->SetBit(kCannotPick); + fMapG.Add(fMapG.GetSize(), (Long_t)k, (Long_t)v); + } + public: MAstroCatalog(); @@ -161,5 +171,5 @@ if (deg>max) deg=max; - if (deg<=0) + if (deg<1) deg=1; @@ -189,4 +199,9 @@ void Draw(Option_t *o=""); + void SetDrawOption(Option_t *option="") + { + TObject::SetDrawOption(option); + Update(kTRUE); + } //*MENU* virtual void EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected=0); Index: trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc (revision 3568) @@ -62,4 +62,5 @@ #include "MAstroSky2Local.h" +#include "MAstro.h" #include "MTime.h" #include "MObservatory.h" @@ -128,18 +129,12 @@ // seen with an Alt/Az telescope. // +// For more information see MAstro::RotationAngle +// Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec) const { - TVector3 loc; - loc.SetMagThetaPhi(1, TMath::Pi()/2-dec, ra); - loc *= *this; + TVector3 v; + v.SetMagThetaPhi(1, TMath::Pi()/2-dec, ra); + v *= *this; - TRotation rot; - rot.RotateZ(-loc.Phi()); - rot.RotateY(-loc.Theta()); - - TVector3 v(1, 0, 0); - v *= *this; - v *= rot; - - return TMath::ATan2(v(1), v(0)); + return MAstro::RotationAngle(ZZ(), XZ(), v.Theta(), v.Phi()); } Index: trunk/MagicSoft/Mars/mastro/MObservatory.cc =================================================================== --- trunk/MagicSoft/Mars/mastro/MObservatory.cc (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MObservatory.cc (revision 3568) @@ -101,97 +101,31 @@ { *fLog << all; - *fLog << fObservatoryName << endl; - *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude*kRad2Deg) : -(fLatitude*kRad2Deg)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl; - *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude*kRad2Deg) : -(fLongitude*kRad2Deg)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl; - *fLog << "Height " << fHeight << "m" << endl; + *fLog << underline << fObservatoryName << ":" << endl; + *fLog << "Latitude: " << TMath::Abs(fLatitude*kRad2Deg) << " deg " << (fLatitude > 0 ? "W" : "E") << endl; + *fLog << "Longitude: " << TMath::Abs(fLongitude*kRad2Deg) << " deg " << (fLongitude < 0 ? "N" : "S") << endl; + *fLog << "Height: " << fHeight << "m" << endl; } // -------------------------------------------------------------------------- // -// RotationAngle -// -// calculates the angle for the rotation of the sky image in the camera; -// this angle is a function of the local coordinates +// Get the corresponding rotation angle of the sky coordinate system +// seen with an Alt/Az telescope. // -// theta [rad]: polar angle/zenith distance -// phi [rad]: rotation angle/azimuth -// -// Return sin/cos component of angle -// -// calculate rotation angle alpha of sky image in camera -// (see TDAS 00-11, eqs. (18) and (20)) +// For more information see MAstro::RotationAngle // void MObservatory::RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const { - const Double_t sint = TMath::Sin(theta); - const Double_t cost = TMath::Cos(theta); - - const Double_t sinl = fSinLatitude*sint; - const Double_t cosl = fCosLatitude*cost; - - const Double_t sinp = TMath::Sin(phi); - const Double_t cosp = TMath::Cos(phi); - - const Double_t v1 = sint*sinp; - const Double_t v2 = cosl - sinl*cosp; - - const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2); - - sin = fCosLatitude*sinp / denom; - cos = (sinl - cosl*cosp) / denom; + MAstro::RotationAngle(fSinLatitude, fCosLatitude, theta, phi, sin, cos); } // -------------------------------------------------------------------------- // -// RotationAngle -// -// calculates the angle for the rotation of the sky image in the camera; -// this angle is a function of the local coordinates +// Get the corresponding rotation angle of the sky coordinate system +// seen with an Alt/Az telescope. // -// theta [rad]: polar angle/zenith distance -// phi [rad]: rotation angle/azimuth -// -// Return RotationAngle in rad -// -// calculate rotation angle alpha of sky image in camera -// (see TDAS 00-11, eqs. (18) and (20)) +// For more information see MAstro::RotationAngle // Double_t MObservatory::RotationAngle(Double_t theta, Double_t phi) const { - const Double_t sint = TMath::Sin(theta); - const Double_t cost = TMath::Cos(theta); - - const Double_t sinp = TMath::Sin(phi); - const Double_t cosp = TMath::Cos(phi); - - const Double_t v1 = sint*sinp; - const Double_t v2 = fCosLatitude*cost - fSinLatitude*sint*cosp; - - const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2); - - return TMath::ASin((fCosLatitude*sinp) / denom); + return MAstro::RotationAngle(fSinLatitude, fCosLatitude, theta, phi); } - -// -------------------------------------------------------------------------- -// -// RotationAngle -// -// calculates the angle for the rotation of the sky image in the camera; -// this angle is a function of the sky coordinates, the observatory -// location and the time -// -// ra [rad]: Right ascension -// dec [rad]: Declination -// -// Return RotationAngle in rad -// -Double_t MObservatory::RotationAngle(Double_t ra, Double_t dec, const MTime &t) const -{ - const Double_t alpha = t.GetGmst() + GetElong(); - - TVector3 v; - v.SetMagThetaPhi(1, TMath::Pi()/2-dec, alpha-ra); - v.RotateY(GetPhi()-TMath::Pi()/2); - - return RotationAngle(v.Theta(), v.Phi()); -} Index: trunk/MagicSoft/Mars/mastro/MObservatory.h =================================================================== --- trunk/MagicSoft/Mars/mastro/MObservatory.h (revision 3567) +++ trunk/MagicSoft/Mars/mastro/MObservatory.h (revision 3568) @@ -58,5 +58,4 @@ void RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const; Double_t RotationAngle(Double_t theta, Double_t phi) const; - Double_t RotationAngle(Double_t ra, Double_t dec, const MTime &t) const; LocationName_t GetObservatoryKey() const { return fObservatoryKey; } Index: trunk/MagicSoft/Mars/mbase/MEvtLoop.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MEvtLoop.cc (revision 3567) +++ trunk/MagicSoft/Mars/mbase/MEvtLoop.cc (revision 3568) @@ -496,14 +496,14 @@ *fLog << all << "Ready!" << endl << endl; - *fLog << dec << endl << "CPU - " - << "Time: " << clock.CpuTime() << "s" - << " for " << numcnts << " Events" - << " --> " << numcnts/clock.CpuTime() << " Events/s" - << endl; - *fLog << "Real - " - << "Time: " << clock.RealTime() << "s" - << " for " << numcnts << " Events" - << " --> " << numcnts/clock.RealTime() << " Events/s" - << endl << endl; + *fLog << dec << endl << "CPU - Time: "; + *fLog << clock.CpuTime() << "s" << " for " << numcnts << " Events"; + if (numcnts>0) + *fLog << " --> " << numcnts/clock.CpuTime() << " Events/s"; + *fLog << endl << "Real - Time: "; + *fLog << clock.RealTime() << "s" << " for " << numcnts << " Events"; + if (numcnts>0) + *fLog << " --> " << numcnts/clock.RealTime() << " Events/s"; + + *fLog << endl << endl; return rc!=kERROR; Index: trunk/MagicSoft/Mars/mbase/MParList.cc =================================================================== --- trunk/MagicSoft/Mars/mbase/MParList.cc (revision 3567) +++ trunk/MagicSoft/Mars/mbase/MParList.cc (revision 3568) @@ -329,4 +329,6 @@ // 'name' is the name of the object you are searching for. // +// In words: Find object name and check whether it inherits from classname +// TObject *MParList::FindObject(const char *name, const char *classname) const { @@ -427,4 +429,7 @@ // "Name;1", "Name;2", ... If a semicolon is detected leading dots // are stripped from the object-name (eg. "name;5.") +// +// In words: Create object of type classname and set its name to objname. +// If an object with objname already exists return it. // MParContainer *MParList::FindCreateObj(const char *classname, const char *objname) Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc (revision 3567) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc (revision 3568) @@ -38,7 +38,7 @@ // // ///////////////////////////////////////////////////////////////////////////// - #include "MCT1ReadAscii.h" +#include #include @@ -132,12 +132,15 @@ const char *expname = gSystem->ExpandPathName(name); fIn = new ifstream(expname); + + const Bool_t noexist = !(*fIn); + if (noexist) + { + *fLog << err << "Cannot open file " << expname << ": "; + *fLog << strerror(errno) << endl; + } + else + *fLog << inf << "Open file: '" << name << "'" << endl; + delete [] expname; - - const Bool_t noexist = !(*fIn); - - if (noexist) - *fLog << dbginf << "Cannot open file '" << name << "'" << endl; - else - *fLog << "Open file: '" << name << "'" << endl; // Index: trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 3567) +++ trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc (revision 3568) @@ -47,8 +47,10 @@ #include "MCT1ReadPreProc.h" +#include #include #include #include +#include #define LINUX @@ -81,5 +83,4 @@ #include "MBinning.h" -#include "TRandom3.h" #include "MParameters.h" @@ -590,4 +591,10 @@ fIn = new ifstream(fname); + if (!*fIn) + { + *fLog << err << "Cannot open file " << fname << ": "; + *fLog << strerror(errno) << endl; + return kFALSE; + } *fLog << inf << "-----------------------------------------------------------------------" << endl; Index: trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc (revision 3567) +++ trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc (revision 3568) @@ -32,4 +32,5 @@ #include "MReadRflFile.h" +#include #include @@ -315,5 +316,6 @@ if (!*fIn) { - cout << "Error openeng file " << name << "." << endl; + *fLog << err << "Cannot open file " << name << ": "; + *fLog << strerror(errno) << endl; return kFALSE; } Index: trunk/MagicSoft/Mars/mfileio/MReadTree.cc =================================================================== --- trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 3567) +++ trunk/MagicSoft/Mars/mfileio/MReadTree.cc (revision 3568) @@ -549,5 +549,5 @@ if (TestBit(kChainWasChanged)) { - *fLog << inf << "Scanning chain... " << flush; + *fLog << inf << "Scanning chain " << fChain->GetName() << "... " << flush; fNumEntries = (UInt_t)fChain->GetEntries(); *fLog << fNumEntries << " events found." << endl; Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.cc =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3567) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.cc (revision 3568) @@ -95,4 +95,8 @@ // - a more clever (and faster) algorithm to fill the histogram, eg by // calculating alpha once and fill the two coils around the mean +// - more drawing options... +// - Move Significance() to a more 'general' place and implement +// also different algorithms like (Li/Ma) +// - implement fit for best alpha distribution -- online (Paint) // ////////////////////////////////////////////////////////////////////////////// @@ -101,7 +105,9 @@ #include #include +#include #include #include #include +#include #include "MGeomCam.h" @@ -236,4 +242,7 @@ // Calculate Significance as // significance = (s-b)/sqrt(s+k*k*b) mit k=s/b +// +// s: total number of events in signal region +// b: number of background events in signal region // Double_t MHFalseSource::Significance(Double_t s, Double_t b) @@ -245,4 +254,37 @@ } +// -------------------------------------------------------------------------- +// +// calculates the significance according to Li & Ma +// ApJ 272 (1983) 317 +// +// s: total number of events in signal region +// +// Double_t fGamma; // Nbg = Non - gamma * Noff +// - the effective number of background events (fNoff), and fGamma : +// +// fGamma = b / fNoff; +// fGamma = fdNbg / sqrt(fNoff); +// fGamma = fdNbg*fdNbg / fNbg; +// fNoff = b*b / (fdNbg*fdNbg); +// Double_t fNbg; // number of background events in signal region +// b = fNOff *fGamma; +/* +Double_t MHFindSignificance::SignificanceLiMa(Double_t non, Double_t noff, Double_t gamma) +{ + if (gamma <= 0.0 || non <= 0.0 || noff <= 0.0) + { + *siglima = 0.0; + return kFALSE; + } + + Double_t help1 = TMath::Log( (gamma+1)*s / (gamma*(s+noff)) ); + Double_t help2 = TMath::Log( (gamma+1)*noff / ( s+noff ) ); + + Double_t siglima = TMath::Sqrt((s*help1+noff*help2)*2); + + return nonGetRa(); - const Double_t dec = fPointPos->GetDec(); - - rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec); - } + rho = fPointPos->RotationAngle(*fObservatory, *fTime); + //if (fPointPos) + // rho = fPointPos->RotationAngle(*fObservatory); MSrcPosCam src; @@ -554,5 +593,4 @@ h4->Draw("colz"); h4->SetBit(kCanDelete); - } @@ -587,5 +625,5 @@ gPad->cd(); } - +/* Double_t fcn(Double_t *arg, Double_t *p) { @@ -595,5 +633,5 @@ const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx); - const Double_t f2 = p[3]*x*x + p[4]; + const Double_t f2 = p[3] + p[5]*x*x; return f1 + f2; @@ -602,5 +640,5 @@ Double_t FcnI1(Double_t x, Double_t *p) { - return (p[3]*x*x/3+p[4])*x; + return (p[5]*x*x/3+p[3])*x; } Double_t FcnI2(Double_t x, Double_t *p) @@ -615,54 +653,212 @@ return f2; } - - -void MHFalseSource::FitSignificance(Float_t sigmax, Float_t bgmin, Float_t bgmax) -{ - TH1D h0a("A", "Parameter A", 50, 0, 4000); - TH1D h0b("a", "Parameter a", 50, 0, 4000); - TH1D h1("mu", "Parameter \\mu", 50, -1, 1); - TH1D h2("sigma", "Parameter \\sigma", 50, 0, 90); - TH1D h3("b", "Parameter b", 50, -0.1, 0.1); - TH1D h4a("chisq1", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35); - TH1D h5a("prob1", "Fit probability", 50, 0, 1.1); - TH1D h4b("chisq2", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35); - TH1D h5b("prob", "Fit probability", 50, 0, 1.1); - TH1D h6("Signif", "Significance", 50, -20, 20); - h0a.SetDirectory(NULL); - h0b.SetDirectory(NULL); - h1.SetDirectory(NULL); - h2.SetDirectory(NULL); - h3.SetDirectory(NULL); - h4a.SetDirectory(NULL); - h5b.SetDirectory(NULL); - h4a.SetDirectory(NULL); - h5b.SetDirectory(NULL); - h6.SetDirectory(NULL); +*/ +/* +class MHSignificance : public MH +{ +private: + TH1D fHist; + + MParameterD *fParam; + +public: + MHSignificance() : fParam(0) + { + fHist.SetName("Alpha"); + fHist.SetTitle("Distribution of \\alpha"); + fHist.SetXTitle("\\alpha [\\circ]"); + fHist.SetYTitle("Counts"); + } + Int_t SetupFill(MParList *p) + { + fHist.Reset(); + + fParam = (MParameterD*)p->FindCreateObj("Significance", "MParameterD"); + if (fParam) + return kFALSE; + + return kTRUE; + } + Int_t Process(MParContainer *p, Double_t w=1) + { + MHillasSrc *hil = dynamic_cast(p); + if (!hil) + { + *fLog << err << dbginf << "Got no MHillasSrc as argument of Fill()..." << endl; + return kFALSE; + } + + fHist->Fill(hil->GetAlpha(), w); + + return kTRUE; + } + Int_t Finalize() + { + if (fHist.GetEntries()==0) + { + *fLog << err << "Histogram empty." << endl; + return kFALSE; + } + + Float_t sigmax=15; + Float_t bgmin =45; + Float_t bgmax =80; + + fHist.SetNameTitle("Significance", + Form("Fit Region: Signal<%.1f\\circ, %.1f\\circFit(&func, "N0Q", "", bgmin, bgmax); + + // Now fit a gaus in the on region on top of the polynom + func.SetParameter(0, alpha0-func.GetParameter(4)); + func.SetParameter(2, sigmax*0.75); + + func.ReleaseParameter(0); + func.ReleaseParameter(2); + func.FixParameter(3, func.GetParameter(3)); + func.FixParameter(5, func.GetParameter(5)); + + func.SetParLimits(2, 0, 80); + h->Fit(&func, "N0Q", "", 0, sigmax); + + TArrayD p(func.GetNpar(), func.GetParameters()); + + Double_t sig=0; + + const Int_t n = hist->GetBin(ix+1, iy+1); + if (!(func.GetParameter(0)>alpha0*2 || + func.GetParameter(2)<2.5 || + func.GetParameter(2)>70)) + { + // Implementing the integral as analytical function + // gives the same result in the order of 10e-5 + // and it is not faster at all... + const Double_t s = func.Integral(0, 15); + + func.SetParameter(0, 0); + func.SetParameter(2, 1); + + const Double_t b = func.Integral(0, 15); + + sig = Significance(s, b); + } + + fParam->SetValue(sig); + fParam->SetReadyToSave(); + return kTRUE; + } + ClassDef(MHSignificance, 0) +}; +*/ + + +// -------------------------------------------------------------------------- +// +// This is a preliminary implementation of a alpha-fit procedure for +// all possible source positions. It will be moved into its own +// more powerfull class soon. +// +// The fit function is "gaus(0)+pol2(3)" which is equivalent to: +// [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2 +// or +// A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2 +// +// Parameter [1] is fixed to 0 while the alpha peak should be +// symmetric around alpha=0. +// +// Parameter [4] is fixed to 0 because the first derivative at +// alpha=0 should be 0, too. +// +// In a first step the background is fitted between bgmin and bgmax, +// while the parameters [0]=0 and [2]=1 are fixed. +// +// In a second step the signal region (alphaSetDirectory(0); + TH1 *hists = fHist.Project3D("xy_fit"); + hists->SetDirectory(0); + TH1 *histb = fHist.Project3D("xy_fit"); + histb->SetDirectory(0); hist->Reset(); + hists->Reset(); + histb->Reset(); hist->SetNameTitle("Significance", Form("Fit Region: Signal<%.1f\\circ, %.1f\\circSetNameTitle("Excess", Form("Number of excess events for \\alpha<%.0f\\circ", sigint)); + histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint)); hist->SetXTitle(fHist.GetXaxis()->GetTitle()); + hists->SetXTitle(fHist.GetXaxis()->GetTitle()); + histb->SetXTitle(fHist.GetXaxis()->GetTitle()); hist->SetYTitle(fHist.GetXaxis()->GetTitle()); + hists->SetYTitle(fHist.GetXaxis()->GetTitle()); + histb->SetYTitle(fHist.GetXaxis()->GetTitle()); // xmin, xmax, npar - TF1 func("MyFunc", fcn, 0, 90, 5); - TArrayD par(5); - - func.SetParName(0, "A"); - func.SetParName(1, "mu"); - func.SetParName(2, "sigma"); - func.SetParName(3, "a"); - func.SetParName(4, "b"); - + //TF1 func("MyFunc", fcn, 0, 90, 6); + // Implementing the function yourself is only about 5% faster + TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90); + TArrayD maxpar(func.GetNpar()); + + /* + func.SetParName(0, "A"); + func.SetParName(1, "mu"); + func.SetParName(2, "sigma"); + func.SetParName(3, "a"); + func.SetParName(4, "b"); + func.SetParName(5, "c"); + */ + + func.FixParameter(1, 0); + func.FixParameter(4, 0); + func.SetParLimits(2, 0, 90); func.SetParLimits(3, -1, 1); @@ -676,11 +872,22 @@ Int_t maxy=0; + TStopwatch clk; + clk.Start(); + + *fLog << inf; + *fLog << "Signal fit: alpha < " << sigmax << endl; + *fLog << "Integration: alpha < " << sigint << endl; + *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl; + *fLog << "Polynom order: " << (int)polynom << endl; + *fLog << "Fitting False Source Plot..." << flush; + TH1 *h=0; - for (int ix=0; ixGetBinContent(1); + const Double_t alphaw = h->GetXaxis()->GetBinWidth(1); // Check for the regios which is not filled... @@ -693,8 +900,9 @@ // First fit a polynom in the off region func.FixParameter(0, 0); - func.FixParameter(1, 0); func.FixParameter(2, 1); func.ReleaseParameter(3); - func.ReleaseParameter(4); + + for (int i=5; iFit(&func, "N0Q", "", bgmin, bgmax); @@ -702,9 +910,5 @@ h4a.Fill(func.GetChisquare()); - h5a.Fill(func.GetProb()); - - // Now fit a gaus in the on region on top of the polynom - func.SetParameter(0, alpha0-func.GetParameter(4)); - func.SetParameter(2, sigmax*0.75); + //h5a.Fill(func.GetProb()); //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1)); @@ -715,34 +919,57 @@ func.ReleaseParameter(2); func.FixParameter(3, func.GetParameter(3)); - func.FixParameter(4, func.GetParameter(4)); - - func.SetParLimits(2, 0, 80); + for (int i=5; iFit(&func, "N0Q", "", 0, sigmax); //*fLog << dbg << " " << func.GetParameter(0) << " " << func.GetParameter(1) << " " << func.GetParameter(2) << endl; + TArrayD p(func.GetNpar(), func.GetParameters()); + // Fill results into some histograms - h0a.Fill(func.GetParameter(0)); - h0b.Fill(func.GetParameter(4)); - h1.Fill(func.GetParameter(1)); - h2.Fill(func.GetParameter(2)); - h3.Fill(func.GetParameter(3)); + h0a.Fill(p[0]); + h0b.Fill(p[3]); + h1.Fill(p[1]); + h2.Fill(p[2]); + if (polynom>1) + h3.Fill(p[5]); h4b.Fill(func.GetChisquare()); - h5b.Fill(func.GetProb()); - - Double_t sig=0; - - const Int_t n = hist->GetBin(ix+1, iy+1); - if (!(func.GetParameter(0)>alpha0*2 || - func.GetParameter(2)<2.5 || - func.GetParameter(2)>70 - /*func.GetProb()<0.005 ||*/)) + //h5b.Fill(func.GetProb()); + + // Implementing the integral as analytical function + // gives the same result in the order of 10e-5 + // and it is not faster at all... + //const Bool_t ok = /*p[0]>=0 && /*p[0]1.75;// && p[2]<88.5; + /* + if (p[0]<-fabs(alpha0-p[3])*7 && p[2]SetBinContent(n, sig==0?1e-15:sig); + */ + + const Double_t s = func.Integral(0, sigint); + + func.SetParameter(0, 0); + func.SetParameter(2, 1); + + const Double_t b = func.Integral(0, sigint); + const Double_t sig = Significance(s, b); + + const Int_t n = hist->GetBin(ix, iy); + hists->SetBinContent(n, s-b); + histb->SetBinContent(n, b); + + hist->SetBinContent(n, sig); if (sig!=0) h6.Fill(sig); @@ -751,14 +978,21 @@ { maxs = sig; - maxx = ix+1; - maxy = iy+1; - for(int i=0; iSetRangeUser(0, maxalpha0*1.5); h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5); + //hists->SetMinimum(GetMinimumGT(*hists)); + histb->SetMinimum(GetMinimumGT(*histb)); + + clk.Stop(); + clk.Print(); + TCanvas *c=new TCanvas; @@ -766,5 +1000,16 @@ c->cd(1); gPad->SetBorderMode(0); - gPad->Divide(1,2, 0, 0); + hists->Draw("colz"); + hists->SetBit(kCanDelete); + c->cd(2); + gPad->SetBorderMode(0); + hist->Draw("colz"); + hist->SetBit(kCanDelete); + c->cd(3); + gPad->SetBorderMode(0); + histb->Draw("colz"); + histb->SetBit(kCanDelete); + c->cd(4); + gPad->Divide(1,3, 0, 0); TVirtualPad *p=gPad; p->SetBorderMode(0); @@ -776,22 +1021,21 @@ gPad->SetBorderMode(0); h3.DrawCopy(); - c->cd(2); - gPad->SetBorderMode(0); - hist->Draw("colz"); - hist->SetDirectory(NULL); - hist->SetBit(kCanDelete); - c->cd(3); + p->cd(3); + gPad->SetBorderMode(0); + h2.DrawCopy(); + c->cd(6); + gPad->Divide(1,2, 0, 0); + TVirtualPad *q=gPad; + q->SetBorderMode(0); + q->cd(1); + gPad->SetBorderMode(0); gPad->SetBorderMode(0); h4b.DrawCopy(); h4a.DrawCopy("same"); h6.DrawCopy("same"); - c->cd(4); - gPad->SetBorderMode(0); - h2.DrawCopy(); - c->cd(6); - gPad->SetBorderMode(0); - h5b.DrawCopy(); - h5a.DrawCopy("same"); - + q->cd(2); + gPad->SetBorderMode(0); + //h5b.DrawCopy(); + //h5a.DrawCopy("same"); c->cd(5); gPad->SetBorderMode(0); @@ -801,5 +1045,7 @@ fHist.GetXaxis()->GetBinCenter(maxx), fHist.GetYaxis()->GetBinCenter(maxy), maxs); + TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy); + result->SetDirectory(NULL); result->SetNameTitle("Result \\alpha", title); result->SetBit(kCanDelete); @@ -808,8 +1054,8 @@ result->Draw(); - TF1 f1("MyFunc1", fcn, 0, 90, 5); - TF1 f2("MyFunc2", fcn, 0, 90, 5); - f1.SetParameters(par.GetArray()); - f2.SetParameters(par.GetArray()); + TF1 f1("", func.GetExpFormula(), 0, 90); + TF1 f2("", func.GetExpFormula(), 0, 90); + f1.SetParameters(maxpar.GetArray()); + f2.SetParameters(maxpar.GetArray()); f2.FixParameter(0, 0); f2.FixParameter(1, 0); @@ -824,14 +1070,52 @@ leg->SetBorderSize(1); leg->SetTextSize(0.04); - leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22); + leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22); + //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22); leg->AddLine(0, 0.65, 0, 0.65); - leg->AddText(0.06, 0.54, Form("A=%.2f", par[0]))->SetTextAlign(11); - leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", par[2]))->SetTextAlign(11); - leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fixed)", par[1]))->SetTextAlign(11); - leg->AddText(0.60, 0.54, Form("a=%.2f", par[3]))->SetTextAlign(11); - leg->AddText(0.60, 0.34, Form("b=%.2f", par[4]))->SetTextAlign(11); + leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11); + leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11); + leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11); + leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11); + leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11); + if (polynom>1) + leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11); leg->SetBit(kCanDelete); leg->Draw(); - } - -} + + q->cd(2); + + TGraph *g = new TGraph; + g->SetPoint(0, 0, 0); + + Int_t max=0; + Float_t maxsig=0; + for (int i=1; i<89; i++) + { + const Double_t s = f1.Integral(0, (float)i); + const Double_t b = f2.Integral(0, (float)i); + + const Double_t sig = Significance(s, b); + + g->SetPoint(g->GetN(), i, sig); + + if (sig>maxsig) + { + max = i; + maxsig = sig; + } + } + + g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha"); + g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]"); + g->GetHistogram()->SetYTitle("Significance"); + g->SetBit(kCanDelete); + g->Draw("AP"); + + leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC"); + leg->SetBorderSize(1); + leg->SetTextSize(0.1); + leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max)); + leg->SetBit(kCanDelete); + leg->Draw(); + } +} Index: trunk/MagicSoft/Mars/mhist/MHFalseSource.h =================================================================== --- trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3567) +++ trunk/MagicSoft/Mars/mhist/MHFalseSource.h (revision 3568) @@ -51,5 +51,5 @@ TH1 *GetHistByName(const TString name) { return &fHist; } - void FitSignificance(Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80); //*MENU* + void FitSignificance(Float_t sigint=15, Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80, Byte_t polynom=2); //*MENU* void FitSignificanceStd() { FitSignificance(); } //*MENU* Index: trunk/MagicSoft/Mars/mimage/MConcentration.cc =================================================================== --- trunk/MagicSoft/Mars/mimage/MConcentration.cc (revision 3567) +++ trunk/MagicSoft/Mars/mimage/MConcentration.cc (revision 3568) @@ -116,64 +116,48 @@ Int_t MConcentration::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hillas) { - const UInt_t npixevt = evt.GetNumPixels(); - Float_t maxpix[9] = {0,0,0,0,0,0,0,0,0}; // [#phot] - for (UInt_t i=0; iGetPixId(); + const Double_t nphot = pix->GetNumPhotons()* geom.GetPixRatio(pixid); // Get number of photons in the 8 most populated pixels + if (maxpix[0]<=nphot) + { + for(int i=0;i<8;i++) + maxpix[8-i]=maxpix[7-i]; + maxpix[0]=nphot; + continue; + } - if(maxpix[0]<=nphot) { - for(int i=0;i<8;i++) - maxpix[8-i]=maxpix[7-i]; - maxpix[0]=nphot; - } - else - for(int i=0;i<8;i++){ - if (nphot=maxpix[7-i]) + continue; + for(int j=0;j ! -! Copyright: MAGIC Software Development, 2000-2003 +! Copyright: MAGIC Software Development, 2000-2004 ! ! @@ -27,10 +27,50 @@ // MPointingPos // -// Store the current pointing position of the telescope... +// In this container we store the corrected pointing position of the +// telscope. The pointing coordinates are read into MReportDrive together +// with its time. +// +// MPointingPosCalc afterwards calculates corrections and checks for the +// cosistency of the coordinates. The result (the real coordinates) +// are stored in this container. No further correction should be necessary +// using MPointingPos. +// +// If you need the rotation angle of the starfield in the camera you can +// get it from here. // ///////////////////////////////////////////////////////////////////////////// #include "MPointingPos.h" +#include "MTime.h" +#include "MObservatory.h" +#include "MAstroSky2Local.h" + ClassImp(MPointingPos); using namespace std; + +// -------------------------------------------------------------------------- +// +// Get the corresponding rotation angle of the sky coordinate system +// seen with an Alt/Az telescope calculated from the stored local +// (Zd/Az) coordinates. +// +// For more information see MAstro::RotationAngle +// +Double_t MPointingPos::RotationAngle(const MObservatory &o) const +{ + return o.RotationAngle(fZd*TMath::DegToRad(), fAz*TMath::DegToRad()); +} + +// -------------------------------------------------------------------------- +// +// Get the corresponding rotation angle of the sky coordinate system +// seen with an Alt/Az telescope calculated from the stored sky +// (Ra/Dec) coordinates. +// +// For more information see MAstro::RotationAngle +// +Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t) const +{ + return MAstroSky2Local(t, o).RotationAngle(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad()); +} Index: trunk/MagicSoft/Mars/mpointing/MPointingPos.h =================================================================== --- trunk/MagicSoft/Mars/mpointing/MPointingPos.h (revision 3567) +++ trunk/MagicSoft/Mars/mpointing/MPointingPos.h (revision 3568) @@ -5,4 +5,7 @@ #include "MParContainer.h" #endif + +class MTime; +class MObservatory; class MPointingPos : public MParContainer @@ -32,4 +35,7 @@ Double_t GetDec() const { return fDec; } + Double_t RotationAngle(const MObservatory &o) const; + Double_t RotationAngle(const MObservatory &o, const MTime &t) const; + ClassDef(MPointingPos, 1) //Container storing the (corrected) telescope pointing position }; Index: trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (revision 3568) +++ trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.cc (revision 3568) @@ -0,0 +1,182 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 3/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSrcPosCalc +// +// Calculate the current source position in the camera from the (already +// corrected) local pointing position of the telescope (MPointingPos) by +// derotating the position (given in x and y coordinates in the camera +// plane) at culmination time by the current rotation angle of the +// starfield dtermined from MObservatory and MPointingPos +// +// The conversion factor between the camera plain (mm) and the +// sky (deg) is taken from MGeomCam. +// +// Input Container: +// MPointingPos +// MObservatory +// MGeomCam +// +// Output Container: +// MSrcPosCam +// +// To be done: +// - wobble mode missing +// +////////////////////////////////////////////////////////////////////////////// +#include "MSrcPosCalc.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MRawRunHeader.h" +#include "MObservatory.h" +#include "MSrcPosCam.h" +#include "MAstroSky2Local.h" +#include "MPointingPos.h" +#include "MGeomCam.h" + +ClassImp(MSrcPosCalc); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Initialize fY and fY with 0 +// +MSrcPosCalc::MSrcPosCalc(const char *name, const char *title) + : fX(0), fY(0) +{ + fName = name ? name : "MSrcPosCalc"; + fTitle = title ? title : "Derotates the source position in the camera"; +} + +// -------------------------------------------------------------------------- +// +// Search and if necessary create MSrcPosCam in the parameter list +// +Int_t MSrcPosCalc::PreProcess(MParList *pList) +{ + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeom) + { + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; + } + + fPointPos = (MPointingPos*)pList->FindObject("MPointingPos"); + if (!fPointPos) + { + *fLog << err << "MPointPos not found... aborting." << endl; + return kFALSE; + } + + fSrcPos = (MSrcPosCam*)pList->FindCreateObj("MSrcPosCam"); + if (!fSrcPos) + return kFALSE; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Search the parameter list for MObservatory and MTime +// +Bool_t MSrcPosCalc::ReInit(MParList *pList) +{ + fObservatory = (MObservatory*)pList->FindObject("MObservatory"); + if (!fObservatory) + { + *fLog << err << "MObservatory not found... aborting." << endl; + return kFALSE; + } + + if (fX!=0 || fY!=0) + return kTRUE; + + *fLog << warn << "fX==0 && fY==0: Using arbitrary source position!" << endl; + + // + // This is a workaround for current crab misspointing - DO NOT USE! + // It will be removed soon! + // + const MRawRunHeader *h = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); + if (!h) + { + *fLog << err << "MRawRunHeader not found... aborting." << endl; + return kFALSE; + } + + const MTime &t = h->GetRunStart(); + + const Double_t rho = fPointPos->RotationAngle(*fObservatory, t); + + // Calculate x, y of Zeta Tauri + + Double_t tm = t.GetAxisTime(); + + Double_t x = 1.7267e6-6.03285e-3*tm; // [mm] + Double_t y = -189.823+974.908*exp(-52.3083*(tm*1e-5-2861.5)); // [mm] + + const Double_t cs = TMath::Cos(rho-fDrho); + const Double_t sn = TMath::Sin(rho-fDrho); + + const Double_t dx = fR*cs; + const Double_t dy = fR*sn; + + fSrcPos->SetXY(x-dx, y-dy); + + *fLog << dbg << t << " - Position: x=" << x-dx << "mm y=" << y-dy << "mm" << endl; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Derotate fX/fY by the current rotation angle, set MSrcPosCam +// +Int_t MSrcPosCalc::Process() +{ + if (fX==0 && fY==0) + return kTRUE; + + // Define source position in the camera plain + TVector2 v(fX, fY); + + // rotate the source position by the current rotation angle + const Double_t rho = fPointPos->RotationAngle(*fObservatory); + v.Rotate(-rho); + + // Convert coordinates into camera plain (mm) + v *= 1./fGeom->GetConvMm2Deg(); + + // Set current source position + fSrcPos->SetXY(v); + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.h =================================================================== --- trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.h (revision 3568) +++ trunk/MagicSoft/Mars/mpointing/MSrcPosCalc.h (revision 3568) @@ -0,0 +1,44 @@ +#ifndef MARS_MSrcPosCalc +#define MARS_MSrcPosCalc + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MObservatory; +class MPointingPos; +class MSrcPosCam; +class MGeomCam; + +class MSrcPosCalc : public MTask +{ +private: + MObservatory *fObservatory; + MPointingPos *fPointPos; + MSrcPosCam *fSrcPos; + MGeomCam *fGeom; + + Double_t fR; // Distance of source to a fitted star + Double_t fDrho; // angle between the two stars + + Double_t fX; // [deg] x-position of source + Double_t fY; // [deg] y-position of source + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Bool_t ReInit(MParList *pList); + +public: + MSrcPosCalc(const char *name=NULL, const char *title=NULL); + + void SetOffset(Double_t r, Double_t drho) + { + fDrho = drho*TMath::DegToRad(); + fR = r; + } + void SetPositionXY(Double_t x, Double_t y) { fX=x; fY=y; } + + ClassDef(MSrcPosCalc, 0) // Derotates the source position in the camera +}; + +#endif Index: trunk/MagicSoft/Mars/mpointing/MSrcPosCam.cc =================================================================== --- trunk/MagicSoft/Mars/mpointing/MSrcPosCam.cc (revision 3568) +++ trunk/MagicSoft/Mars/mpointing/MSrcPosCam.cc (revision 3568) @@ -0,0 +1,120 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Thomas Bretz 12/2000 +! Author(s): Rudolf Bock 10/2001 +! +! Copyright: MAGIC Software Development, 2000-2002 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MSrcPosCam +// +// Storage Container to hold the current position of the source (or +// anti/false source) in the camera plain +// +////////////////////////////////////////////////////////////////////////////// +#include "MSrcPosCam.h" + +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MSrcPosCam); + +using namespace std; + +static const TString gsDefName = "MSrcPosCam"; +static const TString gsDefTitle = "Virtual source position in the camera"; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MSrcPosCam::MSrcPosCam(const char *name, const char *title) : fX(0), fY(0) +{ + fName = name ? name : gsDefName.Data(); + fTitle = title ? title : gsDefTitle.Data(); +} + +// ----------------------------------------------------------------------- +// +void MSrcPosCam::Print(Option_t *) const +{ + *fLog << all; + *fLog << "Source position in the camera plain (" << GetName() << ")" << endl; + *fLog << " - x [mm] = " << fX << endl; + *fLog << " - y [mm] = " << fY << endl; +} + +void MSrcPosCam::SetXY(const TVector2 &v) +{ + fX = v.X(); + fY = v.Y(); +} + +TVector2 MSrcPosCam::GetXY() const +{ + return TVector2(fX, fY); +} + +/* +// ----------------------------------------------------------------------- +// +// overloaded MParContainer to read MSrcPosCam from an ascii file +// +void MSrcPosCam::AsciiRead(ifstream &fin) +{ + fin >> fX; + fin >> fY; +} + +// ----------------------------------------------------------------------- +// +// overloaded MParContainer to write MSrcPosCam to an ascii file +// +void MSrcPosCam::AsciiWrite(ofstream &fout) const +{ + fout << fX << " " << fY; +} +*/ + +// -------------------------------------------------------------------------- +// +// Implementation of SavePrimitive. Used to write the call to a constructor +// to a macro. In the original root implementation it is used to write +// gui elements to a macro-file. +// +void MSrcPosCam::StreamPrimitive(ofstream &out) const +{ + out << " MSrcPosCam " << GetUniqueName(); + if (fName!=gsDefName) + { + out << "(\"" << fName << "\""; + if (fTitle!=gsDefTitle) + out << ", \"" << fTitle << "\""; + out <<")"; + } + out << ";" << endl; + + out << " " << GetUniqueName() << ".SetXY(" << fX << ", " << fY << ");" << endl;} Index: trunk/MagicSoft/Mars/mpointing/MSrcPosCam.h =================================================================== --- trunk/MagicSoft/Mars/mpointing/MSrcPosCam.h (revision 3568) +++ trunk/MagicSoft/Mars/mpointing/MSrcPosCam.h (revision 3568) @@ -0,0 +1,40 @@ +#ifndef MARS_MSrcPosCam +#define MARS_MSrcPosCam + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +class TVector2; + +class MSrcPosCam : public MParContainer +{ +private: + Float_t fX; // [mm] x position of source in camera + Float_t fY; // [mm] y position of source in camera + + void StreamPrimitive(ofstream &out) const; + +public: + MSrcPosCam(const char *name=NULL, const char *title=NULL); + + void Clear(Option_t *) { fX = 0; fY = 0; } + + void SetX(Float_t x) { fX = x; } + void SetY(Float_t y) { fY = y; } + void SetXY(Float_t x, Float_t y) { fX = x; fY = y; } + void SetXY(const TVector2 &v); + + Float_t GetX() const { return fX; } + Float_t GetY() const { return fY; } + TVector2 GetXY() const; + + void Print(Option_t *opt=NULL) const; + + //void AsciiRead(ifstream &fin); + //void AsciiWrite(ofstream &fout) const; + + ClassDef(MSrcPosCam, 1) // container to store source position in the camera plain +}; + +#endif Index: trunk/MagicSoft/Mars/mpointing/Makefile =================================================================== --- trunk/MagicSoft/Mars/mpointing/Makefile (revision 3567) +++ trunk/MagicSoft/Mars/mpointing/Makefile (revision 3568) @@ -22,5 +22,5 @@ # connect the include files defined in the config.mk file # -INCLUDES = -I. -I../mbase -I../mraw -I../mreport -I../mmc +INCLUDES = -I. -I../mbase -I../mraw -I../mreport -I../mmc -I../mastro #------------------------------------------------------------------------------ Index: trunk/MagicSoft/Mars/mraw/MRawFileRead.cc =================================================================== --- trunk/MagicSoft/Mars/mraw/MRawFileRead.cc (revision 3567) +++ trunk/MagicSoft/Mars/mraw/MRawFileRead.cc (revision 3568) @@ -38,7 +38,7 @@ // // ////////////////////////////////////////////////////////////////////////////// - #include "MRawFileRead.h" +#include #include @@ -138,5 +138,6 @@ if (!(*fIn)) { - *fLog << err << "Error: Cannot open file '" << fFileName << "'" << endl; + *fLog << err << "Cannot open file " << fFileName << ": "; + *fLog << strerror(errno) << endl; return kFALSE; } Index: trunk/MagicSoft/Mars/mreport/MReportFileRead.cc =================================================================== --- trunk/MagicSoft/Mars/mreport/MReportFileRead.cc (revision 3567) +++ trunk/MagicSoft/Mars/mreport/MReportFileRead.cc (revision 3568) @@ -40,4 +40,5 @@ #include "MReportFileRead.h" +#include #include @@ -212,7 +213,9 @@ if (!(*fIn)) { - *fLog << err << "Error: Cannot open file '" << fFileName << "'" << endl; + *fLog << err << "Cannot open file " << fFileName << ": "; + *fLog << strerror(errno) << endl; return kFALSE; } + if (TestBit(kHasNoHeader)) return kTRUE;