Changeset 3568


Ignore:
Timestamp:
03/22/04 09:25:49 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
4 added
25 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3567 r3568  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2004/03/22: Thomas Bretz
     22 
     23   * mpointing/MSrcPosCalc.[h,cc], MSrcPosCam.[h,cc]:
     24     - added
     25     
     26   * mastro/MAstro.[h,cc]:
     27     - added code to calculate rotationangle previously in MObservatory
     28     - changed definition of rotation angle such, that it is now
     29       180deg if Ra and Az grid is parallel
     30
     31   * mastro/MAstroCatalog.[h,cc]:
     32     - fixes and enhancements to the display (such as misscalculated
     33       number of grid lines, title display, etc)
     34     - enhancements to the output
     35     - generalized creation of grid - for further usage
     36
     37   * mastro/MAstroSky2Local.[h,cc]:
     38     - replaced calculation of rotation angle by the function in
     39       MAstro
     40
     41   * mastro/MObservatory.[h,cc]:
     42     - small changes to Print output
     43     - moved code for calculation of rotation angle to MAstro
     44
     45   * mbase/MEvtLoop.cc:
     46     - do not output number of events per second if no events processed
     47     
     48   * mbase/MParList.cc:
     49     - updated some comments
     50     
     51   * mfileio/MCT1ReadAscii.cc, mfileio/MCT1ReadPreProc.cc,
     52     mfileio/MReadRflFile.cc, mraw/MRawFileRead.cc,
     53     mreport/MReportFileRead.cc:
     54     - output error string if file cannot be opened
     55     
     56   * mfileio/MReadTree.cc:
     57     - output name of chain which is scanned
     58
     59   * mimage/MConcentration.cc:
     60     - replaced loop by iterator
     61     - removed obsolete (unused) variables
     62     
     63   * mimage/MHNewImagePar.[h,cc]:
     64     - fixed display colors
     65
     66   * mpointing/MPointingPos.[h,cc]:
     67     - added member function to calculate rotation angle
     68     - added comments
     69
     70   * mpointing/Makefile:
     71     - added include MAstro
     72
     73
     74
    2075 2004/03/19: Markus Gaug
    2176
    2277   * mcalib/MHCalibrationChargePix.cc
    23      - added some style sto the default Draw in order to see the
     78     - added some style to the default Draw in order to see the
    2479       label and axis titles better
    2580
    2681   * mcalib/MHCalibrationChargeCam.[h,cc]
    27      - store and display more informaiton on the average pxiels
     82     - store and display more information on the average pxiels
    2883
    2984   * mcalib/MCalibrationCam.cc
    … …  
    4297     - included MCalibrationQECam to be initialized
    4398
    44    * mcalib/MCalibrationChargePix.[h,cc]
    45    * mcalib/MCalibrationQEPix.[h,cc]
     99   * mcalib/MCalibrationChargePix.[h,cc],
     100     mcalib/MCalibrationQEPix.[h,cc]:
    46101     - replace DefinePixId by SetPixId
    47102   
    … …  
    59114     - remove uncommented piece of code
    60115
    61    * msignal/MExtractSignal.cc
    62    * msignal/MExtractSignal2.cc
     116   * msignal/MExtractSignal.cc, msignal/MExtractSignal2.cc:
    63117     - remove warning about pixel with low gain saturation,
    64118       now in MBadPixelsPix
    65119
    66    * mbadpixels/MBadPixelsPix.[h,cc]
    67    * mcalib/MCalibrationChargeCam.cc
     120   * mbadpixels/MBadPixelsPix.[h,cc], mcalib/MCalibrationChargeCam.cc:
    68121     - added new flag: kDeviatingNumPhes
    69122
    70123   * mcalib/MCalibrationChargePix.cc
    71      - check for mean arr. time in last bin replaced by check in last two bins
    72 
    73    * mcalib/MCalibrationChargePix.[h,cc]
    74    * mcalib/MCalibrationChargeCam.cc
    75    * mcalib/MHCalibrationChargeCam.cc
    76      - removed flag kHiGainFitted, kLoGainFitted, since they are available
    77        from MBadPixelsPix
    78 
    79    * macros/calibration.C
    80    * macros/calibrate_data.C
     124     - check for mean arr. time in last bin replaced by check in last
     125       two bins
     126
     127   * mcalib/MCalibrationChargePix.[h,cc],
     128     mcalib/MCalibrationChargeCam.cc,
     129     mcalib/MHCalibrationChargeCam.cc:
     130     - removed flag kHiGainFitted, kLoGainFitted, since they are
     131       available from MBadPixelsPix
     132
     133   * macros/calibration.C, macros/calibrate_data.C
    81134     - a few flags from MCalibrationChargeCam::GetPixelContent were wrong,
    82135       corrected them
     136
    83137
    84138

  • trunk/MagicSoft/Mars/mastro/MAstro.cc

    r3537 r3568  
    310310    return MTime(ut1).GetGmst();
    311311}
     312
     313// --------------------------------------------------------------------------
     314//
     315// RotationAngle
     316//
     317// calculates the angle for the rotation of the sky coordinate system
     318// with respect to the local coordinate system. This is identical
     319// to the rotation angle of the sky image in the camera.
     320//
     321//  sinl  [rad]: sine of observers latitude
     322//  cosl  [rad]: cosine of observers latitude
     323//  theta [rad]: polar angle/zenith distance
     324//  phi   [rad]: rotation angle/azimuth
     325//
     326// Return sin/cos component of angle
     327//
     328// The convention is such, that the rotation angle is -pi/pi if
     329// right ascension and local rotation angle are counted in the
     330// same direction, 0 if counted in the opposite direction.
     331//
     332// (In other words: The rotation angle is 0 when the source culminates)
     333//
     334// Using vectors it can be done like:
     335//    TVector3 v, p;
     336//    v.SetMagThetaPhi(1, theta, phi);
     337//    p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
     338//    v = v.Cross(l));
     339//    v.RotateZ(-phi);
     340//    v.Rotate(-theta)
     341//    rho = TMath::ATan2(v(2), v(1));
     342//
     343// For more information see TDAS 00-11, eqs. (18) and (20)
     344//
     345void MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos)
     346{
     347    const Double_t sint = TMath::Sin(theta);
     348    const Double_t cost = TMath::Cos(theta);
     349
     350    const Double_t snlt = sinl*sint;
     351    const Double_t cslt = cosl*cost;
     352
     353    const Double_t sinp = TMath::Sin(phi);
     354    const Double_t cosp = TMath::Cos(phi);
     355
     356    const Double_t v1 = sint*sinp;
     357    const Double_t v2 = cosl - snlt*cosp;
     358
     359    const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2);
     360
     361    cos =  -cosl*sinp      / denom;
     362    sin = (snlt-cslt*cosp) / denom;
     363}
     364
     365// --------------------------------------------------------------------------
     366//
     367// RotationAngle
     368//
     369// calculates the angle for the rotation of the sky coordinate system
     370// with respect to the local coordinate system. This is identical
     371// to the rotation angle of the sky image in the camera.
     372//
     373//  sinl  [rad]: sine of observers latitude
     374//  cosl  [rad]: cosine of observers latitude
     375//  theta [rad]: polar angle/zenith distance
     376//  phi   [rad]: rotation angle/azimuth
     377//
     378// Return angle [rad] in the range -pi, pi
     379//
     380// The convention is such, that the rotation angle is -pi/pi if
     381// right ascension and local rotation angle are counted in the
     382// same direction, 0 if counted in the opposite direction.
     383//
     384// (In other words: The rotation angle is 0 when the source culminates)
     385//
     386// Using vectors it can be done like:
     387//    TVector3 v, p;
     388//    v.SetMagThetaPhi(1, theta, phi);
     389//    p.SetMagThetaPhi(1, TMath::Pi()/2-latitude, 0);
     390//    v = v.Cross(l));
     391//    v.RotateZ(-phi);
     392//    v.Rotate(-theta)
     393//    rho = TMath::ATan2(v(2), v(1));
     394//
     395// For more information see TDAS 00-11, eqs. (18) and (20)
     396//
     397Double_t MAstro::RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi)
     398{
     399    const Double_t snlt = sinl*TMath::Sin(theta);
     400    const Double_t cslt = cosl*TMath::Cos(theta);
     401
     402    const Double_t sinp = TMath::Sin(phi);
     403    const Double_t cosp = TMath::Cos(phi);
     404
     405    return TMath::ATan2(-cosl*sinp, snlt-cslt*cosp);
     406}

  • trunk/MagicSoft/Mars/mastro/MAstro.h

    r3371 r3568  
    5050    static Double_t UT2GMST(Double_t ut1);
    5151
     52    // Rotation angle between local and sky coordinate system
     53    static void     RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi, Double_t &sin, Double_t &cos);
     54    static Double_t RotationAngle(Double_t sinl, Double_t cosl, Double_t theta, Double_t phi);
     55
    5256    ClassDef(MAstro, 0)
    5357};

  • trunk/MagicSoft/Mars/mastro/MAstroCamera.cc

    r3537 r3568  
    5757using namespace std;
    5858
     59// --------------------------------------------------------------------------
    5960MAstroCamera::MAstroCamera() : fGeom(0), fMirrors(0)
    6061{
    … …  
    6364}
    6465
     66// --------------------------------------------------------------------------
    6567MAstroCamera::~MAstroCamera()
    6668{
    … …  
    7375}
    7476
    75 void MAstroCamera::SetMirrors(TClonesArray *arr)
    76 {
    77     if (!arr || arr->GetClass()!=MGeomMirror::Class())
     77// --------------------------------------------------------------------------
     78void MAstroCamera::SetMirrors(TClonesArray &arr)
     79{
     80    if (arr.GetClass()!=MGeomMirror::Class())
    7881        return;
    7982
    80     const Int_t n = arr->GetSize();
     83    const Int_t n = arr.GetSize();
    8184
    8285    if (!fMirrors)
    … …  
    8689
    8790    for (int i=0; i<n; i++)
    88         memcpy((*fMirrors)[i], (*arr)[i], sizeof(MGeomMirror));
    89 
    90 }
    91 
     91        memcpy((*fMirrors)[i], arr[i], sizeof(MGeomMirror));
     92
     93}
     94
     95// --------------------------------------------------------------------------
    9296void MAstroCamera::SetGeom(const MGeomCam &cam)
    9397{
    … …  
    98102}
    99103
    100 Int_t MAstroCamera::ConvertToPad(const TVector3 &w, TVector2 &v)
    101 {
    102     const TVector3 spot = fMirror0->GetReflection(w, fGeom->GetCameraDist())*1000;
    103 
     104// --------------------------------------------------------------------------
     105Int_t MAstroCamera::ConvertToPad(const TVector3 &w, TVector2 &v) const
     106{
    104107    /*
    105108     --- Use this to plot the 'mean grid' instead of the grid of a
    … …  
    113116        */
    114117
    115 
     118    const TVector3 spot = fMirror0->GetReflection(w, fGeom->GetCameraDist())*1000;
    116119    v.Set(spot(0), spot(1));
    117120
    … …  
    120123}
    121124
    122 void MAstroCamera::DrawNet(const TRotation &trans)
    123 {
    124     const TRotation rot(MAstroSky2Local(*fTime, *fObservatory));
    125 
    126     TVector2 radec(fRaDec.Phi(), fRaDec.Theta());
    127     MAstroCatalog::DrawNet(radec, trans*rot, 2);
    128 
    129     const TVector3 zdaz0 = MAstroSky2Local(*fTime, *fObservatory)*fRaDec;
    130     TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta());
    131     MAstroCatalog::DrawNet(zdaz, trans, 1);
    132 }
    133 
     125// --------------------------------------------------------------------------
    134126TObject *FindObjectInPad(const char *name, TVirtualPad *pad)
    135127{
    … …  
    155147}
    156148
     149// --------------------------------------------------------------------------
    157150void MAstroCamera::AddPrimitives(Option_t *o)
    158151{
    … …  
    171164    const Bool_t hasdot  = opt.Contains(".", TString::kIgnoreCase);
    172165    const Bool_t usecam  = opt.Contains("c", TString::kIgnoreCase);
    173 
    174     MAstroSky2Local rot(*fTime, *fObservatory);
    175 
    176     const Float_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta());
    177 
    178     TString str = fTime->GetSqlDateTime();
    179     str += Form("  (\\alpha=%.1fh \\delta=%.1f\\circ)  \\rho=%.1f\\circ",
    180                 fRaDec.Phi()/TMath::Pi()*12, 90-fRaDec.Theta()*TMath::RadToDeg(),
    181                 rho *TMath::RadToDeg());
    182166
    183167    // Get camera
    … …  
    198182    }
    199183
    200     camera->SetTitle(str);
     184    camera->SetTitle(GetPadTitle());
    201185
    202186    gPad->cd(1);
    … …  
    221205    }
    222206
    223     TVector3 zdaz0 = fRaDec;
    224     zdaz0 *= rot;
    225 
    226     cout << zdaz0.Phi()*TMath::RadToDeg() << " " << zdaz0.Theta()*TMath::RadToDeg() << endl;
    227 
    228     TVector3 test = zdaz0;
    229     test *= rot.Inverse();
    230 
    231     cout << test.Phi()*TMath::RadToDeg()/15 << " " << test.Theta()*TMath::RadToDeg() << endl;
    232 
    233     TRotation rot2;
    234     rot2.RotateZ(-zdaz0.Phi());
    235     rot2.RotateY(-zdaz0.Theta());
    236     rot2.RotateZ(-TMath::Pi()/2); // align coordinate system
    237 
    238     DrawNet(rot2);
     207    const TRotation rot(GetGrid(kTRUE));
    239208
    240209    MVector3 *radec;
    … …  
    247216        TVector3 star(*radec);
    248217
    249         // Calculate local coordinates
     218        // Rotate Star into telescope system
    250219        star *= rot;
    251         // Rotate Star into telescope system
    252         star *= rot2;
    253220
    254221        TVector3 mean;
    … …  
    270237            {
    271238                TMarker *m=new TMarker(spot(0), spot(1), 1);
    272                 m->SetBit(kCannotPick);
    273                 m->SetBit(kCanDelete);
    274239                m->SetMarkerColor(kMagenta);
    275240                m->SetMarkerStyle(kDot);
    … …  
    303268        case kKey_Plus:
    304269            fTime->SetMjd(fTime->GetMjd()+0.25/24);
    305             SetBit(kHasChanged);
    306             gPad->Modified();
    307             gPad->Update();
     270            Update(kTRUE);
    308271            return;
    309272
    310273        case kKey_Minus:
    311274            fTime->SetMjd(fTime->GetMjd()-0.25/24);
    312             SetBit(kHasChanged);
    313             gPad->Modified();
    314             gPad->Update();
     275            Update(kTRUE);
    315276            return;
    316277        }

  • trunk/MagicSoft/Mars/mastro/MAstroCamera.h

    r3537 r3568  
    2121    MGeomMirror  *fMirror0;     //!
    2222
    23     Int_t  ConvertToPad(const TVector3 &w, TVector2 &v);
     23    Int_t  ConvertToPad(const TVector3 &w, TVector2 &v) const;
    2424    void   AddPrimitives(Option_t *o);
    2525    void   SetRangePad() { }
    26     void   DrawNet(const TRotation &rot);
    2726    void   ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2);
    2827
    … …  
    3130    ~MAstroCamera();
    3231
    33     void SetMirrors(TClonesArray *arr);
     32    void SetMirrors(TClonesArray &arr);
    3433    void SetGeom(const MGeomCam &cam);
    3534

  • trunk/MagicSoft/Mars/mastro/MAstroCatalog.cc

    r3537 r3568  
    3333#include "MAstroCatalog.h"
    3434
     35#include <errno.h>
    3536#include <fstream>
    3637#include <stdlib.h>
     38#include <limits.h> // INT_MAX (Suse 7.3/gcc 2.95)
    3739
    3840#include <KeySymbols.h>
    … …  
    4547#include <TGToolTip.h>
    4648#include <TRotation.h>
     49#include <TPaveText.h>
    4750#include <TStopwatch.h>
    4851
    … …  
    5962
    6063using namespace std;
    61 
     64/*
    6265class MRotation : public TRotation
    6366{
    … …  
    7679    }
    7780};
    78 
     81*/
    7982/*
    8083MVector3 MVector3::GetZdAz(const MObservatory &obs, Double_t gmst) const
    … …  
    134137}
    135138*/
    136 MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(99), fToolTip(0), fObservatory(0), fTime(0)
     139// --------------------------------------------------------------------------
     140MAstroCatalog::MAstroCatalog() : fLimMag(99), fRadiusFOV(90), fToolTip(0), fObservatory(0), fTime(0)
    137141{
    138142    fList.SetOwner();
    … …  
    140144}
    141145
     146// --------------------------------------------------------------------------
    142147MAstroCatalog::~MAstroCatalog()
    143148{
    … …  
    161166}
    162167
     168// --------------------------------------------------------------------------
    163169TString MAstroCatalog::FindToken(TString &line, Char_t tok)
    164170{
    … …  
    176182}
    177183
     184// --------------------------------------------------------------------------
    178185Int_t MAstroCatalog::atoi(const TSubString &sub)
    179186{
    … …  
    181188}
    182189
     190// --------------------------------------------------------------------------
    183191Float_t MAstroCatalog::atof(const TSubString &sub)
    184192{
    … …  
    186194}
    187195
     196// --------------------------------------------------------------------------
    188197Int_t MAstroCatalog::atoi(const TString &s)
    189198{
    … …  
    191200}
    192201
     202// --------------------------------------------------------------------------
    193203Float_t MAstroCatalog::atof(const TString &s)
    194204{
    … …  
    196206}
    197207
     208// --------------------------------------------------------------------------
    198209Int_t MAstroCatalog::ReadXephem(TString catalog)
    199210{
    … …  
    203214
    204215    ifstream fin(catalog);
     216    if (!fin)
     217    {
     218        gLog << err << "Cannot open file " << catalog << ": ";
     219        gLog << strerror(errno) << endl;
     220        return 0;
     221    }
    205222
    206223    Int_t add =0;
    … …  
    276293}
    277294
     295// --------------------------------------------------------------------------
    278296Int_t MAstroCatalog::ReadNGC2000(TString catalog)
    279297{
    … …  
    283301
    284302    ifstream fin(catalog);
     303    if (!fin)
     304    {
     305        gLog << err << "Cannot open file " << catalog << ": ";
     306        gLog << strerror(errno) << endl;
     307        return 0;
     308    }
    285309
    286310    Int_t add=0;
    … …  
    339363}
    340364
     365// --------------------------------------------------------------------------
    341366Int_t MAstroCatalog::ReadBSC(TString catalog)
    342367{
    … …  
    346371
    347372    ifstream fin(catalog);
     373    if (!fin)
     374    {
     375        gLog << err << "Cannot open file " << catalog << ": ";
     376        gLog << strerror(errno) << endl;
     377        return 0;
     378    }
    348379
    349380    Int_t add=0;
    … …  
    404435}
    405436
     437// --------------------------------------------------------------------------
    406438void MAstroCatalog::Paint(Option_t *o)
    407439{
    408 //    if (!gPad->IsBatch())
    409 //        gVirtualX->ClearWindow();
     440    SetRangePad();
    410441
    411442    if (TestBit(kHasChanged))
    … …  
    413444}
    414445
     446// --------------------------------------------------------------------------
    415447void MAstroCatalog::DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t transparent, const char *txt)
    416448{
    … …  
    428460    TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotLarge);;
    429461    tip->SetMarkerColor(kBlack);
    430     tip->SetBit(kCanDelete);
    431     tip->SetBit(kCannotPick);
    432462    AddMap(tip, new TString(str));
    433463}
    434464
    435 void MAstroCatalog::Update()
     465// --------------------------------------------------------------------------
     466void MAstroCatalog::Update(Bool_t upd)
    436467{
    437468    if (gPad && TestBit(kMustCleanup))
    … …  
    439470        SetBit(kHasChanged);
    440471        gPad->Modified();
    441     }
    442 }
    443 
     472        if (upd)
     473            gPad->Update();
     474    }
     475}
     476
     477// --------------------------------------------------------------------------
     478//
     479// Set the observation time. Necessary to use local coordinate
     480// system.
     481//
    444482void MAstroCatalog::SetTime(const MTime &time)
    445483{
    … …  
    449487}
    450488
     489// --------------------------------------------------------------------------
     490//
     491// Set the observatory location. Necessary to use local coordinate
     492// system.
     493//
    451494void MAstroCatalog::SetObservatory(const MObservatory &obs)
    452495{
    … …  
    456499}
    457500
    458 Int_t MAstroCatalog::ConvertToPad(const TVector3 &w0, TVector2 &v)
     501// --------------------------------------------------------------------------
     502//
     503// Convert the vector to pad coordinates. After conversion
     504// the x- coordinate of the vector must be the x coordinate
     505// of the pad - the same for y. If the coordinate is inside
     506// the current draw area return kTRUE, otherwise kFALSE.
     507// If it is an invalid coordinate return kERROR
     508//
     509Int_t MAstroCatalog::ConvertToPad(const TVector3 &w0, TVector2 &v) const
    459510{
    460511    TVector3 w(w0);
    461512
    462     // Stretch such, that the X-component is alwas the same. Now
    463     // Y and Z contains the crossing point between the star-light
    464     // and the plain of a virtual screen (ccd...)
     513    // Stretch such, that the Z-component is alwas the same. Now
     514    // X and Y contains the intersection point between the star-light
     515    // and the plain of a virtual plain screen (ccd...)
    465516    if (TestBit(kPlainScreen))
    466         w *= TMath::RadToDeg()/w.X();
    467     else
    468         w.SetMag(TMath::RadToDeg());
    469 
    470     v.Set(w(1), w(2));
    471 
    472     if (w(0)<0)
     517        w *= 1./w(2);
     518
     519    w *= TMath::RadToDeg();
     520    v.Set(w(0), w(1));
     521
     522    if (w(2)<0)
    473523        return kERROR;
    474524
    475     return w(1)>gPad->GetX1() && w(2)>gPad->GetY1() &&
    476            w(1)<gPad->GetX2() && w(2)<gPad->GetY2();
    477 }
    478 
    479 Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v)
     525    return w(0)>gPad->GetX1() && w(1)>gPad->GetY1() &&
     526           w(0)<gPad->GetX2() && w(1)<gPad->GetY2();
     527}
     528
     529// --------------------------------------------------------------------------
     530Int_t MAstroCatalog::Convert(const TRotation &rot, TVector2 &v) const
    480531{
    481532    MVector3 w;
    … …  
    486537}
    487538
     539// --------------------------------------------------------------------------
    488540Bool_t MAstroCatalog::DrawLine(const TVector2 &v, Double_t dx, Double_t dy, const TRotation &rot, Int_t type)
    489541{
    … …  
    501553
    502554    TLine *line = new TLine(v0.X(), v0.Y(), v1.X(), v1.Y());
    503     line->SetBit(kCanDelete);
    504555    line->SetLineStyle(kDashDotted); //kDashed, kDotted, kDashDotted
    505556    line->SetLineColor(kWhite+type*2);
    506     line->SetBit(kCannotPick);
    507557    AddMap(line);
    508558
    509559    const TVector2 deg = v*TMath::RadToDeg();
    510560    TString txt = type==1 ?
    511         Form("Ra=%.1fh  Dec=%.1fd", fmod(deg.X()/15+48, 24),  fmod(90-deg.Y()+270,180)-90) :
     561        Form("Ra=%.2fh  Dec=%.1fd", fmod(deg.X()/15+48, 24),  fmod(90-deg.Y()+270,180)-90) :
    512562        Form("Zd=%.1fd  Az=%.1fd",  fmod(deg.Y()+270,180)-90, fmod(deg.X()+720, 360));
    513563
    514564    TMarker *tip=new TMarker(v0.X(), v0.Y(), kDot);
    515     tip->SetBit(kCanDelete);
    516     tip->SetBit(kCannotPick);
    517565    tip->SetMarkerColor(kWhite+type*2);
    518566    AddMap(tip, new TString(txt));
    … …  
    521569}
    522570
    523 
     571// --------------------------------------------------------------------------
     572//
     573// Use "local" draw option to align the display to the local
     574// coordinate system instead of the sky coordinate system.
     575//
    524576void MAstroCatalog::Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type)
    525577{
    … …  
    570622}
    571623
    572 void MAstroCatalog::DrawNet(const TVector2 &v0, const TRotation &rot, Int_t type)
     624// --------------------------------------------------------------------------
     625void MAstroCatalog::DrawGrid(const TVector3 &v0, const TRotation &rot, Int_t type)
    573626{
    574627    TArrayI dx(1);
    … …  
    576629
    577630    // align to 1deg boundary
    578     TVector2 v = v0*TMath::RadToDeg();
     631    TVector2 v(v0.Phi()*TMath::RadToDeg(), v0.Theta()*TMath::RadToDeg());
    579632    v.Set((Float_t)TMath::Nint(v.X()), (Float_t)TMath::Nint(v.Y()));
    580633
    581634    // calculate stepsizes based on visible FOV
    582     Int_t stepx=1;
    583 
    584     if (fabs(90-v.Y())>90-fRadiusFOV || fabs(90-v.Y())<fRadiusFOV)
    585         stepx = 180/10;
     635    Int_t stepx = 1;
     636
     637    if (v.Y()<fRadiusFOV || v.Y()>180-fRadiusFOV)
     638        stepx=36;
    586639    else
    587640    {
    588641        // This is a rough estimate how many degrees are visible
    589         const Float_t m = log(fRadiusFOV/180.)/log(90./fRadiusFOV-1);
     642        const Float_t m = log(fRadiusFOV/180.)/log(90./(fRadiusFOV+1)+1);
    590643        const Float_t t = log(180.)-m*log(fRadiusFOV);
    591         const Int_t n = (Int_t)(exp(m*log(90-fabs(90-v.Y()))+t)+0.5);
    592         stepx = n<6 ? 1 : n/6;
    593     }
    594 
    595     const Int_t n = (Int_t)(fRadiusFOV+1);
    596     Int_t stepy = n<4 ? 1 : n/4;
     644        const Float_t f = m*log(90-fabs(90-v.Y()))+t;
     645        const Int_t nx = (Int_t)(exp(f)+0.5);
     646        stepx = nx<4 ? 1 : nx/4;
     647        if (stepx>36)
     648            stepx=36;
     649    }
     650
     651    const Int_t ny = (Int_t)(fRadiusFOV+1);
     652    Int_t stepy = ny<4 ? 1 : ny/4;
    597653
    598654    // align stepsizes to be devisor or 180 and 90
    … …  
    602658        stepy++;
    603659
    604     // align to step-size boundary
     660    // align to step-size boundary (search for the nearest one)
     661    Int_t dv = 1;
    605662    while ((int)(v.X())%stepx)
    606         v.Set(v.X()+1, v.Y());
    607 
     663    {
     664        v.Set(v.X()+dv, v.Y());
     665        dv = -TMath::Sign(TMath::Abs(dv)+1, dv);
     666    }
     667
     668    dv = 1;
    608669    while ((int)(v.Y())%stepy)
    609         v.Set(v.X(), v.Y()+1);
     670    {
     671        v.Set(v.X(), v.Y()+dv);
     672        dv = -TMath::Sign(TMath::Abs(dv)+1, dv);
     673    }
    610674
    611675    // draw...
    … …  
    615679}
    616680
     681// --------------------------------------------------------------------------
     682//
     683// Get a rotation matrix which aligns the pointing position
     684// to the center of the x,y plain
     685//
     686TRotation MAstroCatalog::AlignCoordinates(const TVector3 &v) const
     687{
     688    TRotation trans;
     689    trans.RotateZ(-v.Phi());
     690    trans.RotateY(-v.Theta());
     691    trans.RotateZ(-TMath::Pi()/2);
     692    return trans;
     693}
     694
     695// --------------------------------------------------------------------------
     696//
     697// Return the rotation matrix which converts either sky or
     698// local coordinates to coordinates which pole is the current
     699// pointing direction.
     700//
     701TRotation MAstroCatalog::GetGrid(Bool_t local)
     702{
     703    const Bool_t enable = fTime && fObservatory;
     704
     705    if (!local)
     706    {
     707        const TRotation trans(AlignCoordinates(fRaDec));
     708
     709        DrawGrid(fRaDec, trans, 1);
     710
     711        if (enable)
     712        {
     713            const MAstroSky2Local rot(*fTime, *fObservatory);
     714            DrawGrid(rot*fRaDec, trans*rot.Inverse(), 2);
     715        }
     716
     717        return trans;
     718    }
     719
     720    if (local && enable)
     721    {
     722        const MAstroSky2Local rot(*fTime, *fObservatory);
     723
     724        const TRotation trans(AlignCoordinates(rot*fRaDec));
     725
     726        DrawGrid(fRaDec,     trans*rot, 1);
     727        DrawGrid(rot*fRaDec, trans,     2);
     728
     729        return trans*rot;
     730    }
     731
     732    return TRotation();
     733}
     734
     735// --------------------------------------------------------------------------
     736TString MAstroCatalog::GetPadTitle() const
     737{
     738    const Double_t ra  = fRaDec.Phi()*TMath::RadToDeg();
     739    const Double_t dec = (TMath::Pi()/2-fRaDec.Theta())*TMath::RadToDeg();
     740
     741    TString txt;
     742    txt += Form("\\alpha=%.2fh ",      fmod(ra/15+48, 24));
     743    txt += Form("\\delta=%.1f\\circ ", fmod(dec+270,180)-90);
     744    txt += Form("/ FOV=%.1f\\circ",    fRadiusFOV);
     745
     746    if (!fTime || !fObservatory)
     747        return txt;
     748
     749    const MAstroSky2Local rot(*fTime, *fObservatory);
     750    const TVector3 loc = rot*fRaDec;
     751
     752    const Double_t rho = rot.RotationAngle(fRaDec.Phi(), TMath::Pi()/2-fRaDec.Theta());
     753
     754    const Double_t zd = TMath::RadToDeg()*loc.Theta();
     755    const Double_t az = TMath::RadToDeg()*loc.Phi();
     756
     757    txt.Prepend("#splitline{");
     758
     759    txt += Form("  \\theta=%.1fh ",    fmod(zd+270,180)-90);
     760    txt += Form("\\phi=%.1f\\circ ",   fmod(az+720, 360));
     761    txt += Form(" / \\rho=%.1f\\circ", rho*TMath::RadToDeg());
     762    txt += "}{<";
     763    txt += fTime->GetSqlDateTime();
     764    txt += ">}";
     765    return txt;
     766}
     767
     768// --------------------------------------------------------------------------
    617769void MAstroCatalog::AddPrimitives(Option_t *o)
    618770{
    619     // Precalc Sin/Cos...
    620     TRotation trans;
    621     trans.RotateZ(-fRaDec.Phi());
    622     trans.Rotate(TMath::Pi()/2-fRaDec.Theta(), TVector3(0, 1, 0));
    623 
    624     if (fTime && fObservatory)
    625     {
    626         const TRotation rot(MAstroSky2Local(*fTime, *fObservatory));
    627         const TVector3 zdaz0 = rot*fRaDec;
    628         const TVector2 zdaz(zdaz0.Phi(), zdaz0.Theta());
    629         DrawNet(zdaz, trans*rot.Inverse(), 2);
    630     }
    631 
    632     const TVector2 radec(fRaDec.Phi(), fRaDec.Theta());
    633     DrawNet(radec, trans, 1);
     771    const Bool_t local = TString(o).Contains("local", TString::kIgnoreCase);
     772
     773    cout << "Opt: " << o << endl;
     774
     775    const TRotation rot(GetGrid(local));
    634776
    635777    TIter Next(&fList);
    … …  
    639781        // FIXME: Check Magnitude!
    640782        TVector2 s(v->Phi(), v->Theta());
    641         if (Convert(trans, s)==kTRUE)
     783        if (Convert(rot, s)==kTRUE)
    642784            DrawStar(s.X(), TMath::Pi()/2-s.Y(), *v, kFALSE);
    643785    }
    644 }
    645 
     786
     787    TPaveText *pv = new TPaveText(0.01, 0.90, 0.63, 0.99, "brNDC");
     788    pv->AddText(GetPadTitle());
     789    AddMap(pv);
     790
     791    TMarker *mk=new TMarker(0, 0, kMultiply);
     792    mk->SetMarkerColor(kBlack);
     793    mk->SetMarkerSize(1.5);
     794    AddMap(mk);
     795}
     796
     797// --------------------------------------------------------------------------
    646798void MAstroCatalog::SetRangePad()
    647799{
    … …  
    658810}
    659811
     812// --------------------------------------------------------------------------
    660813void MAstroCatalog::DrawPrimitives(Option_t *o)
    661814{
    … …  
    680833}
    681834
     835// --------------------------------------------------------------------------
    682836void MAstroCatalog::DrawMap()
    683837{
    … …  
    688842}
    689843
     844// --------------------------------------------------------------------------
    690845void MAstroCatalog::Draw(Option_t *o)
    691846{
    … …  
    770925}
    771926
     927// --------------------------------------------------------------------------
    772928void MAstroCatalog::ExecuteEventKbd(Int_t keycode, Int_t keysym)
    773929{
    … …  
    821977}
    822978
     979// --------------------------------------------------------------------------
    823980Int_t MAstroCatalog::DistancetoPrimitive(Int_t px, Int_t py)
    824981{
    … …  
    8431000}
    8441001
     1002// --------------------------------------------------------------------------
    8451003void MAstroCatalog::ShowToolTip(Int_t px, Int_t py, const char *txt)
    8461004{
    … …  
    8581016    fToolTip->Show(x+4, y+4);
    8591017}
    860 

  • trunk/MagicSoft/Mars/mastro/MAstroCatalog.h

    r3537 r3568  
    7474{
    7575private:
    76     Double_t fLimMag;    // [1]   Limiting Magnitude
    77     Double_t fRadiusFOV; // [deg] Radius of Field of View
    78 
    79     TGToolTip *fToolTip; //!
     76    Double_t   fLimMag;    // [1]   Limiting Magnitude
     77    Double_t   fRadiusFOV; // [deg] Radius of Field of View
     78
     79    TExMap     fMapG;      //! A map with all gui primitives and tooltips
     80    TGToolTip *fToolTip;   //! The tooltip currently displayed
    8081
    8182    void ShowToolTip(Int_t px, Int_t py, const char *txt);
    … …  
    9293//#endif
    9394
    94 protected:
    95     enum {
    96         kHasChanged  = BIT(15),
    97         kGuiActive   = BIT(16),
    98         kPlainScreen = BIT(17)
    99     };
    100 
    101     TExMap   fMapG;
    102     TList    fList;      // List of stars loaded
    103     MVector3 fRaDec;     // pointing position
    104 
    105     MObservatory *fObservatory; // Possible obervatora location
    106     MTime        *fTime;        // Possible observation time
    107 
    108     virtual Int_t  ConvertToPad(const TVector3 &w, TVector2 &v);
    109     virtual Int_t  Convert(const TRotation &rot, TVector2 &v);
    110     virtual Bool_t DrawLine(const TVector2 &v0, Double_t dx, Double_t dy, const TRotation &rot, Int_t type);
    111     virtual void   AddPrimitives(Option_t *o);
    112     virtual void   DrawPrimitives(Option_t *o);
    113     virtual void   SetRangePad();
    114     void  Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type);
    115     void  DrawNet(const TVector2 &v0, const TRotation &rot, Int_t type);
    116     void  DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t t, const char *txt=0);
    117     void  Paint(Option_t *o="");
    118     Int_t DistancetoPrimitive(Int_t px, Int_t py);
    119     void  Update();
    120 
    121     void ExecuteEventKbd(Int_t keycode, Int_t keysym);
    122     void ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2);
    123 
    124     void DrawMap();
    125     void AddMap(void *k, void *v=0) { fMapG.Add(fMapG.GetSize(), (Long_t)k, (Long_t)v); }
    126     void DeleteMap()
     95    virtual Int_t ConvertToPad(const TVector3 &w, TVector2 &v) const;
     96    virtual void  AddPrimitives(Option_t *o);
     97    virtual void  SetRangePad();
     98
     99    Int_t     Convert(const TRotation &rot, TVector2 &v) const;
     100    void      Draw(const TVector2 &v0, const TRotation &rot, TArrayI &dx, TArrayI &dy, Int_t stepx, Int_t stepy, Int_t type);
     101    void      DrawPrimitives(Option_t *o);
     102    Bool_t    DrawLine(const TVector2 &v0, Double_t dx, Double_t dy, const TRotation &rot, Int_t type);
     103    void      DrawGrid(const TVector3 &v0, const TRotation &rot, Int_t type);
     104    TRotation AlignCoordinates(const TVector3 &v) const;
     105    void      Paint(Option_t *o="");
     106    Int_t     DistancetoPrimitive(Int_t px, Int_t py);
     107    void      DrawMap();
     108    void      DeleteMap()
    127109    {
    128110        Long_t key, val;
    … …  
    148130    }
    149131
     132protected:
     133    enum {
     134        kHasChanged  = BIT(15),
     135        kGuiActive   = BIT(16),
     136        kPlainScreen = BIT(17)
     137    };
     138
     139    TList    fList;      // List of stars loaded
     140    MVector3 fRaDec;     // pointing position
     141
     142    MObservatory *fObservatory; // Possible obervatora location
     143    MTime        *fTime;        // Possible observation time
     144
     145    virtual TString GetPadTitle() const;
     146    TRotation GetGrid(Bool_t local);
     147    void      DrawStar(Double_t x, Double_t y, const TVector3 &v, Bool_t t, const char *txt=0);
     148    void      Update(Bool_t upd=kFALSE);
     149
     150    void ExecuteEventKbd(Int_t keycode, Int_t keysym);
     151    void ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2);
     152
     153    void AddMap(TObject *k, void *v=0)
     154    {
     155        k->SetBit(kCanDelete);
     156        k->SetBit(kCannotPick);
     157        fMapG.Add(fMapG.GetSize(), (Long_t)k, (Long_t)v);
     158    }
     159
    150160public:
    151161    MAstroCatalog();
    … …  
    161171        if (deg>max)
    162172            deg=max;
    163         if (deg<=0)
     173        if (deg<1)
    164174            deg=1;
    165175
    … …  
    189199
    190200    void Draw(Option_t *o="");
     201    void SetDrawOption(Option_t *option="")
     202    {
     203        TObject::SetDrawOption(option);
     204        Update(kTRUE);
     205    } //*MENU*
    191206
    192207    virtual void EventInfo(Int_t event, Int_t px, Int_t py, TObject *selected=0);

  • trunk/MagicSoft/Mars/mastro/MAstroSky2Local.cc

    r3537 r3568  
    6262#include "MAstroSky2Local.h"
    6363
     64#include "MAstro.h"
    6465#include "MTime.h"
    6566#include "MObservatory.h"
    … …  
    128129// seen with an Alt/Az telescope.
    129130//
     131// For more information see MAstro::RotationAngle
     132//
    130133Double_t MAstroSky2Local::RotationAngle(Double_t ra, Double_t dec) const
    131134{
    132     TVector3 loc;
    133     loc.SetMagThetaPhi(1, TMath::Pi()/2-dec, ra);
    134     loc *= *this;
     135    TVector3 v;
     136    v.SetMagThetaPhi(1, TMath::Pi()/2-dec, ra);
     137    v *= *this;
    135138
    136     TRotation rot;
    137     rot.RotateZ(-loc.Phi());
    138     rot.RotateY(-loc.Theta());
    139 
    140     TVector3 v(1, 0, 0);
    141     v *= *this;
    142     v *= rot;
    143 
    144     return TMath::ATan2(v(1), v(0));
     139    return MAstro::RotationAngle(ZZ(), XZ(), v.Theta(), v.Phi());
    145140}

  • trunk/MagicSoft/Mars/mastro/MObservatory.cc

    r3498 r3568  
    101101{
    102102    *fLog << all;
    103     *fLog << fObservatoryName << endl;
    104     *fLog << "Latitude " << (fLatitude > 0 ? (fLatitude*kRad2Deg) : -(fLatitude*kRad2Deg)) << " deg " << (fLatitude > 0 ? "W" : "E") << endl;
    105     *fLog << "Longitude " << (fLongitude > 0 ? (fLongitude*kRad2Deg) : -(fLongitude*kRad2Deg)) <<" deg " << (fLongitude < 0 ? "N" : "S") << endl;
    106     *fLog << "Height " << fHeight << "m" << endl;
     103    *fLog << underline << fObservatoryName << ":" << endl;
     104    *fLog << "Latitude:  " << TMath::Abs(fLatitude*kRad2Deg)  << " deg " << (fLatitude > 0 ? "W" : "E") << endl;
     105    *fLog << "Longitude: " << TMath::Abs(fLongitude*kRad2Deg) << " deg " << (fLongitude < 0 ? "N" : "S") << endl;
     106    *fLog << "Height:    " << fHeight << "m" << endl;
    107107}
    108108
    109109// --------------------------------------------------------------------------
    110110//
    111 // RotationAngle
    112 //
    113 // calculates the angle for the rotation of the sky image in the camera;
    114 // this angle is a function of the local coordinates
     111// Get the corresponding rotation angle of the sky coordinate system
     112// seen with an Alt/Az telescope.
    115113//
    116 //  theta [rad]: polar angle/zenith distance
    117 //  phi   [rad]: rotation angle/azimuth
    118 //
    119 // Return sin/cos component of angle
    120 //
    121 // calculate rotation angle alpha of sky image in camera
    122 // (see TDAS 00-11, eqs. (18) and (20))
     114// For more information see MAstro::RotationAngle
    123115//
    124116void MObservatory::RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const
    125117{
    126     const Double_t sint = TMath::Sin(theta);
    127     const Double_t cost = TMath::Cos(theta);
    128 
    129     const Double_t sinl = fSinLatitude*sint;
    130     const Double_t cosl = fCosLatitude*cost;
    131 
    132     const Double_t sinp = TMath::Sin(phi);
    133     const Double_t cosp = TMath::Cos(phi);
    134 
    135     const Double_t v1 = sint*sinp;
    136     const Double_t v2 = cosl - sinl*cosp;
    137 
    138     const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2);
    139 
    140     sin =  fCosLatitude*sinp / denom;
    141     cos = (sinl - cosl*cosp) / denom;
     118    MAstro::RotationAngle(fSinLatitude, fCosLatitude, theta, phi, sin, cos);
    142119}
    143120
    144121// --------------------------------------------------------------------------
    145122//
    146 // RotationAngle
    147 //
    148 // calculates the angle for the rotation of the sky image in the camera;
    149 // this angle is a function of the local coordinates
     123// Get the corresponding rotation angle of the sky coordinate system
     124// seen with an Alt/Az telescope.
    150125//
    151 //  theta [rad]: polar angle/zenith distance
    152 //  phi   [rad]: rotation angle/azimuth
    153 //
    154 // Return RotationAngle in rad
    155 //
    156 // calculate rotation angle alpha of sky image in camera
    157 // (see TDAS 00-11, eqs. (18) and (20))
     126// For more information see MAstro::RotationAngle
    158127//
    159128Double_t MObservatory::RotationAngle(Double_t theta, Double_t phi) const
    160129{
    161     const Double_t sint = TMath::Sin(theta);
    162     const Double_t cost = TMath::Cos(theta);
    163 
    164     const Double_t sinp = TMath::Sin(phi);
    165     const Double_t cosp = TMath::Cos(phi);
    166 
    167     const Double_t v1 = sint*sinp;
    168     const Double_t v2 = fCosLatitude*cost - fSinLatitude*sint*cosp;
    169 
    170     const Double_t denom = TMath::Sqrt(v1*v1 + v2*v2);
    171 
    172     return TMath::ASin((fCosLatitude*sinp) / denom);
     130    return MAstro::RotationAngle(fSinLatitude, fCosLatitude, theta, phi);
    173131}
    174 
    175 // --------------------------------------------------------------------------
    176 //
    177 // RotationAngle
    178 //
    179 // calculates the angle for the rotation of the sky image in the camera;
    180 // this angle is a function of the sky coordinates, the observatory
    181 // location and the time
    182 //
    183 //  ra  [rad]: Right ascension
    184 //  dec [rad]: Declination
    185 //
    186 // Return RotationAngle in rad
    187 //
    188 Double_t MObservatory::RotationAngle(Double_t ra, Double_t dec, const MTime &t) const
    189 {
    190     const Double_t alpha = t.GetGmst() + GetElong();
    191 
    192     TVector3 v;
    193     v.SetMagThetaPhi(1, TMath::Pi()/2-dec, alpha-ra);
    194     v.RotateY(GetPhi()-TMath::Pi()/2);
    195 
    196     return RotationAngle(v.Theta(), v.Phi());
    197 }

  • trunk/MagicSoft/Mars/mastro/MObservatory.h

    r3497 r3568  
    5858    void RotationAngle(Double_t theta, Double_t phi, Double_t &sin, Double_t &cos) const;
    5959    Double_t RotationAngle(Double_t theta, Double_t phi) const;
    60     Double_t RotationAngle(Double_t ra, Double_t dec, const MTime &t) const;
    6160
    6261    LocationName_t GetObservatoryKey() const { return fObservatoryKey; }

  • trunk/MagicSoft/Mars/mbase/MEvtLoop.cc

    r2728 r3568  
    496496    *fLog << all << "Ready!" << endl << endl;
    497497
    498     *fLog << dec << endl << "CPU  - "
    499         << "Time: " << clock.CpuTime() << "s"
    500         << " for " << numcnts << " Events"
    501         << " --> " << numcnts/clock.CpuTime() << " Events/s"
    502         << endl;
    503     *fLog << "Real - "
    504         << "Time: " << clock.RealTime() << "s"
    505         << " for " << numcnts << " Events"
    506         << " --> " << numcnts/clock.RealTime() << " Events/s"
    507         << endl << endl;
     498    *fLog << dec << endl << "CPU  - Time: ";
     499    *fLog << clock.CpuTime() << "s" << " for " << numcnts << " Events";
     500    if (numcnts>0)
     501        *fLog << " --> " << numcnts/clock.CpuTime() << " Events/s";
     502    *fLog << endl << "Real - Time: ";
     503    *fLog << clock.RealTime() << "s" << " for " << numcnts << " Events";
     504    if (numcnts>0)
     505        *fLog << " --> " << numcnts/clock.RealTime() << " Events/s";
     506
     507    *fLog << endl << endl;
    508508
    509509    return rc!=kERROR;

  • trunk/MagicSoft/Mars/mbase/MParList.cc

    r2556 r3568  
    329329//  'name' is the name of the object you are searching for.
    330330//
     331// In words: Find object name and check whether it inherits from classname
     332//
    331333TObject *MParList::FindObject(const char *name, const char *classname) const
    332334{
    … …  
    427429//  "Name;1", "Name;2", ... If a semicolon is detected leading dots
    428430//  are stripped from the object-name (eg. "name;5.")
     431//
     432// In words: Create object of type classname and set its name to objname.
     433//           If an object with objname already exists return it.
    429434//
    430435MParContainer *MParList::FindCreateObj(const char *classname, const char *objname)

  • trunk/MagicSoft/Mars/mfileio/MCT1ReadAscii.cc

    r2781 r3568  
    3838//                                                                         //
    3939/////////////////////////////////////////////////////////////////////////////
    40 
    4140#include "MCT1ReadAscii.h"
    4241
     42#include <errno.h>
    4343#include <fstream>
    4444
    … …  
    132132    const char *expname = gSystem->ExpandPathName(name);
    133133    fIn = new ifstream(expname);
     134
     135    const Bool_t noexist = !(*fIn);
     136    if (noexist)
     137    {
     138        *fLog << err << "Cannot open file " << expname << ": ";
     139        *fLog << strerror(errno) << endl;
     140    }
     141    else
     142        *fLog << inf << "Open file: '" << name << "'" << endl;
     143
    134144    delete [] expname;
    135 
    136     const Bool_t noexist = !(*fIn);
    137 
    138     if (noexist)
    139         *fLog << dbginf << "Cannot open file '" << name << "'" << endl;
    140     else
    141         *fLog << "Open file: '" << name << "'" << endl;
    142145
    143146    //

  • trunk/MagicSoft/Mars/mfileio/MCT1ReadPreProc.cc

    r2781 r3568  
    4747#include "MCT1ReadPreProc.h"
    4848
     49#include <errno.h>
    4950#include <fstream>
    5051
    5152#include <TList.h>
    5253#include <TSystem.h>
     54#include <TRandom3.h>
    5355
    5456#define LINUX
    … …  
    8183#include "MBinning.h"
    8284
    83 #include "TRandom3.h"
    8485#include "MParameters.h"
    8586
    … …  
    590591
    591592    fIn = new ifstream(fname);
     593    if (!*fIn)
     594    {
     595        *fLog << err << "Cannot open file " << fname << ": ";
     596        *fLog << strerror(errno) << endl;
     597        return kFALSE;
     598    }
    592599
    593600    *fLog << inf << "-----------------------------------------------------------------------" << endl;

  • trunk/MagicSoft/Mars/mfileio/MReadRflFile.cc

    r2296 r3568  
    3232#include "MReadRflFile.h"
    3333
     34#include <errno.h>
    3435#include <fstream>
    3536
    … …  
    315316    if (!*fIn)
    316317    {
    317         cout << "Error openeng file " << name << "." << endl;
     318        *fLog << err << "Cannot open file " << name << ": ";
     319        *fLog << strerror(errno) << endl;
    318320        return kFALSE;
    319321    }

  • trunk/MagicSoft/Mars/mfileio/MReadTree.cc

    r3503 r3568  
    549549    if (TestBit(kChainWasChanged))
    550550    {
    551         *fLog << inf << "Scanning chain... " << flush;
     551        *fLog << inf << "Scanning chain " << fChain->GetName() << "... " << flush;
    552552        fNumEntries = (UInt_t)fChain->GetEntries();
    553553        *fLog << fNumEntries << " events found." << endl;

  • trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

    r3557 r3568  
    9595//  - a more clever (and faster) algorithm to fill the histogram, eg by
    9696//    calculating alpha once and fill the two coils around the mean
     97//  - more drawing options...
     98//  - Move Significance() to a more 'general' place and implement
     99//    also different algorithms like (Li/Ma)
     100//  - implement fit for best alpha distribution -- online (Paint)
    97101//
    98102//////////////////////////////////////////////////////////////////////////////
    … …  
    101105#include <TF1.h>
    102106#include <TH2.h>
     107#include <TGraph.h>
    103108#include <TStyle.h>
    104109#include <TCanvas.h>
    105110#include <TPaveText.h>
     111#include <TStopwatch.h>
    106112
    107113#include "MGeomCam.h"
    … …  
    236242// Calculate Significance as
    237243// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
     244//
     245// s: total number of events in signal region
     246// b: number of background events in signal region
    238247//
    239248Double_t MHFalseSource::Significance(Double_t s, Double_t b)
    … …  
    245254}
    246255
     256// --------------------------------------------------------------------------
     257//
     258//  calculates the significance according to Li & Ma
     259//  ApJ 272 (1983) 317
     260//
     261//    s: total number of events in signal region
     262//
     263//    Double_t fGamma;   // Nbg = Non - gamma * Noff
     264//  - the effective number of background events (fNoff), and fGamma :
     265//
     266//    fGamma = b / fNoff;
     267//    fGamma = fdNbg / sqrt(fNoff);
     268//    fGamma = fdNbg*fdNbg / fNbg;
     269//    fNoff  =  b*b  / (fdNbg*fdNbg);
     270//     Double_t fNbg;     // number of background events in signal region
     271//    b = fNOff *fGamma;
     272/*
     273Double_t MHFindSignificance::SignificanceLiMa(Double_t non, Double_t noff, Double_t gamma)
     274{
     275    if (gamma <= 0.0  ||  non <= 0.0  ||  noff <= 0.0)
     276    {
     277        *siglima = 0.0;
     278        return kFALSE;
     279    }
     280
     281    Double_t help1 = TMath::Log( (gamma+1)*s  / (gamma*(s+noff)) );
     282    Double_t help2 = TMath::Log( (gamma+1)*noff / (       s+noff ) );
     283
     284    Double_t siglima = TMath::Sqrt((s*help1+noff*help2)*2);
     285
     286    return non<gamma*noff ? -siglima : siglima;
     287}
     288*/
    247289// --------------------------------------------------------------------------
    248290//
    … …  
    276318
    277319        MBinning b;
    278         b.SetEdges(100, -r, r);
     320        b.SetEdges(20, -r, r);
    279321        SetBinning(&fHist, &b, &b, &binsa);
    280322    }
    … …  
    311353    Double_t rho = 0;
    312354    if (fTime && fObservatory && fPointPos)
    313     {
    314         const Double_t ra  = fPointPos->GetRa();
    315         const Double_t dec = fPointPos->GetDec();
    316 
    317         rho = MAstroSky2Local(*fTime, *fObservatory).RotationAngle(ra, dec);
    318     }
     355        rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     356    //if (fPointPos)
     357    //    rho = fPointPos->RotationAngle(*fObservatory);
    319358
    320359    MSrcPosCam src;
    … …  
    554593    h4->Draw("colz");
    555594    h4->SetBit(kCanDelete);
    556 
    557595}
    558596
    … …  
    587625    gPad->cd();
    588626}
    589 
     627/*
    590628Double_t fcn(Double_t *arg, Double_t *p)
    591629{
    … …  
    595633
    596634    const Double_t f1 = p[0]*TMath::Exp(-0.5*dx*dx);
    597     const Double_t f2 = p[3]*x*x + p[4];
     635    const Double_t f2 = p[3] + p[5]*x*x;
    598636
    599637    return f1 + f2;
    … …  
    602640Double_t FcnI1(Double_t x, Double_t *p)
    603641{
    604     return (p[3]*x*x/3+p[4])*x;
     642    return (p[5]*x*x/3+p[3])*x;
    605643}
    606644Double_t FcnI2(Double_t x, Double_t *p)
    … …  
    615653    return f2;
    616654}
    617 
    618 
    619 void MHFalseSource::FitSignificance(Float_t sigmax, Float_t bgmin, Float_t bgmax)
    620 {
    621     TH1D h0a("A",     "Parameter A",       50,  0, 4000);
    622     TH1D h0b("a",     "Parameter a",       50,  0, 4000);
    623     TH1D h1("mu",     "Parameter \\mu",    50, -1, 1);
    624     TH1D h2("sigma",  "Parameter \\sigma", 50,  0, 90);
    625     TH1D h3("b",      "Parameter b",       50, -0.1, 0.1);
    626     TH1D h4a("chisq1", "\\chi^{2} (red, green) / significance (black)",       50,  0, 35);
    627     TH1D h5a("prob1",  "Fit probability",  50,  0, 1.1);
    628     TH1D h4b("chisq2", "\\chi^{2} (red, green) / significance (black)",       50,  0, 35);
    629     TH1D h5b("prob",  "Fit probability",  50,  0, 1.1);
    630     TH1D h6("Signif",  "Significance",     50,  -20, 20);
    631     h0a.SetDirectory(NULL);
    632     h0b.SetDirectory(NULL);
    633     h1.SetDirectory(NULL);
    634     h2.SetDirectory(NULL);
    635     h3.SetDirectory(NULL);
    636     h4a.SetDirectory(NULL);
    637     h5b.SetDirectory(NULL);
    638     h4a.SetDirectory(NULL);
    639     h5b.SetDirectory(NULL);
    640     h6.SetDirectory(NULL);
     655*/
     656/*
     657class MHSignificance : public MH
     658{
     659private:
     660    TH1D fHist;
     661
     662    MParameterD *fParam;
     663
     664public:
     665    MHSignificance() : fParam(0)
     666    {
     667        fHist.SetName("Alpha");
     668        fHist.SetTitle("Distribution of \\alpha");
     669        fHist.SetXTitle("\\alpha [\\circ]");
     670        fHist.SetYTitle("Counts");
     671    }
     672    Int_t SetupFill(MParList *p)
     673    {
     674        fHist.Reset();
     675
     676        fParam = (MParameterD*)p->FindCreateObj("Significance", "MParameterD");
     677        if (fParam)
     678            return kFALSE;
     679
     680        return kTRUE;
     681    }
     682    Int_t Process(MParContainer *p, Double_t w=1)
     683    {
     684        MHillasSrc *hil = dynamic_cast<MHillasSrc*>(p);
     685        if (!hil)
     686        {
     687            *fLog << err << dbginf << "Got no MHillasSrc as argument of Fill()..." << endl;
     688            return kFALSE;
     689        }
     690
     691        fHist->Fill(hil->GetAlpha(), w);
     692
     693        return kTRUE;
     694    }
     695    Int_t Finalize()
     696    {
     697        if (fHist.GetEntries()==0)
     698        {
     699            *fLog << err << "Histogram empty." << endl;
     700            return kFALSE;
     701        }
     702
     703        Float_t sigmax=15;
     704        Float_t bgmin =45;
     705        Float_t bgmax =80;
     706
     707        fHist.SetNameTitle("Significance",
     708                           Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
     709                                sigmax, bgmin, bgmax));
     710
     711        // Implementing the function yourself is not faster at all!
     712        TF1 func("gaus(0) + pol2(3)", fcn, 0, 90, 6);
     713        TArrayD maxpar(func.GetNpar());
     714
     715        func.FixParameter(1, 0);
     716        func.FixParameter(4, 0);
     717        func.SetParLimits(3, -1, 1);
     718
     719        const Double_t alpha0 = fHist.GetBinContent(1);
     720
     721        // First fit a polynom in the off region
     722        func.FixParameter(0, 0);
     723        func.FixParameter(2, 1);
     724        func.ReleaseParameter(3);
     725        func.ReleaseParameter(5);
     726
     727        h->Fit(&func, "N0Q", "", bgmin, bgmax);
     728
     729        // Now fit a gaus in the on region on top of the polynom
     730        func.SetParameter(0, alpha0-func.GetParameter(4));
     731        func.SetParameter(2, sigmax*0.75);
     732
     733        func.ReleaseParameter(0);
     734        func.ReleaseParameter(2);
     735        func.FixParameter(3, func.GetParameter(3));
     736        func.FixParameter(5, func.GetParameter(5));
     737
     738        func.SetParLimits(2, 0, 80);
     739        h->Fit(&func, "N0Q", "", 0, sigmax);
     740
     741        TArrayD p(func.GetNpar(), func.GetParameters());
     742
     743        Double_t sig=0;
     744
     745        const Int_t n = hist->GetBin(ix+1, iy+1);
     746        if (!(func.GetParameter(0)>alpha0*2 ||
     747              func.GetParameter(2)<2.5      ||
     748              func.GetParameter(2)>70))
     749        {
     750            // Implementing the integral as analytical function
     751            // gives the same result in the order of 10e-5
     752            // and it is not faster at all...
     753            const Double_t s = func.Integral(0, 15);
     754
     755            func.SetParameter(0, 0);
     756            func.SetParameter(2, 1);
     757
     758            const Double_t b = func.Integral(0, 15);
     759
     760            sig = Significance(s, b);
     761        }
     762
     763        fParam->SetValue(sig);
     764        fParam->SetReadyToSave();
     765        return kTRUE;
     766    }
     767    ClassDef(MHSignificance, 0)
     768};
     769*/
     770
     771
     772// --------------------------------------------------------------------------
     773//
     774// This is a preliminary implementation of a alpha-fit procedure for
     775// all possible source positions. It will be moved into its own
     776// more powerfull class soon.
     777//
     778// The fit function is "gaus(0)+pol2(3)" which is equivalent to:
     779//   [0]*exp(-0.5*((x-[1])/[2])^2) + [3] + [4]*x + [5]*x^2
     780// or
     781//   A*exp(-0.5*((x-mu)/sigma)^2) + a + b*x + c*x^2
     782//
     783// Parameter [1] is fixed to 0 while the alpha peak should be
     784// symmetric around alpha=0.
     785//
     786// Parameter [4] is fixed to 0 because the first derivative at
     787// alpha=0 should be 0, too.
     788//
     789// In a first step the background is fitted between bgmin and bgmax,
     790// while the parameters [0]=0 and [2]=1 are fixed.
     791//
     792// In a second step the signal region (alpha<sigmax) is fittet using
     793// the whole function with parameters [1], [3], [4] and [5] fixed.
     794//
     795// The number of excess and background events are calculated as
     796//   s = int(0, sigint, gaus(0)+pol2(3))
     797//   b = int(0, sigint,         pol2(3))
     798//
     799// The Significance is calculated using the Significance() member
     800// function.
     801//
     802void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
     803{
     804    TH1D h0a("A",          "", 50,   0, 4000);
     805    TH1D h4a("chisq1",     "", 50,   0,   35);
     806    //TH1D h5a("prob1",      "", 50,   0,  1.1);
     807    TH1D h6("signifcance", "", 50, -20,   20);
     808
     809    TH1D h1("mu",    "Parameter \\mu",    50,   -1,    1);
     810    TH1D h2("sigma", "Parameter \\sigma", 50,    0,   90);
     811    TH1D h3("b",     "Parameter b",       50, -0.1,  0.1);
     812
     813    TH1D h0b("a",         "Parameter a (red), A (blue)", 50, 0, 4000);
     814    TH1D h4b("\\chi^{2}", "\\chi^{2} (red, green) / significance (black)", 50, 0, 35);
     815    //TH1D h5b("prob",      "Fit probability: Bg(red), F(blue)", 50, 0, 1.1);
    641816
    642817    h0a.SetLineColor(kBlue);
    643818    h4a.SetLineColor(kBlue);
    644     h5a.SetLineColor(kBlue);
     819    //h5a.SetLineColor(kBlue);
    645820    h0b.SetLineColor(kRed);
    646821    h4b.SetLineColor(kRed);
    647     h5b.SetLineColor(kRed);
    648 
    649     TH1 *hist = fHist.Project3D("xy_fit");
     822    //h5b.SetLineColor(kRed);
     823
     824    TH1 *hist  = fHist.Project3D("xy_fit");
     825    hist->SetDirectory(0);
     826    TH1 *hists = fHist.Project3D("xy_fit");
     827    hists->SetDirectory(0);
     828    TH1 *histb = fHist.Project3D("xy_fit");
     829    histb->SetDirectory(0);
    650830    hist->Reset();
     831    hists->Reset();
     832    histb->Reset();
    651833    hist->SetNameTitle("Significance",
    652834                       Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
    653835                            sigmax, bgmin, bgmax));
     836    hists->SetNameTitle("Excess",     Form("Number of excess events for \\alpha<%.0f\\circ", sigint));
     837    histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint));
    654838    hist->SetXTitle(fHist.GetXaxis()->GetTitle());
     839    hists->SetXTitle(fHist.GetXaxis()->GetTitle());
     840    histb->SetXTitle(fHist.GetXaxis()->GetTitle());
    655841    hist->SetYTitle(fHist.GetXaxis()->GetTitle());
     842    hists->SetYTitle(fHist.GetXaxis()->GetTitle());
     843    histb->SetYTitle(fHist.GetXaxis()->GetTitle());
    656844
    657845    //                      xmin, xmax, npar
    658     TF1 func("MyFunc", fcn, 0,    90,   5);
    659     TArrayD par(5);
    660 
    661     func.SetParName(0, "A");
    662     func.SetParName(1, "mu");
    663     func.SetParName(2, "sigma");
    664     func.SetParName(3, "a");
    665     func.SetParName(4, "b");
    666 
     846    //TF1 func("MyFunc", fcn, 0, 90, 6);
     847    // Implementing the function yourself is only about 5% faster
     848    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
     849    TArrayD maxpar(func.GetNpar());
     850
     851    /*
     852     func.SetParName(0, "A");
     853     func.SetParName(1, "mu");
     854     func.SetParName(2, "sigma");
     855     func.SetParName(3, "a");
     856     func.SetParName(4, "b");
     857     func.SetParName(5, "c");
     858    */
     859
     860    func.FixParameter(1, 0);
     861    func.FixParameter(4, 0);
     862    func.SetParLimits(2, 0, 90);
    667863    func.SetParLimits(3, -1, 1);
    668864
    … …  
    676872    Int_t maxy=0;
    677873
     874    TStopwatch clk;
     875    clk.Start();
     876
     877    *fLog << inf;
     878    *fLog << "Signal fit:     alpha < " << sigmax << endl;
     879    *fLog << "Integration:    alpha < " << sigint << endl;
     880    *fLog << "Background fit: " << bgmin << " < alpha < " << bgmax << endl;
     881    *fLog << "Polynom order:  " << (int)polynom << endl;
     882    *fLog << "Fitting False Source Plot..." << flush;
     883
    678884    TH1 *h=0;
    679     for (int ix=0; ix<nx; ix++)
    680         for (int iy=0; iy<ny; iy++)
     885    for (int ix=1; ix<=nx; ix++)
     886        for (int iy=1; iy<=ny; iy++)
    681887        {
    682             h = fHist.ProjectionZ("AlphaFit", ix+1, ix+1, iy+1, iy+1);
     888            h = fHist.ProjectionZ("AlphaFit", ix, ix, iy, iy);
    683889
    684890            const Double_t alpha0 = h->GetBinContent(1);
     891            const Double_t alphaw = h->GetXaxis()->GetBinWidth(1);
    685892
    686893            // Check for the regios which is not filled...
    … …  
    693900            // First fit a polynom in the off region
    694901            func.FixParameter(0, 0);
    695             func.FixParameter(1, 0);
    696902            func.FixParameter(2, 1);
    697903            func.ReleaseParameter(3);
    698             func.ReleaseParameter(4);
     904
     905            for (int i=5; i<func.GetNpar(); i++)
     906                func.ReleaseParameter(i);
    699907
    700908            h->Fit(&func, "N0Q", "", bgmin, bgmax);
    … …  
    702910
    703911            h4a.Fill(func.GetChisquare());
    704             h5a.Fill(func.GetProb());
    705 
    706             // Now fit a gaus in the on region on top of the polynom
    707             func.SetParameter(0, alpha0-func.GetParameter(4));
    708             func.SetParameter(2, sigmax*0.75);
     912            //h5a.Fill(func.GetProb());
    709913
    710914            //func.SetParLimits(0, 0.5*h->GetBinContent(1), 1.5*h->GetBinContent(1));
    … …  
    715919            func.ReleaseParameter(2);
    716920            func.FixParameter(3, func.GetParameter(3));
    717             func.FixParameter(4, func.GetParameter(4));
    718 
    719             func.SetParLimits(2, 0, 80);
     921            for (int i=5; i<func.GetNpar(); i++)
     922                func.FixParameter(i, func.GetParameter(i));
     923
     924            // Do not allow signals smaller than the background
     925            const Double_t A  = alpha0-func.GetParameter(3);
     926            const Double_t dA = TMath::Abs(A);
     927            func.SetParLimits(0, -dA*4, dA*4);
     928            func.SetParLimits(2, 0, 90);
     929
     930            // Now fit a gaus in the on region on top of the polynom
     931            func.SetParameter(0, A);
     932            func.SetParameter(2, sigmax*0.75);
     933
    720934            h->Fit(&func, "N0Q", "", 0, sigmax);
    721935            //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
    722936
     937            TArrayD p(func.GetNpar(), func.GetParameters());
     938
    723939            // Fill results into some histograms
    724             h0a.Fill(func.GetParameter(0));
    725             h0b.Fill(func.GetParameter(4));
    726             h1.Fill(func.GetParameter(1));
    727             h2.Fill(func.GetParameter(2));
    728             h3.Fill(func.GetParameter(3));
     940            h0a.Fill(p[0]);
     941            h0b.Fill(p[3]);
     942            h1.Fill(p[1]);
     943            h2.Fill(p[2]);
     944            if (polynom>1)
     945                h3.Fill(p[5]);
    729946            h4b.Fill(func.GetChisquare());
    730             h5b.Fill(func.GetProb());
    731 
    732             Double_t sig=0;
    733 
    734             const Int_t n = hist->GetBin(ix+1, iy+1);
    735             if (!(func.GetParameter(0)>alpha0*2 ||
    736                   func.GetParameter(2)<2.5      ||
    737                   func.GetParameter(2)>70
    738                 /*func.GetProb()<0.005 ||*/))
     947            //h5b.Fill(func.GetProb());
     948
     949            // Implementing the integral as analytical function
     950            // gives the same result in the order of 10e-5
     951            // and it is not faster at all...
     952            //const Bool_t ok = /*p[0]>=0 && /*p[0]<alpha0*2 &&*/ p[2]>1.75;// && p[2]<88.5;
     953            /*
     954            if (p[0]<-fabs(alpha0-p[3])*7 && p[2]<alphaw/3)
    739955            {
    740                 const Double_t b = FcnI1(15, func.GetParameters());
    741                 const Double_t s = FcnI2(15, func.GetParameters());
    742 
    743                 sig = Significance(s+b, b);
     956                func.SetParameter(0, alpha0-p[3]);
     957                cout << p[2] << " " << p[0] << " " << alpha0-p[3] << endl;
    744958            }
    745 
    746             hist->SetBinContent(n, sig==0?1e-15:sig);
     959            */
     960
     961            const Double_t s = func.Integral(0, sigint);
     962
     963            func.SetParameter(0, 0);
     964            func.SetParameter(2, 1);
     965
     966            const Double_t b   = func.Integral(0, sigint);
     967            const Double_t sig = Significance(s, b);
     968
     969            const Int_t n = hist->GetBin(ix, iy);
     970            hists->SetBinContent(n, s-b);
     971            histb->SetBinContent(n, b);
     972
     973            hist->SetBinContent(n, sig);
    747974            if (sig!=0)
    748975                h6.Fill(sig);
    … …  
    751978            {
    752979                maxs = sig;
    753                 maxx = ix+1;
    754                 maxy = iy+1;
    755                 for(int i=0; i<par.GetSize(); i++)
    756                     par[i] = func.GetParameter(i);
     980                maxx = ix;
     981                maxy = iy;
     982                maxpar = p;
    757983            }
    758984        }
    759985
     986    *fLog << inf << "done." << endl;
     987
    760988    h0a.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    761989    h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
    762990
     991    //hists->SetMinimum(GetMinimumGT(*hists));
     992    histb->SetMinimum(GetMinimumGT(*histb));
     993
     994    clk.Stop();
     995    clk.Print();
     996
    763997    TCanvas *c=new TCanvas;
    764998
    … …  
    7661000    c->cd(1);
    7671001    gPad->SetBorderMode(0);
    768     gPad->Divide(1,2, 0, 0);
     1002    hists->Draw("colz");
     1003    hists->SetBit(kCanDelete);
     1004    c->cd(2);
     1005    gPad->SetBorderMode(0);
     1006    hist->Draw("colz");
     1007    hist->SetBit(kCanDelete);
     1008    c->cd(3);
     1009    gPad->SetBorderMode(0);
     1010    histb->Draw("colz");
     1011    histb->SetBit(kCanDelete);
     1012    c->cd(4);
     1013    gPad->Divide(1,3, 0, 0);
    7691014    TVirtualPad *p=gPad;
    7701015    p->SetBorderMode(0);
    … …  
    7761021    gPad->SetBorderMode(0);
    7771022    h3.DrawCopy();
    778     c->cd(2);
    779     gPad->SetBorderMode(0);
    780     hist->Draw("colz");
    781     hist->SetDirectory(NULL);
    782     hist->SetBit(kCanDelete);
    783     c->cd(3);
     1023    p->cd(3);
     1024    gPad->SetBorderMode(0);
     1025    h2.DrawCopy();
     1026    c->cd(6);
     1027    gPad->Divide(1,2, 0, 0);
     1028    TVirtualPad *q=gPad;
     1029    q->SetBorderMode(0);
     1030    q->cd(1);
     1031    gPad->SetBorderMode(0);
    7841032    gPad->SetBorderMode(0);
    7851033    h4b.DrawCopy();
    7861034    h4a.DrawCopy("same");
    7871035    h6.DrawCopy("same");
    788     c->cd(4);
    789     gPad->SetBorderMode(0);
    790     h2.DrawCopy();
    791     c->cd(6);
    792     gPad->SetBorderMode(0);
    793     h5b.DrawCopy();
    794     h5a.DrawCopy("same");
    795 
     1036    q->cd(2);
     1037    gPad->SetBorderMode(0);
     1038    //h5b.DrawCopy();
     1039    //h5a.DrawCopy("same");
    7961040    c->cd(5);
    7971041    gPad->SetBorderMode(0);
    … …  
    8011045                                 fHist.GetXaxis()->GetBinCenter(maxx),
    8021046                                 fHist.GetYaxis()->GetBinCenter(maxy), maxs);
     1047
    8031048        TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
     1049        result->SetDirectory(NULL);
    8041050        result->SetNameTitle("Result \\alpha", title);
    8051051        result->SetBit(kCanDelete);
    … …  
    8081054        result->Draw();
    8091055
    810         TF1 f1("MyFunc1", fcn, 0, 90, 5);
    811         TF1 f2("MyFunc2", fcn, 0, 90, 5);
    812         f1.SetParameters(par.GetArray());
    813         f2.SetParameters(par.GetArray());
     1056        TF1 f1("", func.GetExpFormula(), 0, 90);
     1057        TF1 f2("", func.GetExpFormula(), 0, 90);
     1058        f1.SetParameters(maxpar.GetArray());
     1059        f2.SetParameters(maxpar.GetArray());
    8141060        f2.FixParameter(0, 0);
    8151061        f2.FixParameter(1, 0);
    … …  
    8241070        leg->SetBorderSize(1);
    8251071        leg->SetTextSize(0.04);
    826         leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
     1072        leg->AddText(0.5, 0.82, Form("A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + pol%d", polynom))->SetTextAlign(22);
     1073        //leg->AddText(0.5, 0.82, "A * exp(-(\\frac{x-\\mu}{\\sigma})^{2}/2) + b*x^{2} + a")->SetTextAlign(22);
    8271074        leg->AddLine(0, 0.65, 0, 0.65);
    828         leg->AddText(0.06, 0.54, Form("A=%.2f", par[0]))->SetTextAlign(11);
    829         leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", par[2]))->SetTextAlign(11);
    830         leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fixed)", par[1]))->SetTextAlign(11);
    831         leg->AddText(0.60, 0.54, Form("a=%.2f", par[3]))->SetTextAlign(11);
    832         leg->AddText(0.60, 0.34, Form("b=%.2f", par[4]))->SetTextAlign(11);
     1075        leg->AddText(0.06, 0.54, Form("A=%.2f", maxpar[0]))->SetTextAlign(11);
     1076        leg->AddText(0.06, 0.34, Form("\\sigma=%.2f", maxpar[2]))->SetTextAlign(11);
     1077        leg->AddText(0.06, 0.14, Form("\\mu=%.2f (fix)", maxpar[1]))->SetTextAlign(11);
     1078        leg->AddText(0.60, 0.54, Form("a=%.2f", maxpar[3]))->SetTextAlign(11);
     1079        leg->AddText(0.60, 0.34, Form("b=%.2f (fix)", maxpar[4]))->SetTextAlign(11);
     1080        if (polynom>1)
     1081            leg->AddText(0.60, 0.14, Form("c=%.2f", maxpar[5]))->SetTextAlign(11);
    8331082        leg->SetBit(kCanDelete);
    8341083        leg->Draw();
    835     }
    836 
    837 }
     1084
     1085        q->cd(2);
     1086
     1087        TGraph *g = new TGraph;
     1088        g->SetPoint(0, 0, 0);
     1089
     1090        Int_t max=0;
     1091        Float_t maxsig=0;
     1092        for (int i=1; i<89; i++)
     1093        {
     1094            const Double_t s = f1.Integral(0, (float)i);
     1095            const Double_t b = f2.Integral(0, (float)i);
     1096
     1097            const Double_t sig = Significance(s, b);
     1098
     1099            g->SetPoint(g->GetN(), i, sig);
     1100
     1101            if (sig>maxsig)
     1102            {
     1103                max = i;
     1104                maxsig = sig;
     1105            }
     1106        }
     1107
     1108        g->SetNameTitle("SigVs\\alpha", "Significance vs \\alpha");
     1109        g->GetHistogram()->SetXTitle("\\alpha_{0} [\\circ]");
     1110        g->GetHistogram()->SetYTitle("Significance");
     1111        g->SetBit(kCanDelete);
     1112        g->Draw("AP");
     1113
     1114        leg = new TPaveText(0.35, 0.10, 0.90, 0.25, "brNDC");
     1115        leg->SetBorderSize(1);
     1116        leg->SetTextSize(0.1);
     1117        leg->AddText(Form("\\sigma_{max}=%.1f at \\alpha_{max}=%d\\circ", maxsig, max));
     1118        leg->SetBit(kCanDelete);
     1119        leg->Draw();
     1120    }
     1121}

  • trunk/MagicSoft/Mars/mhist/MHFalseSource.h

    r3557 r3568  
    5151    TH1 *GetHistByName(const TString name) { return &fHist; }
    5252
    53     void FitSignificance(Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80); //*MENU*
     53    void FitSignificance(Float_t sigint=15, Float_t sigmax=35, Float_t bgmin=40, Float_t bgmax=80, Byte_t polynom=2); //*MENU*
    5454    void FitSignificanceStd() { FitSignificance(); } //*MENU*
    5555

  • trunk/MagicSoft/Mars/mimage/MConcentration.cc

    r3543 r3568  
    116116Int_t MConcentration::Calc(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hillas)
    117117{
    118     const UInt_t npixevt = evt.GetNumPixels();
    119 
    120118    Float_t maxpix[9] = {0,0,0,0,0,0,0,0,0};             // [#phot]
    121119
    122     for (UInt_t i=0; i<npixevt; i++)
     120    TIter Next(evt);
     121    MCerPhotPix *pix = 0;
     122    while ((pix=(MCerPhotPix*)Next()))
    123123    {
    124         const MCerPhotPix &pix = evt[i];
    125 
    126         // skip unused pixels
    127         if (!pix.IsPixelUsed())
    128             continue;
    129 
    130         const Int_t pixid = pix.GetPixId();
    131 
    132         const MGeomPix &gpix = geom[pixid];
    133 
    134         Double_t nphot = pix.GetNumPhotons();
    135 
    136         //
    137         // Now we are working on absolute values of nphot, which
    138         // must take pixel size into account
    139         //
    140         nphot *= geom.GetPixRatio(pixid);
     124        const Int_t    pixid = pix->GetPixId();
     125        const Double_t nphot = pix->GetNumPhotons()* geom.GetPixRatio(pixid);
    141126
    142127        // Get number of photons in the 8 most populated pixels
     128        if (maxpix[0]<=nphot)
     129        {
     130            for(int i=0;i<8;i++)
     131                maxpix[8-i]=maxpix[7-i];
     132            maxpix[0]=nphot;
     133            continue;
     134        }
    143135
    144         if(maxpix[0]<=nphot) {
    145           for(int i=0;i<8;i++)
    146             maxpix[8-i]=maxpix[7-i];
    147           maxpix[0]=nphot;
    148         }
    149         else
    150           for(int i=0;i<8;i++){
    151             if (nphot<maxpix[7-i]){
    152               for(int j=0;j<i-1;j++){ 
    153                 maxpix[7-j]=maxpix[6-j];                 // [#phot]
    154               }
    155               maxpix[8-i]=nphot;
    156               break;
    157             }
    158           }
     136        // Check if the latest value is 'somewhere in between'
     137        for (int i=0; i<8; i++)
     138        {
     139            if (nphot>=maxpix[7-i])
     140                continue;
    159141
     142            for(int j=0;j<i-1;j++)
     143                maxpix[7-j]=maxpix[6-j];                 // [#phot]
     144
     145            maxpix[8-i]=nphot;
     146            break;
     147        }
    160148    }
    161149
    162150    // Compute concentrations from the 8 pixels with higher signal
    163    
    164151    fConc[0]=maxpix[0];
    165     for(int i=1;i<8;i++)
    166       {
    167         fConc[i]=fConc[i-1]+maxpix[i];
    168       }
    169    
    170     for(int i=0;i<8;i++)
    171       {
    172         fConc[i]/=hillas.GetSize();                       // [ratio]
    173       }
     152
     153    // No calculate the integral of the n highest pixels
     154    for(int i=1; i<8; i++)
     155        fConc[i] = fConc[i-1]+maxpix[i];
     156
     157    for(int i=0; i<8; i++)
     158        fConc[i] /= hillas.GetSize();                    // [ratio]
    174159
    175160    SetReadyToSave();
    176161
     162    return 0;
    177163}
    178 
    179 

  • trunk/MagicSoft/Mars/mimage/MHNewImagePar.cc

    r2414 r3568  
    6565    fHistLeakage1.SetYTitle("Counts");
    6666    fHistLeakage1.SetDirectory(NULL);
     67    fHistLeakage1.UseCurrentStyle();
    6768    fHistLeakage1.SetFillStyle(4000);
    68     fHistLeakage1.UseCurrentStyle();
    6969
    7070    fHistLeakage2.SetName("Leakage2");
    … …  
    7373    fHistLeakage2.SetYTitle("Counts");
    7474    fHistLeakage2.SetDirectory(NULL);
     75    fHistLeakage2.UseCurrentStyle();
    7576    fHistLeakage2.SetLineColor(kBlue);
    7677    fHistLeakage2.SetFillStyle(4000);
    77     fHistLeakage2.UseCurrentStyle();
    7878 
    7979    fHistUsedPix.SetName("UsedPix");
    … …  
    8282    fHistUsedPix.SetYTitle("Counts");
    8383    fHistUsedPix.SetDirectory(NULL);
    84     fHistUsedPix.SetLineColor(kGreen);
     84    fHistUsedPix.UseCurrentStyle();
     85    fHistUsedPix.SetLineColor(kBlue);
    8586    fHistUsedPix.SetFillStyle(4000);
    86     fHistUsedPix.UseCurrentStyle();
    8787
    8888    fHistCorePix.SetName("CorePix");
    … …  
    9191    fHistCorePix.SetYTitle("Counts");
    9292    fHistCorePix.SetDirectory(NULL);
    93     fHistCorePix.SetLineColor(kRed);
     93    fHistCorePix.UseCurrentStyle();
     94    fHistCorePix.SetLineColor(kBlack);
    9495    fHistCorePix.SetFillStyle(4000);
    95     fHistCorePix.UseCurrentStyle();
    9696
    9797    fHistConc.SetDirectory(NULL);
    … …  
    105105    fHistConc.SetYTitle("Counts");
    106106    fHistConc1.SetYTitle("Counts");
     107    fHistConc1.UseCurrentStyle();
     108    fHistConc.UseCurrentStyle();
    107109    fHistConc.SetFillStyle(4000);
    108110    fHistConc1.SetFillStyle(4000);
    109111    fHistConc1.SetLineColor(kBlue);
    110112    fHistConc.SetFillStyle(0);
    111     fHistConc1.UseCurrentStyle();
    112     fHistConc.UseCurrentStyle();
    113113
    114114

  • trunk/MagicSoft/Mars/mpointing/MPointingPos.cc

    r2601 r3568  
    1818!   Author(s): Thomas Bretz, 11/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
    1919!
    20 !   Copyright: MAGIC Software Development, 2000-2003
     20!   Copyright: MAGIC Software Development, 2000-2004
    2121!
    2222!
    … …  
    2727// MPointingPos
    2828//
    29 // Store the current pointing position of the telescope...
     29// In this container we store the corrected pointing position of the
     30// telscope. The pointing coordinates are read into MReportDrive together
     31// with its time.
     32//
     33// MPointingPosCalc afterwards calculates corrections and checks for the
     34// cosistency of the coordinates. The result (the real coordinates)
     35// are stored in this container. No further correction should be necessary
     36// using MPointingPos.
     37//
     38// If you need the rotation angle of the starfield in the camera you can
     39// get it from here.
    3040//
    3141/////////////////////////////////////////////////////////////////////////////
    3242#include "MPointingPos.h"
    3343
     44#include "MTime.h"
     45#include "MObservatory.h"
     46#include "MAstroSky2Local.h"
     47
    3448ClassImp(MPointingPos);
    3549
    3650using namespace std;
     51
     52// --------------------------------------------------------------------------
     53//
     54// Get the corresponding rotation angle of the sky coordinate system
     55// seen with an Alt/Az telescope calculated from the stored local
     56// (Zd/Az) coordinates.
     57//
     58// For more information see MAstro::RotationAngle
     59//
     60Double_t MPointingPos::RotationAngle(const MObservatory &o) const
     61{
     62    return o.RotationAngle(fZd*TMath::DegToRad(), fAz*TMath::DegToRad());
     63}
     64
     65// --------------------------------------------------------------------------
     66//
     67// Get the corresponding rotation angle of the sky coordinate system
     68// seen with an Alt/Az telescope calculated from the stored sky
     69// (Ra/Dec) coordinates.
     70//
     71// For more information see MAstro::RotationAngle
     72//
     73Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t) const
     74{
     75    return MAstroSky2Local(t, o).RotationAngle(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
     76}

  • trunk/MagicSoft/Mars/mpointing/MPointingPos.h

    r3544 r3568  
    55#include "MParContainer.h"
    66#endif
     7
     8class MTime;
     9class MObservatory;
    710
    811class MPointingPos : public MParContainer
    … …  
    3235    Double_t GetDec() const { return fDec; }
    3336
     37    Double_t RotationAngle(const MObservatory &o) const;
     38    Double_t RotationAngle(const MObservatory &o, const MTime &t) const;
     39
    3440    ClassDef(MPointingPos, 1) //Container storing the (corrected) telescope pointing position
    3541};

  • trunk/MagicSoft/Mars/mpointing/Makefile

    r2800 r3568  
    2222#  connect the include files defined in the config.mk file
    2323#
    24 INCLUDES = -I. -I../mbase -I../mraw -I../mreport -I../mmc
     24INCLUDES = -I. -I../mbase -I../mraw -I../mreport -I../mmc -I../mastro
    2525
    2626#------------------------------------------------------------------------------

  • trunk/MagicSoft/Mars/mraw/MRawFileRead.cc

    r3226 r3568  
    3838//                                                                          //
    3939//////////////////////////////////////////////////////////////////////////////
    40 
    4140#include "MRawFileRead.h"
    4241
     42#include <errno.h>
    4343#include <fstream>
    4444
    … …  
    138138    if (!(*fIn))
    139139    {
    140         *fLog << err << "Error: Cannot open file '" << fFileName << "'" << endl;
     140        *fLog << err << "Cannot open file " << fFileName << ": ";
     141        *fLog << strerror(errno) << endl;
    141142        return kFALSE;
    142143    }

  • trunk/MagicSoft/Mars/mreport/MReportFileRead.cc

    r3324 r3568  
    4040#include "MReportFileRead.h"
    4141
     42#include <errno.h>
    4243#include <fstream>
    4344
    … …  
    212213    if (!(*fIn))
    213214    {
    214         *fLog << err << "Error: Cannot open file '" << fFileName << "'" << endl;
     215        *fLog << err << "Cannot open file " << fFileName << ": ";
     216        *fLog << strerror(errno) << endl;
    215217        return kFALSE;
    216218    }
     219
    217220    if (TestBit(kHasNoHeader))
    218221        return kTRUE;
Note: See TracChangeset for help on using the changeset viewer.