Changeset 3666 for trunk/MagicSoft/Mars


Ignore:
Timestamp:
04/06/04 13:41:56 (21 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
2 added
45 edited

Legend:

Unmodified
Added
Removed

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

    r3596 r3666  
    8585Double_t MAstro::Dms2Hor(Int_t deg, UInt_t min, Double_t sec, Char_t sgn)
    8686{
    87     return Hms2Sec(deg, min, sec, sgn)/15.;
     87    return Hms2Sec(deg, min, sec, sgn)/54000.;
    8888}
    8989

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

    r3568 r3666  
    148148
    149149// --------------------------------------------------------------------------
    150 void MAstroCamera::AddPrimitives(Option_t *o)
     150//
     151// Options:
     152//
     153//  '*' Draw the mean of the reflections on all mirrors
     154//  '.' Draw a dot for the reflection on each mirror
     155//
     156void MAstroCamera::AddPrimitives(TString o)
    151157{
    152158    if (!fTime || !fObservatory || !fMirrors)
    … …  
    156162    }
    157163
    158     TString opt(o);
    159     if (opt.IsNull())
    160         opt = "*.";
    161 
    162     const Bool_t hashist = opt.Contains("h", TString::kIgnoreCase);
    163     const Bool_t hasmean = opt.Contains("*", TString::kIgnoreCase);
    164     const Bool_t hasdot  = opt.Contains(".", TString::kIgnoreCase);
    165     const Bool_t usecam  = opt.Contains("c", TString::kIgnoreCase);
     164    if (o.IsNull())
     165        o = "*.";
     166
     167    const Bool_t hashist = o.Contains("h", TString::kIgnoreCase);
     168    const Bool_t hasmean = o.Contains("*", TString::kIgnoreCase);
     169    const Bool_t hasdot  = o.Contains(".", TString::kIgnoreCase);
     170    const Bool_t usecam  = o.Contains("c", TString::kIgnoreCase);
    166171
    167172    // Get camera

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

    r3568 r3666  
    2222
    2323    Int_t  ConvertToPad(const TVector3 &w, TVector2 &v) const;
    24     void   AddPrimitives(Option_t *o);
    25     void   SetRangePad() { }
     24    void   AddPrimitives(TString o);
     25    void   SetRangePad(Option_t *o) { }
    2626    void   ExecuteEvent(Int_t event, Int_t mp1, Int_t mp2);
    2727

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

    r3571 r3666  
    147147MAstroCatalog::~MAstroCatalog()
    148148{
    149     if (fTime)
    150         delete fTime;
    151     if (fObservatory)
    152         delete fObservatory;
    153 
    154     fToolTip->Hide();
    155     delete fToolTip;
    156 
    157     DeleteMap();
    158 
     149    // First disconnect the EventInfo...
    159150    // FIXME: There must be an easier way!
    160151    TIter Next(gROOT->GetListOfCanvases());
    … …  
    164155                      "EventInfo(Int_t,Int_t,Int_t,TObject*)");
    165156
     157    // Now delete the data members
     158    if (fTime)
     159        delete fTime;
     160    if (fObservatory)
     161        delete fObservatory;
     162
     163    fToolTip->Hide();
     164    delete fToolTip;
     165
     166    DeleteMap();
    166167}
    167168
    … …  
    438439void MAstroCatalog::Paint(Option_t *o)
    439440{
    440     SetRangePad();
     441    SetRangePad(o);
    441442
    442443    if (TestBit(kHasChanged))
    443444        DrawPrimitives(o);
     445
     446    PaintMap();
    444447}
    445448
    … …  
    458461
    459462    // draw star on the camera display
    460     TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotLarge);;
     463    TMarker *tip=new TMarker(x, y, transparent ? kDot : kFullDotMedium);;
    461464    tip->SetMarkerColor(kBlack);
    462465    AddMap(tip, new TString(str));
    … …  
    517520        w *= 1./w(2);
    518521
    519     w *= TMath::RadToDeg();
    520     v.Set(w(0), w(1));
     522    w *= TMath::RadToDeg(); // FIXME: *conversion factor?
     523    v.Set(TestBit(kMirrorX) ? -w(0) : w(0),
     524          TestBit(kMirrorY) ? -w(1) : w(1));
    521525
    522526    if (w(2)<0)
    523527        return kERROR;
    524528
    525     return w(0)>gPad->GetX1() && w(1)>gPad->GetY1() &&
    526            w(0)<gPad->GetX2() && w(1)<gPad->GetY2();
     529    return v.X()>gPad->GetX1() && v.Y()>gPad->GetY1() &&
     530           v.X()<gPad->GetX2() && v.Y()<gPad->GetY2();
    527531}
    528532
    … …  
    610614            dx[0] = d[0]+dirs[i][0];
    611615            dy[0] = d[1]+dirs[i][1];
    612 
    613             //cout << dx[0] << " " << dy[0] << endl;
    614616
    615617            // Draw corresponding line and iterate through grid
    … …  
    756758
    757759    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());
     760    txt += Form("  \\theta=%.1f\\circ ", fmod(zd+270,180)-90);
     761    txt += Form("\\phi=%.1f\\circ ",     fmod(az+720, 360));
     762    txt += Form(" / \\rho=%.1f\\circ",   rho*TMath::RadToDeg());
    762763    txt += "}{<";
    763764    txt += fTime->GetSqlDateTime();
    … …  
    767768
    768769// --------------------------------------------------------------------------
    769 void MAstroCatalog::AddPrimitives(Option_t *o)
    770 {
    771     const Bool_t local = TString(o).Contains("local", TString::kIgnoreCase);
    772 
    773     cout << "Opt: " << o << endl;
     770//
     771// To overlay the catalog make sure, that in any case you are using
     772// the 'same' option.
     773//
     774// If you want to overlay this on top of any picture which is created
     775// by derotation of the camera plain you have to use the 'mirror' option
     776// the compensate the mirroring of the image in the camera plain.
     777//
     778// If you have already compensated this by x=-x and y=-y when creating
     779// the histogram you can simply overlay the catalog.
     780//
     781// To overlay the catalog on a 2D histogram the histogram must have
     782// units of degrees (which are plain, like you directly convert the
     783// camera units by multiplication to degrees)
     784//
     785// To be 100% exact you must use the option 'plain' which assumes a plain
     786// screen. This is not necessary for the MAGIC-camera because the
     787// difference between both is less than 1e-3.
     788//
     789// You should always be aware of the fact, that the shown stars and the
     790// displayed grid is the ideal case, like a reflection on a virtual
     791// perfectly aligned central mirror. In reality the star-positions are
     792// smeared to the edge of the camera the more the distance to the center
     793// is, such that the center of gravity of the light distribution might
     794// be more far away from the center than the display shows.
     795//
     796//
     797void MAstroCatalog::AddPrimitives(TString o)
     798{
     799    const Bool_t same   = o.Contains("same",    TString::kIgnoreCase);
     800    const Bool_t local  = o.Contains("local",   TString::kIgnoreCase);
     801    const Bool_t mirx   = o.Contains("mirrorx", TString::kIgnoreCase);
     802    const Bool_t miry   = o.Contains("mirrory", TString::kIgnoreCase);
     803    const Bool_t mirror = o.Contains("mirror",  TString::kIgnoreCase) && !mirx && !miry;
     804
     805    // X is vice versa, because ra is defined anti-clockwise
     806    mirx || mirror ? ResetBit(kMirrorX) : SetBit(kMirrorX);
     807    miry || mirror ? SetBit(kMirrorY) : ResetBit(kMirrorY);
    774808
    775809    const TRotation rot(GetGrid(local));
    … …  
    782816        TVector2 s(v->Phi(), v->Theta());
    783817        if (Convert(rot, s)==kTRUE)
    784             DrawStar(s.X(), -(TMath::Pi()/2-s.Y()), *v, kFALSE);
    785     }
    786 
    787     TPaveText *pv = new TPaveText(0.01, 0.90, 0.63, 0.99, "brNDC");
    788     pv->AddText(GetPadTitle());
    789     AddMap(pv);
     818            DrawStar(s.X(), s.Y(), *v, kFALSE);
     819    }
     820
     821    if (!same)
     822    {
     823        TPaveText *pv = new TPaveText(0.01, 0.90, 0.63, 0.99, "brNDC");
     824        pv->AddText(GetPadTitle());
     825        AddMap(pv);
     826    }
    790827
    791828    TMarker *mk=new TMarker(0, 0, kMultiply);
    … …  
    796833
    797834// --------------------------------------------------------------------------
    798 void MAstroCatalog::SetRangePad()
    799 {
    800      const Double_t edge = fRadiusFOV/TMath::Sqrt(2.);
    801      gPad->Range(-edge, -edge, edge, edge);
    802 
    803      const Float_t w = gPad->GetWw();
    804      const Float_t h = gPad->GetWh();
    805 
    806      if (w<h)
    807          gPad->Range(-edge, -edge*h/w, edge, edge*h/w);
    808      else
    809          gPad->Range(-edge*w/h, -edge, edge*w/h, edge);
     835void MAstroCatalog::SetRangePad(Option_t *o)
     836{
     837    if (TString(o).Contains("same", TString::kIgnoreCase))
     838        return;
     839
     840    const Double_t edge = fRadiusFOV/TMath::Sqrt(2.);
     841    gPad->Range(-edge, -edge, edge, edge);
     842
     843    const Float_t w = gPad->GetWw();
     844    const Float_t h = gPad->GetWh();
     845
     846    if (w<h)
     847        gPad->Range(-edge, -edge*h/w, edge, edge*h/w);
     848    else
     849        gPad->Range(-edge*w/h, -edge, edge*w/h, edge);
    810850}
    811851
    … …  
    815855    DeleteMap();
    816856
    817     SetRangePad();
     857    SetRangePad(o);
    818858
    819859    TStopwatch clk;
    … …  
    827867        AppendPad(o);
    828868
    829     // Append all objects to pad
    830     DrawMap();
    831 
    832869    ResetBit(kHasChanged);
    833870}
    834871
    835872// --------------------------------------------------------------------------
    836 void MAstroCatalog::DrawMap()
     873void MAstroCatalog::PaintMap()
    837874{
    838875    Long_t key, val;
    839876    TExMapIter map(&fMapG);
    840877    while (map.Next(key, val))
    841         ((TObject*)key)->Draw();
     878        ((TObject*)key)->Paint();
    842879}
    843880
    … …  
    10161053    fToolTip->Show(x+4, y+4);
    10171054}
     1055/*
     1056void MAstroCatalog::RecursiveRemove(TObject *obj)
     1057{
     1058    ULong_t hash;
     1059    Long_t key, val;
     1060
     1061    TExMapIter map(&fMapG);
     1062    while (map.Next(hash, key, val))
     1063    {
     1064        if (key != (Long_t)obj)
     1065            continue;
     1066
     1067        fMapG.Remove(hash, key);
     1068        delete (TObject*)(key);
     1069        if (val)
     1070            delete (TString*)(val);
     1071        break;
     1072    }
     1073}
     1074*/

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

    r3568 r3666  
    9494
    9595    virtual Int_t ConvertToPad(const TVector3 &w, TVector2 &v) const;
    96     virtual void  AddPrimitives(Option_t *o);
    97     virtual void  SetRangePad();
     96    virtual void  AddPrimitives(TString o);
     97    virtual void  SetRangePad(Option_t *o);
    9898
    9999    Int_t     Convert(const TRotation &rot, TVector2 &v) const;
    … …  
    105105    void      Paint(Option_t *o="");
    106106    Int_t     DistancetoPrimitive(Int_t px, Int_t py);
    107     void      DrawMap();
     107    //void      RecursiveRemove(TObject *obj);
     108    void      PaintMap();
    108109    void      DeleteMap()
    109110    {
    … …  
    117118
    118119            delete (TString*)(val);
    119          /*
    120             Long_t key2, val2;
    121             TExMapIter map2(&fMapG);
    122             while (map2.Next(key2, val2))
    123                 if (val==val2)
    124                 {
    125                     delete (TObject*)key;
    126                     fMapG.Remove(key);
    127                 }*/
     120            /*
     121              Long_t key2, val2;
     122              TExMapIter map2(&fMapG);
     123              while (map2.Next(key2, val2))
     124                  if (val==val2)
     125                  {
     126                      delete (TObject*)key;
     127                      fMapG.Remove(key);
     128                  }
     129             */
    128130        }
    129131        fMapG.Delete();
    … …  
    134136        kHasChanged  = BIT(15),
    135137        kGuiActive   = BIT(16),
    136         kPlainScreen = BIT(17)
     138        kPlainScreen = BIT(17),
     139        kMirrorX     = BIT(18),
     140        kMirrorY     = BIT(19)
    137141    };
    138142
    … …  
    140144    MVector3 fRaDec;     // pointing position
    141145
    142     MObservatory *fObservatory; // Possible obervatora location
     146    MObservatory *fObservatory; // Possible obervatory location
    143147    MTime        *fTime;        // Possible observation time
    144148
    … …  
    153157    void AddMap(TObject *k, void *v=0)
    154158    {
    155         k->SetBit(kCanDelete);
    156         k->SetBit(kCannotPick);
     159        //k->SetBit(kCanDelete);
     160        //k->SetBit(kCannotPick);
    157161        fMapG.Add(fMapG.GetSize(), (Long_t)k, (Long_t)v);
    158162    }
    … …  
    164168    void SetTime(const MTime &time);
    165169    void SetObservatory(const MObservatory &obs);
    166     void SetLimMag(Double_t mag) { fLimMag=mag;    Update(); } // *MENU* *ARGS={mag=>fLimMag}
     170    void SetLimMag(Double_t mag) { fLimMag=mag; Update(); } // *MENU* *ARGS={mag=>fLimMag}
    167171    void SetRadiusFOV(Double_t deg)
    168172    {

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

    r3568 r3666  
    333333
    334334    //
    335     // Update current speed each second
    336     //
    337     if (fDisplay && t0-t2>(TTime)1000)
    338     {
    339         const Int_t speed = 1000*(num-start)/(long int)(t0-t2);
     335    // Update current speed each 1.5 second
     336    //
     337    if (fDisplay && t0-t2>(TTime)1500)
     338    {
     339        const Float_t speed = 1000.*(num-start)/(long int)(t0-t2);
    340340        TString txt = "Processing...";
    341341        if (speed>0)
    342342        {
    343343            txt += " (";
    344             txt += speed;
     344            txt += (Int_t)speed;
    345345            txt += "Evts/s";
    346346            if (fNumEvents>0)
    347347            {
    348348                txt += ", est: ";
    349                 txt += (int)((long int)(t0-t2)*fNumEvents/(60000.*(num-start)));
    350                 txt += "m";
     349                txt += (int)((fNumEvents-num)/speed/60)+1;
     350                txt += "min";
    351351            }
    352352            //txt += (int)fmod(entries/(1000.*(num-start)/(long int)(t0-t2)), 60);
    … …  
    356356        fDisplay->SetStatusLine1(txt);
    357357        start = num;
    358         t2 = t1;
     358        t2 = t0;
    359359    }
    360360

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

    r3597 r3666  
    5252// --------------------------------------------------------------------------
    5353//
     54// Symmetrized significance - this is somehow analog to
     55// SignificanceLiMaSigned
     56//
     57// Returns Significance(s,b) if s>b otherwise -Significance(b, s);
     58//
     59Double_t MMath::SignificanceSym(Double_t s, Double_t b)
     60{
     61    return s>b ? Significance(s, b) : -Significance(b, s);
     62}
     63
     64// --------------------------------------------------------------------------
     65//
    5466//  calculates the significance according to Li & Ma
    5567//  ApJ 272 (1983) 317, Formula 17
    … …  
    5870//  b                    // b: number of off events
    5971//  alpha = t_on/t_off;  // t: observation time
     72//
     73//  The significance has the same (positive!) value for s>b and b>s.
    6074//
    6175//  Returns -1 if sum<0 or alpha<0 or the argument of sqrt<0
    … …  
    7791    return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
    7892}
     93
     94// --------------------------------------------------------------------------
     95//
     96// Calculates MMath::SignificanceLiMa(s, b, alpha). Returns 0 if the
     97// calculation has failed. Otherwise the Li/Ma significance which was
     98// calculated. If s<b a negative value is returned.
     99//
     100Double_t MMath::SignificanceLiMaSigned(Double_t s, Double_t b, Double_t alpha)
     101{
     102    const Double_t sig = SignificanceLiMa(s, b, alpha);
     103    if (sig<=0)
     104        return 0;
     105
     106    return TMath::Sign(sig, s-b);
     107}

  • trunk/MagicSoft/Mars/mbase/MMath.h

    r3582 r3666  
    1010public:
    1111    static Double_t Significance(Double_t s, Double_t b);
     12    static Double_t SignificanceSym(Double_t s, Double_t b);
    1213    static Double_t SignificanceLiMa(Double_t s, Double_t b, Double_t alpha=1);
     14    static Double_t SignificanceLiMaSigned(Double_t s, Double_t b, Double_t alpha=1);
    1315
    1416    ClassDef(MMath, 0)

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

    r3574 r3666  
    407407TMethodCall *MParContainer::GetterMethod(const char *name) const
    408408{
    409     TClass *cls = IsA()->GetBaseDataMember(name);
     409    const TString n(name);
     410    const Int_t pos1 = n.First('.');
     411
     412    const TString part1 = pos1<0 ? n : n(0, pos1);
     413    const TString part2 = pos1<0 ? TString("") : n(pos1+1, n.Length());
     414
     415    TClass *cls = IsA()->GetBaseDataMember(part1);
    410416    if (cls)
    411417    {
    412         TDataMember *member = cls->GetDataMember(name);
     418        TDataMember *member = cls->GetDataMember(part1);
    413419        if (!member)
    414420        {
    415             *fLog << err << "Datamember '" << name << "' not in " << GetDescriptor() << endl;
     421            *fLog << err << "Datamember '" << part1 << "' not in " << GetDescriptor() << endl;
    416422            return NULL;
     423        }
     424
     425        // This handles returning references of contained objects, eg
     426        // class X { TObject fO; TObject &GetO() { return fO; } };
     427        if (!member->IsBasic() && !part2.IsNull())
     428        {
     429            cls = gROOT->GetClass(member->GetTypeName());
     430            if (!cls)
     431            {
     432                *fLog << err << "Datamember " << part1 << " [";
     433                *fLog << member->GetTypeName() << "] not in dictionary." << endl;
     434                return NULL;
     435            }
     436            if (!cls->InheritsFrom(MParContainer::Class()))
     437            {
     438                *fLog << err << "Datamember " << part1 << " [";
     439                *fLog << member->GetTypeName() << "] does not inherit from ";
     440                *fLog << "MParContainer." << endl;
     441                return NULL;
     442            }
     443
     444            const MParContainer *sub = (MParContainer*)((ULong_t)this+member->GetOffset());
     445            return sub->GetterMethod(part2);
     446        }
     447
     448        if (member->IsaPointer())
     449        {
     450            *fLog << warn << "Data-member " << part1 << " is a pointer..." << endl;
     451            *fLog << dbginf << "Not yet implemented!" << endl;
     452            //TClass *test = gROOT->GetClass(member->GetTypeName());
     453            return 0;
    417454        }
    418455
    … …  
    422459    }
    423460
    424     *fLog << warn << "No standard access for '" << name << "' in ";
     461    *fLog << warn << "No standard access for '" << part1 << "' in ";
    425462    *fLog << GetDescriptor() << " or one of its base classes." << endl;
    426463
    … …  
    428465
    429466    *fLog << warn << "Trying to find MethodCall '" << ClassName();
    430     *fLog << "::Get" << name << "' instead <LEAKS MEMORY>" << endl;
    431     call = new TMethodCall(IsA(), (TString)"Get"+name, "");
     467    *fLog << "::Get" << part1 << "' instead <LEAKS MEMORY>" << endl;
     468    call = new TMethodCall(IsA(), (TString)"Get"+part1, "");
    432469    if (call->GetMethod())
    433470        return call;
    … …  
    436473
    437474    *fLog << warn << "Trying to find MethodCall '" << ClassName();
    438     *fLog << "::" << name << "' instead <LEAKS MEMORY>" << endl;
    439     call = new TMethodCall(IsA(), name, "");
     475    *fLog << "::" << part1 << "' instead <LEAKS MEMORY>" << endl;
     476    call = new TMethodCall(IsA(), part1, "");
    440477    if (call->GetMethod())
    441478        return call;
    … …  
    443480    delete call;
    444481
    445     *fLog << err << "Sorry, no getter method found for " << name << endl;
     482    *fLog << err << "Sorry, no getter method found for " << part1 << endl;
    446483    return NULL;
    447484}

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

    r3512 r3666  
    14331433    //        where (eg Eventloop) first!
    14341434
    1435     gLog << dbg << fName << " is on heap: " << (int)IsOnHeap() << endl;
     1435    //gLog << dbg << fName << " is on heap: " << (int)IsOnHeap() << endl;
    14361436
    14371437    if (TestBit(kExitLoopOnExit) || TestBit(kExitLoopOnClose))
    … …  
    14431443    if (fIsLocked<=0 && IsOnHeap())
    14441444    {
    1445         gLog << dbg << "delete " << fName << ";" << endl;
     1445        //gLog << dbg << "delete " << fName << ";" << endl;
    14461446        delete this;
    14471447    }
    14481448    fStatus = kFileExit;
    1449     gLog << dbg << fName << ".fStatus=kFileExit;" << endl;
     1449    //gLog << dbg << fName << ".fStatus=kFileExit;" << endl;
    14501450}
    14511451

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

    r3497 r3666  
    106106MTask::MTask(const char *name, const char *title)
    107107    : fFilter(NULL), fSerialNumber(0), fIsPreprocessed(kFALSE),
    108     fNumExecutions(0), fStopwatch(0)
     108    fStopwatch(0)
    109109{
    110110    fName  = name  ? name  : "MTask";
    … …  
    194194Int_t MTask::CallPreProcess(MParList *plist)
    195195{
    196     fNumExecutions = 0;
    197196    fStopwatch->Reset();
    198197
    … …  
    240239        return kTRUE;
    241240
    242     fNumExecutions++;
    243 
    244241    fStopwatch->Start(kFALSE);
    245242    const Int_t rc = Process();
    … …  
    346343// --------------------------------------------------------------------------
    347344//
    348 //  Return total CPU execution time of task
     345//  Return the total number of calls to Process(). If Process() was not
     346//  called due to a set filter this is not counted.
     347//
     348UInt_t MTask::GetNumExecutions() const
     349{
     350    return (UInt_t)fStopwatch->Counter();
     351}
     352
     353// --------------------------------------------------------------------------
     354//
     355//  Return total CPU execution time in seconds of calls to Process().
     356//  If Process() was not called due to a set filter this is not counted.
    349357//
    350358Double_t MTask::GetCpuTime() const
    … …  
    355363// --------------------------------------------------------------------------
    356364//
    357 //  Return total real execution time of task
     365//  Return total real execution time in seconds of calls to Process().
     366//  If Process() was not called due to a set filter this is not counted.
    358367//
    359368Double_t MTask::GetRealTime() const
    … …  
    366375//  Prints the relative time spent in Process() (relative means relative to
    367376//  its parent Tasklist) and the number of times Process() was executed.
     377//  Don't wonder if the sum of the tasks in a tasklist is not 100%,
     378//  because only the call to Process() of the task is measured. The
     379//  time of the support structure is ignored. The faster your analysis is
     380//  the more time is 'wasted' in the support structure.
     381//  Only the CPU time is displayed. This means that exspecially task
     382//  which have a huge part of file i/o will be underestimated in their
     383//  relative wasted time.
    368384//  For convinience the lvl argument results in a number of spaces at the
    369385//  beginning of the line. So that the structur of a tasklist can be
    … …  
    379395
    380396    *fLog << all << setfill(' ') << setw(lvl) << " ";
    381     if (fStopwatch->CpuTime()>0 && time>0)
    382         *fLog << setw(3) << (Int_t)(fStopwatch->CpuTime()/time*100+.5) << "% ";
     397
     398    if (GetCpuTime()>0 && time>0 && GetCpuTime()>=0.001*time)
     399        *fLog << Form("%5.1f", GetCpuTime()/time*100) << "% ";
    383400    else
    384         *fLog << "     ";
     401        *fLog << "       ";
    385402    *fLog << GetDescriptor() << "\t";
    386     *fLog << dec << fNumExecutions;
     403    *fLog << dec << GetNumExecutions();
    387404    if (fFilter)
    388405        *fLog << " <" << fFilter->GetName() << ">";

  • trunk/MagicSoft/Mars/mbase/MTask.h

    r3497 r3666  
    2929
    3030    Bool_t fIsPreprocessed; //! Indicates the success of the PreProcessing (set by MTaskList)
    31     UInt_t fNumExecutions;  //! Number of Excutions
    3231
    33     TStopwatch *fStopwatch; //!
     32    TStopwatch *fStopwatch; //! Count the execution time and number of executions
    3433
    3534    virtual Int_t PreProcess(MParList *pList);
    … …  
    7473    // Filter functions
    7574    virtual void SetFilter(MFilter *filter) { fFilter=filter; }
    76     const MFilter *GetFilter() const      { return fFilter; }
    77     MFilter *GetFilter()  { return fFilter; } // for MContinue only
     75    const MFilter *GetFilter() const        { return fFilter; }
     76    MFilter *GetFilter()                    { return fFilter; } // for MContinue only
    7877
    7978    // Display functions
    … …  
    8685    TString AddSerialNumber(const TString &str) const { return AddSerialNumber(str, fSerialNumber); }
    8786
    88     virtual void SetSerialNumber(Byte_t num) { fSerialNumber = num; }
    89     Byte_t  GetSerialNumber() const { return fSerialNumber; }
     87    virtual void SetSerialNumber(Byte_t num) { fSerialNumber = num;  }
     88    Byte_t GetSerialNumber() const          { return fSerialNumber; }
    9089
    9190    const char *GetDescriptor() const;
    9291
    9392    // Task execution statistics
    94     UInt_t GetNumExecutions() const { return fNumExecutions; }
     93    UInt_t   GetNumExecutions() const;
    9594    Double_t GetCpuTime() const;
    9695    Double_t GetRealTime() const;

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

    r3497 r3666  
    635635        if (title)
    636636            *fLog << "\t" << fTitle;
     637        if (time>=0)
     638            *fLog << Form(" %5.1f", GetCpuTime()/time*100) << "%";
     639        else
     640            *fLog << " 100.0%";
    637641        *fLog << endl;
    638642    }

  • trunk/MagicSoft/Mars/mbase/MTaskList.h

    r3497 r3666  
    7070
    7171    void Print(Option_t *opt = "") const;
    72     void PrintStatistics(const Int_t lvl=0, Bool_t title=kFALSE, Double_t time=0) const;
     72    void PrintStatistics(const Int_t lvl=0, Bool_t title=kFALSE, Double_t time=-1) const;
    7373    void SetOwner(Bool_t enable=kTRUE);
    7474

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

    r3371 r3666  
    362362    const Double_t sum = (r1+r2)/(24*3600);
    363363
    364     return fmod(sum, 1)*TMath::TwoPi()+TMath::TwoPi();
     364    return fmod(sum, 1)*TMath::TwoPi();//+TMath::TwoPi();
    365365}
    366366

  • trunk/MagicSoft/Mars/mcamera/MCameraAUX.h

    r2593 r3666  
    1414    Bool_t fStatusCaosLEDs;  // Monitored status: o=off, 1=on, Cam.CaOs.LED_state
    1515    Bool_t fStatusFansFADC;  // Monitored status: o=off, 1=on, Cam.FADC.Fans_state
     16
    1617public:
    1718    MCameraAUX()
    … …  
    2021        fTitle = "Container storing information about the Camera auxiliary system";
    2122    }
     23
     24    Bool_t GetRequestCaosLEDs() const { return fRequestCaosLEDs; }
     25    Bool_t GetRequestFansFADC() const { return fRequestFansFADC; }
     26    Bool_t GetStatusCaosLEDs() const  { return fStatusCaosLEDs;  }
     27    Bool_t GetStatusFansFADC() const  { return fStatusFansFADC;  }
     28
    2229    ClassDef(MCameraAUX, 1) // Container storing information about the Camera auxiliary system
    2330};

  • trunk/MagicSoft/Mars/mcamera/MCameraCalibration.h

    r2864 r3666  
    3131    Byte_t GetStatusIO() const         { return fStatusIO; }
    3232    Byte_t GetStatusLoVoltage() const  { return fStatusLoVoltage; }
     33
    3334    Bool_t GetRequestHiVoltage() const { return fRequestHiVoltage; }
    3435    Bool_t GetRequestLoVoltage() const { return fRequestLoVoltage; }

  • trunk/MagicSoft/Mars/mcamera/MCameraHV.h

    r3216 r3666  
    2828    TArrayS fHV;            // [V] Measured high Voltages for all PMTs
    2929
    30     Float_t fMean;          // [V] Mean high voltage of the camera
    31 
    3230public:
    3331    MCameraHV() : fHV(577)
    … …  
    3937    Byte_t  GetStatus() const { return fStatus; }
    4038    Bool_t  GetStatusRamping() const { return fStatusRamping; }
     39
    4140    Short_t GetVoltageA() const { return fVoltageA; }
    4241    Short_t GetVoltageB() const { return fVoltageB; }
    43     Byte_t GetCurrentA() const { return fCurrentA; }
    44     Byte_t GetCurrentB() const { return fCurrentB; }
    4542
    46     Float_t GetMean() { fMean = fHV.GetSum()/fHV.GetSize(); return fMean; }
     43    Byte_t  GetCurrentA() const { return fCurrentA; }
     44    Byte_t  GetCurrentB() const { return fCurrentB; }
     45
     46    Float_t GetMean() const { return fHV.GetSum()/fHV.GetSize(); }
    4747
    4848    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const
    49       {
    50         val = fHV[idx];
    51         return val>0;
    52       }
     49    {
     50        val = fHV[idx];
     51        return val>0;
     52    }
    5353    void DrawPixelContent(Int_t num) const
    54       {
    55       }
     54    {
     55    }
    5656
    5757    ClassDef(MCameraHV, 1) // Container storing information about the Camera HV

  • trunk/MagicSoft/Mars/mcamera/MCameraLV.h

    r2895 r3666  
    1414    friend class MReportCamera;
    1515private:
    16     Byte_t  fStatus;                    // CaCo monitored LV PS status:  , Cam.LV_state
    17     Bool_t  fRequestPowerSupply;        // Requested status: o=off, 1=on, blv_ps_status
     16    Byte_t  fStatus;                  // CaCo monitored LV PS status:  , Cam.LV_state
     17    Bool_t  fRequestPowerSupply;      // Requested status: o=off, 1=on, blv_ps_status
    1818
    19     Float_t fTemp;                     // Measured status: o=off, 1=on, blv_temp
    20     Byte_t  fHumidity;                 // Measured status: o=off, 1=on, blv_RelativeHumidity
     19    Float_t fTemp;                    // Measured status: o=off, 1=on, blv_temp
     20    Byte_t  fHumidity;                // Measured status: o=off, 1=on, blv_RelativeHumidity
    2121
    22     MCameraPowerSupply fPowerSupplyA;   // power supply camera part A
    23     MCameraPowerSupply fPowerSupplyB;   // power supply camera part B
     22    MCameraPowerSupply fPowerSupplyA; // power supply camera part A
     23    MCameraPowerSupply fPowerSupplyB; // power supply camera part B
    2424
    2525public:
    … …  
    3232    Byte_t  GetStatus() const { return fStatus; }
    3333    Bool_t  GetRequestPowerSupply() const { return fRequestPowerSupply; }
     34
    3435    Float_t GetTemp() const { return fTemp; }
    3536    Byte_t  GetHumidity() const { return fHumidity; }
     37
     38    const MCameraPowerSupply &GetPowerSupplyA() const { return fPowerSupplyA; }
     39    const MCameraPowerSupply &GetPowerSupplyB() const { return fPowerSupplyB; }
    3640
    3741    ClassDef(MCameraLV, 1) // Container storing information about the Camera LV

  • trunk/MagicSoft/Mars/mcamera/MCameraLid.h

    r2593 r3666  
    1414     Bool_t fLimitOpen;        // 0=not active, 1= active
    1515     Bool_t fLimitClose;       // 0=not active, 1= active
     16
    1617     Bool_t fSafetyLimitOpen;  // 0=not active, 1= active
    1718     Bool_t fSafetyLimitClose; // 0=not active, 1= active
     19
    1820     Byte_t fStatusLid;        // 0=positioning, 1=open, 2=closed
    1921     Byte_t fStatusMotor;      // 0=stopped, 1=opening, 2=closing
     22
    2023public:
    2124    MCameraLid()
    … …  
    2427        fTitle = "Container storing information about a Camera lid";
    2528    }
     29
     30    Bool_t GetLimitOpen() const        { return fLimitOpen;        }
     31    Bool_t GetLimitClose() const       { return fLimitClose;       }
     32
     33    Bool_t GetSafetyLimitOpen() const  { return fSafetyLimitOpen;  }
     34    Bool_t GetSafetyLimitClose() const { return fSafetyLimitClose; }
     35
     36    Byte_t GetStatusLid() const        { return fStatusLid;        }
     37    Byte_t GetStatusMotor() const      { return fStatusMotor;      }
     38
    2639    ClassDef(MCameraLid, 1) // Container storing information about a Camera lid
    2740};

  • trunk/MagicSoft/Mars/mcamera/MCameraLids.h

    r3066 r3666  
    2727    Byte_t GetStatus() const { return fStatus; }
    2828
     29    const MCameraLid &GetLidA() const { return fLidA; }
     30    const MCameraLid &GetLidB() const { return fLidB; }
     31
    2932    ClassDef(MCameraLids, 1) // Container storing information about the Camera lids
    3033};

  • trunk/MagicSoft/Mars/mcamera/MCameraPowerSupply.h

    r2593 r3666  
    1010    friend class MReportCamera;
    1111private:
    12      Float_t fVoltagePos5V;         // [V] voltage_pos5  (+5V)
    13      Float_t fVoltagePos12V;        // [V] voltage_pos12 (+12V)
    14      Float_t fVoltageNeg12V;        // [V] voltage_neg12 (-12V)
    15      Float_t fVoltageOptLinkPos12V; // [V] volatge_opt_link_pos12 (+12V)
     12    Float_t fVoltagePos5V;         // [V] voltage_pos5  (+5V)
     13    Float_t fVoltagePos12V;        // [V] voltage_pos12 (+12V)
     14    Float_t fVoltageNeg12V;        // [V] voltage_neg12 (-12V)
     15    Float_t fVoltageOptLinkPos12V; // [V] volatge_opt_link_pos12 (+12V)
    1616
    17      Float_t fCurrentPos5V;         // [A] current_pos5  (+5V)
    18      Float_t fCurrentPos12V;        // [A] current_pos12 (+12V)
    19      Float_t fCurrentNeg12V;        // [A] current_neg12 (-12V)
    20      Float_t fCurrentOptLinkPos12V; // [A] current_opt_link_pos12 (+12V)
     17    Float_t fCurrentPos5V;         // [A] current_pos5  (+5V)
     18    Float_t fCurrentPos12V;        // [A] current_pos12 (+12V)
     19    Float_t fCurrentNeg12V;        // [A] current_neg12 (-12V)
     20    Float_t fCurrentOptLinkPos12V; // [A] current_opt_link_pos12 (+12V)
    2121
    2222public:
    … …  
    2626        fTitle = "Container storing information about the Camera power supply";
    2727    }
     28
     29    Float_t GetVoltagePos5V() const         { return fVoltagePos5V; }
     30    Float_t GetVoltagePos12V() const        { return fVoltagePos12V; }
     31    Float_t GetVoltageNeg12V() const        { return fVoltageNeg12V; }
     32    Float_t GetVoltageOptLinkPos12V() const { return fVoltageOptLinkPos12V; }
     33
     34    Float_t GetCurrentPos5V() const         { return fCurrentPos5V; }
     35    Float_t GetCurrentPos12V() const        { return fCurrentPos12V; }
     36    Float_t GetCurrentNeg12V() const        { return fCurrentNeg12V; }
     37    Float_t GetCurrentOptLinkPos12V() const { return fCurrentOptLinkPos12V; }
     38
    2839    ClassDef(MCameraPowerSupply, 1) // Container storing information about the Camera power supply
    2940};

  • trunk/MagicSoft/Mars/mdata/DataLinkDef.h

    r1524 r3666  
    1212#pragma link C++ class MDataMember+;
    1313#pragma link C++ class MDataChain+;
     14#pragma link C++ class MDataFormula+;
    1415
    1516#endif

  • trunk/MagicSoft/Mars/mdata/MDataChain.cc

    r3572 r3666  
    2525/////////////////////////////////////////////////////////////////////////////
    2626//
    27 //   MDataChain
     27// MDataChain
     28// ==========
    2829//
    2930// With this chain you can concatenate simple mathematical operations on
    3031// members of mars containers.
     32//
     33//
     34// Rules
     35// -----
    3136//
    3237// In the constructor you can give rule, like
    … …  
    3641// in the containers. The result will be fDist divided by fLength.
    3742//
     43// In case you want to access a data-member which is a data member object
     44// you can acces it with (Remark: it must derive from MParContainer):
     45//   "MCameraLV.fPowerSupplyA.fVoltagePos5V"
     46//
    3847// You can also use brackets:
    3948//   "HillasDource.fDist / (MHillas.fLength + MHillas.fWidth)"
    4049//
    41 // The allowed operations are: +, -, *, /, %
     50//
     51// Operators
     52// ---------
     53//
     54// The allowed operations are: +, -, *, /, %, ^
    4255//
    4356// While a%b returns the floating point reminder of a/b.
     57// While a^b returns a to the power of b
    4458//
    4559// Warning: There is no priority rule build in. So better use brackets
    … …  
    98112// inf is the symbol for an infinite number.
    99113//
     114//
     115// Constants
     116// ---------
     117//
    100118// Constants are implemented in ParseDataMember, namely:
    101 //   kPi:      TMath::Pi()
    102 //   kRad2Deg: 180/kPi
    103 //   kDeg2Rad: kPi/180
     119//   kPi:       TMath::Pi()
     120//   kRad2Deg:  180/kPi
     121//   kDeg2Rad:  kPi/180
    104122//
    105123// You can also defined constants which are defined in TMath by:
    106 //   kLn10      for   static Double_t TMath::Ln10();
    107 //   kLogE      for   static Double_t TMath::LogE();
    108 //   kRadToDeg  for   static Double_t TMath::RadToDeg();
    109 //   kDegToRad  for   static Double_t TMath::DegToRad();
    110 //
    111 // Remark: In older root versions only Pi() and E() are implemented in
    112 //         TMath.
    113 //
     124//   kLn10       for   static Double_t TMath::Ln10();
     125//   kLogE       for   static Double_t TMath::LogE();
     126//   kRadToDeg   for   static Double_t TMath::RadToDeg();
     127//   kDegToRad   for   static Double_t TMath::DegToRad();
     128//   ...
     129//
     130// Remark:
     131//  In older root versions only Pi() and E() are implemented
     132//  in TMath.
     133//
     134//
     135// Variable Parameters
     136// ------------------------
    114137// If you want to use variables, eg for fits you can use [0], [1], ...
    115138// These values are initialized with 0 and set by calling
    116139// SetVariables().
    117140//
     141//
     142// Multi-argument functions
     143// ------------------------
     144// You can use multi-argument functions, too. The access is implemented
     145// via TFormula, which slows down thing a lot. If you can avoid usage
     146// of such expression you accelerate the access a lot. Example:
     147//   "TMath::Hypot(MHillas.fMeanX, MHillas.MeanY)"
     148//   "pow(MHillas.fMeanX*MHillas.MeanY, -1.2)"
     149// It should be possible to use all functions which are accessible
     150// via the root-dictionary.
     151//
     152//
    118153// REMARK:
    119 //         - All the random functions are returning 0 if gRandom==0
    120 //         - You may get better results if you are using a TRandom3
    121 //
    122 // FIXME: The possibility to use other objects inheriting from MData
    123 //        is missing.
    124 //        Maybe we can use gInterpreter->Calc("") for this.
    125 //        gROOT->ProcessLineFast("line");
     154//  - All the random functions are returning 0 if gRandom==0
     155//  - You may get better results if you are using a TRandom3
     156//
     157// To Do:
     158//  - The possibility to use other objects inheriting from MData
     159//    is missing.
     160//  - By automatic pre-adding brackets to the rule it would be possible
     161//    to implement priorities for operators.
    126162//
    127163/////////////////////////////////////////////////////////////////////////////
    … …  
    141177#include "MDataValue.h"
    142178#include "MDataMember.h"
     179#include "MDataFormula.h"
    143180#include "MDataElement.h"
    144181
    … …  
    222259    for (int i=0; i<l; i++)
    223260    {
    224         if (!isalnum(txt[i]) && txt[i]!='.' && txt[i]!=';' &&
     261        if (!isalnum(txt[i]) && txt[i]!='.' && txt[i]!=':' && txt[i]!=';' &&
    225262            /*txt[i]!='['&&txt[i]!=']'&&*/
    226263            ((txt[i]!='-' && txt[i]!='+') || i!=0))
    … …  
    394431        case '/':
    395432        case '%':
     433        case '^':
    396434            if (member0)
    397435            {
    … …  
    507545
    508546            OperatorType_t op = ParseOperator(text);
     547
     548            Int_t first = GetBracket(txt, '(', ')');
     549            TString sub = op==kENegative || op==kEPositive ? text.Remove(0,1) + txt : txt(1, first-1);
     550            txt.Remove(0, first+1);
     551
    509552            if (op==kENoop)
    510553            {
     554                newmember = new MDataFormula(Form("%s(%s)", (const char*)text, (const char*)sub));
     555                if (newmember->IsValid())
     556                    break;
     557
    511558                *fLog << err << dbginf << "Syntax Error: Operator '" << text << "' unknown." << endl;
    512559                if (member0)
    … …  
    514561                return NULL;
    515562            }
    516 
    517             Int_t first = GetBracket(txt, '(', ')');
    518             TString sub = op==kENegative || op==kEPositive ? text.Remove(0,1) + txt : txt(1, first-1);
    519             txt.Remove(0, first+1);
    520563
    521564            newmember = new MDataChain(sub, op);
    … …  
    604647    return 0;
    605648}
    606 
    607     /*
    608 void MDataChain::Print(Option_t *opt) const
    609 {
    610     *fLog << GetRule() << flush;
    611     Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
    612 
    613     switch (fOperatorType)
    614     {
    615     case kEAbs:      *fLog << "abs"   << flush; break;
    616     case kELog:      *fLog << "log"   << flush; break;
    617     case kELog10:    *fLog << "log10" << flush; break;
    618     case kESin:      *fLog << "sin"   << flush; break;
    619     case kECos:      *fLog << "cos"   << flush; break;
    620     case kETan:      *fLog << "tan"   << flush; break;
    621     case kESinH:     *fLog << "sinh"  << flush; break;
    622     case kECosH:     *fLog << "cosh"  << flush; break;
    623     case kETanH:     *fLog << "tanh"  << flush; break;
    624     case kEASin:     *fLog << "asin"  << flush; break;
    625     case kEACos:     *fLog << "acos"  << flush; break;
    626     case kEATan:     *fLog << "atan"  << flush; break;
    627     case kESqrt:     *fLog << "sqrt"  << flush; break;
    628     case kEExp:      *fLog << "exp"   << flush; break;
    629     case kEPow10:    *fLog << "pow10" << flush; break;
    630     case kESgn:      *fLog << "sgn"   << flush; break;
    631     case kENegative: *fLog << "-" << flush; break;
    632     case kEPositive: *fLog << "+" << flush; break;
    633     case kENoop:
    634         break;
    635     }
    636 
    637     if (bracket)
    638         *fLog << "(" << flush;
    639 
    640     fMember->Print();
    641 
    642     if (bracket)
    643         *fLog << ")" << flush;
    644         }
    645         */
    646649
    647650// --------------------------------------------------------------------------

  • trunk/MagicSoft/Mars/mdata/MDataChain.h

    r3572 r3666  
    1717{
    1818private:
    19     MData *fMember; // Filter
     19    MData    *fMember; // Filter
    2020
    2121    // PLEASE, always add new enums to the end of the enumeration,

  • trunk/MagicSoft/Mars/mdata/MDataList.cc

    r3572 r3666  
    7979        fSign = kEModul;
    8080        return;
     81    case '^':
     82        fSign = kEPow;
     83        return;
    8184    default:
    8285        fSign = kENone;
    … …  
    142145        while ((member=(MData*)Next()))
    143146        {
    144             Double_t d = member->GetValue();
     147            const Double_t d = member->GetValue();
    145148            if (d==0)
    146149            {
    147                 *fLog << warn << "Warning: Division by zero (" << member->GetName() << ")" << endl;
     150                *fLog << warn << "Warning: Division by zero: " << member->GetRule() << endl;
    148151                return 0;
    149152            }
    … …  
    155158        while ((member=(MData*)Next()))
    156159        {
    157             Double_t d = member->GetValue();
     160            const Double_t d = member->GetValue();
    158161            if (d==0)
    159162            {
    160                 *fLog << warn << "Warning: Modulo division by zero (" << member->GetName() << ")" << endl;
     163                *fLog << warn << "Warning: Modulo division by zero: " << member->GetRule() << endl;
    161164                return 0;
    162165            }
    … …  
    164167        }
    165168        break;
     169    case kEPow:
     170        while ((member=(MData*)Next()))
     171            val = pow(val, member->GetValue());
     172        break;
    166173    }
    167174    return val;
    … …  
    185192}
    186193
     194Bool_t MDataList::IsValid() const
     195{
     196    TIter Next(&fMembers);
     197
     198    MData *data = NULL;
     199    while ((data=(MData*)Next()))
     200        if (!data->IsValid())
     201            return kFALSE;
     202
     203    return kTRUE;
     204}
     205
    187206// --------------------------------------------------------------------------
    188207//
    189208// If you want to add a new member to the list call this function with the
    190 // pointer to the member to be added. 
     209// pointer to the member to be added.
    191210//
    192211Bool_t MDataList::AddToList(MData *member)
    … …  
    335354            str += "%";
    336355            break;
     356
     357        case kEPow:
     358            str += "^";
     359            break;
    337360        }
    338361

  • trunk/MagicSoft/Mars/mdata/MDataList.h

    r3572 r3666  
    2323    TOrdCollection fMembers;    // Container for the filters
    2424
    25     typedef enum { kENone, kEPlus, kEMinus, kEMult, kEDiv, kEModul } SignType_t;
     25    typedef enum { kENone, kEPlus, kEMinus, kEMult, kEDiv, kEModul, kEPow } SignType_t;
    2626    SignType_t fSign;
    2727
    … …  
    4141    void SetOwner(Bool_t enable=kTRUE) { enable ? SetBit(kIsOwner) : ResetBit(kIsOwner); }
    4242
    43     Bool_t IsValid() const { return fMembers.GetSize() ? kTRUE : kFALSE; }
     43    Bool_t IsValid() const;// { return fMembers.GetSize() ? kTRUE : kFALSE; }
    4444    Bool_t IsReadyToSave() const;
    4545

  • trunk/MagicSoft/Mars/mdata/Makefile

    r2800 r3666  
    3737           MDataValue.cc \
    3838           MDataList.cc \
    39            MDataChain.cc
     39           MDataChain.cc \
     40           MDataFormula.cc
    4041
    4142SRCS    = $(SRCFILES)

  • trunk/MagicSoft/Mars/mfbase/MFilterList.cc

    r3573 r3666  
    188188
    189189    if (fFilters.FindObject(name))
    190         *fLog << warn << "MFilterList::AddToList - '" << name << "' exists in List already..." << endl;
     190        *fLog << inf << "MFilterList::AddToList - '" << name << "' exists in List already..." << endl;
    191191
    192192    *fLog << inf << "Adding " << name << " to " << GetName() << "... " << flush;

  • trunk/MagicSoft/Mars/mgeom/MGeomCam.cc

    r3547 r3666  
    5353
    5454#include <TClass.h>     // IsA()->New()
     55#include <TVector2.h>   // TVector2
    5556
    5657#include "MLog.h"
    … …  
    9293
    9394    for (UInt_t i=0; i<npix; i++)
    94       fPixels[i] = new MGeomPix;
     95        fPixels[i] = new MGeomPix;
    9596
    9697    SetReadyToSave();
    … …  
    195196void MGeomCam::CalcMaxRadius()
    196197{
    197 
    198   fMaxRadius.Set(fNumAreas+1);
    199   fMinRadius.Set(fNumAreas+1); 
    200 
    201   for (Int_t i=0; i<fNumAreas+1; i++)
    202     {
    203       fMaxRadius[i] = 0.;
    204       fMinRadius[i] = FLT_MAX;
    205     }
    206  
    207   for (UInt_t i=0; i<fNumPixels; i++)
    208     {
    209 
    210       const MGeomPix &pix = (*this)[i];
    211 
    212       const UInt_t  s = pix.GetAidx();     
    213       const Float_t x = pix.GetX();
    214       const Float_t y = pix.GetY();
    215       const Float_t d = pix.GetD();
    216 
    217       const Float_t r = TMath::Hypot(x, y);
    218 
    219       const Float_t maxr = r + d;
    220       const Float_t minr = r>d ? r-d : 0;
    221      
    222       if (maxr>fMaxRadius[s+1])
    223         fMaxRadius[s+1] = maxr;
    224 
    225       if (minr<fMinRadius[s+1])
    226         fMinRadius[s+1] = minr;
    227 
    228       if (minr<fMinRadius[0])
    229         fMinRadius[0] = minr;
    230 
    231       if (maxr>fMaxRadius[0])
    232         fMaxRadius[0] = maxr;
    233 
    234     }
    235 }
    236 
     198    fMaxRadius.Set(fNumAreas+1);
     199    fMinRadius.Set(fNumAreas+1);
     200
     201    for (Int_t i=0; i<fNumAreas+1; i++)
     202    {
     203        fMaxRadius[i] = 0.;
     204        fMinRadius[i] = FLT_MAX;
     205    }
     206
     207    for (UInt_t i=0; i<fNumPixels; i++)
     208    {
     209        const MGeomPix &pix = (*this)[i];
     210
     211        const UInt_t  s = pix.GetAidx();
     212        const Float_t x = pix.GetX();
     213        const Float_t y = pix.GetY();
     214        const Float_t d = pix.GetD();
     215
     216        const Float_t r = TMath::Hypot(x, y);
     217
     218        const Float_t maxr = r + d;
     219        const Float_t minr = r>d ? r-d : 0;
     220
     221        if (maxr>fMaxRadius[s+1])
     222            fMaxRadius[s+1] = maxr;
     223
     224        if (minr<fMinRadius[s+1])
     225            fMinRadius[s+1] = minr;
     226
     227        if (minr<fMinRadius[0])
     228            fMinRadius[0] = minr;
     229
     230        if (maxr>fMaxRadius[0])
     231            fMaxRadius[0] = maxr;
     232    }
     233}
     234
     235// --------------------------------------------------------------------------
    237236//
    238237// Have to call the radii of the subcameras starting to count from 1
    … …  
    240239Float_t MGeomCam::GetMaxRadius(const Int_t i) const
    241240{
    242   if (i==-1) return fMaxRadius[0];
    243   return i>fNumAreas ? -1 : fMaxRadius[i+1];
    244 }
    245 
     241    return i<-1 || i>fNumAreas ? -1 : fMaxRadius[i+1];
     242}
     243
     244// --------------------------------------------------------------------------
    246245//
    247246// Have to call the radii of the subcameras starting to count from 1
    … …  
    249248Float_t MGeomCam::GetMinRadius(const Int_t i) const
    250249{
    251   if (i==-1) return fMinRadius[0];
    252   return i>fNumAreas ? -1 : fMinRadius[i+1];
     250    return i<-1 || i>fNumAreas ? -1 : fMinRadius[i+1];
    253251}
    254252
    … …  
    302300    return (TObject*)IsA()->New();
    303301}
     302
     303// --------------------------------------------------------------------------
     304//
     305//  Return the pixel index corresponding to the coordinates given in x, y.
     306//  The coordinates are given in pixel units (millimeters)
     307//  If no pixel exists return -1;
     308//
     309Int_t MGeomCam::GetPixelIdxXY(Float_t x, Float_t y) const
     310{
     311    for (unsigned int i=0; i<fNumPixels; i++)
     312        if ((*this)[i].IsInside(x, y))
     313            return i;
     314
     315    return -1;
     316}
     317
     318Int_t MGeomCam::GetPixelIdx(const TVector2 &v) const
     319{
     320    return GetPixelIdxXY(v.X(), v.Y());
     321}
     322
     323Int_t MGeomCam::GetPixelIdxDeg(const TVector2 &v) const
     324{
     325    return GetPixelIdxXYdeg(v.X(), v.Y());
     326}
     327
    304328/*
    305329void MGeomCam::Streamer(TBuffer &R__b)

  • trunk/MagicSoft/Mars/mgeom/MGeomCam.h

    r3547 r3666  
    1212#endif
    1313
     14class TVector2;
    1415class MGeomPix;
    1516
    … …  
    7374
    7475    MGeomPix &operator[](Int_t i);
    75     MGeomPix &operator[](Int_t i)          const;
     76    MGeomPix &operator[](Int_t i) const;
     77
     78    Int_t GetPixelIdx(const TVector2 &v) const;
     79    Int_t GetPixelIdxDeg(const TVector2 &v) const;
     80    Int_t GetPixelIdxXY(Float_t x, Float_t y) const;
     81    Int_t GetPixelIdxXYdeg(Float_t x, Float_t y) const
     82    {
     83        return GetPixelIdxXY(x/fConvMm2Deg, y/fConvMm2Deg);
     84    }
    7685
    7786    virtual void Print(Option_t *opt=NULL)   const;

  • trunk/MagicSoft/Mars/mgeom/MGeomPix.cc

    r3392 r3666  
    133133    *fLog << "d= " << fD << "mm  A= " << fA << "mm²" << endl;
    134134}
     135
     136// ------------------------------------------------------------------------
     137//
     138// compute the distance of a point (px,py) to the Hexagon center in
     139// MGeomPix coordinates. Return kTRUE if inside.
     140//
     141Bool_t MGeomPix::IsInside(Float_t px, Float_t py) const
     142{
     143    //
     144    //  compute the distance of the Point to the center of the Hexagon
     145    //
     146    const Double_t dx = px-fX;
     147
     148    //
     149    // Now check if point is outside of hexagon; just check x coordinate
     150    // in three coordinate systems: the default one, in which two sides of
     151    // the hexagon are paralel to the y axis (see camera displays) and two
     152    // more, rotated with respect to that one by +- 60 degrees.
     153    //
     154    if (TMath::Abs(dx)*2>fD)
     155        return kFALSE;
     156
     157    const Double_t dy = py-fY;
     158
     159    const static Double_t cos60 = TMath::Cos(60/kRad2Deg);
     160    const static Double_t sin60 = TMath::Sin(60/kRad2Deg);
     161
     162    const Double_t dx2 = dx*cos60 + dy*sin60;
     163    if  (TMath::Abs(dx2)*2>fD)
     164        return kFALSE;
     165
     166    const Double_t dx3 = dx*cos60 - dy*sin60;
     167    if (TMath::Abs(dx3)*2>fD)
     168        return kFALSE;
     169
     170    return kTRUE;
     171}

  • trunk/MagicSoft/Mars/mgeom/MGeomPix.h

    r3392 r3666  
    5656    Bool_t IsInOuterRing() const     { return TestBit(kIsInOuterRing); }
    5757
     58    Bool_t IsInside(Float_t px, Float_t py) const;
     59
    5860    /*
    5961     //

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

    r3597 r3666  
    9999//    also different algorithms like (Li/Ma)
    100100//  - implement fit for best alpha distribution -- online (Paint)
     101//  - currently a constant pointing position is assumed in Fill()
     102//  - the center of rotation need not to be the camera center
    101103//
    102104//////////////////////////////////////////////////////////////////////////////
    … …  
    117119#include "MObservatory.h"
    118120#include "MPointingPos.h"
     121#include "MAstroCatalog.h"
    119122#include "MAstroSky2Local.h"
    120123#include "MStatusDisplay.h"
    … …  
    136139//
    137140MHFalseSource::MHFalseSource(const char *name, const char *title)
    138     : fMm2Deg(-1), fUseMmScale(kTRUE), fAlphaCut(12.5), fBgMean(55)
     141    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1),
     142      fAlphaCut(12.5), fBgMean(55), fDistMin(-1), fDistMax(-1)
    139143{
    140144    //
    … …  
    186190// --------------------------------------------------------------------------
    187191//
    188 // Use this function to setup your own conversion factor between degrees
    189 // and millimeters. The conversion factor should be the one calculated in
    190 // MGeomCam. Use this function with Caution: You could create wrong values
    191 // by setting up your own scale factor.
    192 //
    193 void MHFalseSource::SetMm2Deg(Float_t mmdeg)
    194 {
    195     if (mmdeg<0)
    196     {
    197         *fLog << warn << dbginf << "Warning - Conversion factor < 0 - nonsense. Ignored." << endl;
    198         return;
    199     }
    200 
    201     if (fMm2Deg>=0)
    202         *fLog << warn << dbginf << "Warning - Conversion factor already set. Overwriting" << endl;
    203 
    204     fMm2Deg = mmdeg;
    205 }
    206 
    207 // --------------------------------------------------------------------------
    208 //
    209 // With this function you can convert the histogram ('on the fly') between
    210 // degrees and millimeters.
    211 //
    212 void MHFalseSource::SetMmScale(Bool_t mmscale)
    213 {
    214     if (fUseMmScale == mmscale)
    215         return;
    216 
    217     if (fMm2Deg<0)
    218     {
    219         *fLog << warn << dbginf << "Warning - Sorry, no conversion factor for conversion available." << endl;
    220         return;
    221     }
    222 
    223     if (fUseMmScale)
    224     {
    225         fHist.SetXTitle("x [mm]");
    226         fHist.SetYTitle("y [mm]");
    227 
    228         fHist.Scale(1./fMm2Deg);
    229     }
    230     else
    231     {
    232         fHist.SetXTitle("x [\\circ]");
    233         fHist.SetYTitle("y [\\circ]");
    234 
    235         fHist.Scale(1./fMm2Deg);
    236     }
    237 
    238     fUseMmScale = mmscale;
    239 }
    240 
    241 // --------------------------------------------------------------------------
    242 //
    243192// Calculate Significance as
    244193// significance = (s-b)/sqrt(s+k*k*b) mit k=s/b
    … …  
    249198Double_t MHFalseSource::Significance(Double_t s, Double_t b)
    250199{
    251     return MMath::Significance(s, b);
     200    return MMath::SignificanceSym(s, b);
    252201}
    253202
    … …  
    263212Double_t MHFalseSource::SignificanceLiMa(Double_t s, Double_t b, Double_t alpha)
    264213{
    265     const Double_t lima = MMath::SignificanceLiMa(s, b);
    266     return lima<0 ? 0 : lima;
     214    return MMath::SignificanceLiMaSigned(s, b);
     215    /*
     216     const Double_t lima = MMath::SignificanceLiMa(s, b);
     217     return lima<0 ? 0 : lima;
     218     */
    267219}
    268220
    … …  
    277229{
    278230    const MGeomCam *geom = (MGeomCam*)plist->FindObject("MGeomCam");
    279     if (geom)
     231    if (!geom)
    280232    {
    281         fMm2Deg = geom->GetConvMm2Deg();
    282         fUseMmScale = kFALSE;
    283 
    284         fHist.SetXTitle("x [\\circ]");
    285         fHist.SetYTitle("y [\\circ]");
     233        *fLog << err << "MGeomCam not found... aborting." << endl;
     234        return kFALSE;
    286235    }
     236
     237    fMm2Deg = geom->GetConvMm2Deg();
    287238
    288239    MBinning binsa;
    … …  
    292243    if (!bins)
    293244    {
    294         Float_t r = geom ? geom->GetMaxRadius() : 600;
    295         r /= 3;
    296         if (!fUseMmScale)
    297             r *= fMm2Deg;
     245        const Float_t r = (geom ? geom->GetMaxRadius()/3 : 200)*fMm2Deg;
    298246
    299247        MBinning b;
    … …  
    314262    fObservatory = (MObservatory*)plist->FindObject(AddSerialNumber("MObservatory"));
    315263    if (!fObservatory)
    316         *fLog << err << "MObservatory not found...  no derotation." << endl;
    317 
     264        *fLog << warn << "MObservatory not found... no derotation." << endl;
     265
     266    // FIXME: Because the pointing position could change we must check
     267    // for the current pointing position and add a offset in the
     268    // Fill function!
     269    fRa  = fPointPos->GetRa();
     270    fDec = fPointPos->GetDec();
    318271
    319272    return kTRUE;
    … …  
    326279Bool_t MHFalseSource::Fill(const MParContainer *par, const Stat_t w)
    327280{
    328     MHillas *hil = (MHillas*)par;
    329 
     281    const MHillas *hil = dynamic_cast<const MHillas*>(par);
     282    if (!hil)
     283    {
     284        *fLog << err << "MHFalseSource::Fill: No container specified!" << endl;
     285        return kFALSE;
     286    }
     287
     288    // Get max radius...
    330289    const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
    331290
    … …  
    337296    //    rho = fPointPos->RotationAngle(*fObservatory);
    338297
     298    // Create necessary containers for calculation
    339299    MSrcPosCam src;
    340300    MHillasSrc hsrc;
    341301    hsrc.SetSrcPos(&src);
    342302
     303    // Get number of bins and bin-centers
    343304    const Int_t nx = fHist.GetNbinsX();
    344305    const Int_t ny = fHist.GetNbinsY();
    … …  
    352313        for (int iy=0; iy<ny; iy++)
    353314        {
     315            // check distance... to get a circle plot
    354316            if (TMath::Hypot(cx[ix], cy[iy])>maxr)
    355317                continue;
    356318
     319            // rotate center of bin
     320            // precalculation of sin/cos doesn't accelerate
    357321            TVector2 v(cx[ix], cy[iy]);
    358322            if (rho!=0)
    359323                v=v.Rotate(rho);
    360324
    361             if (!fUseMmScale)
    362                 v *= 1./fMm2Deg;
    363 
     325            // convert degrees to millimeters
     326            v *= 1./fMm2Deg;
    364327            src.SetXY(v);
    365328
     329            // Source dependant hillas parameters
    366330            if (!hsrc.Calc(hil))
    367331            {
    … …  
    370334            }
    371335
     336            // FIXME: This should be replaced by an external MFilter
     337            //        and/or MTaskList
     338            // Source dependant distance cut
     339            if (fDistMin>0 && hsrc.GetDist()*fMm2Deg<fDistMin)
     340                continue;
     341
     342            if (fDistMax>0 && hil->GetLength()>fDistMax*hsrc.GetDist())
     343                continue;
     344
     345            // Fill histogram
    372346            const Double_t alpha = hsrc.GetAlpha();
    373 
    374347            fHist.Fill(cx[ix], cy[iy], TMath::Abs(alpha), w);
    375348        }
    … …  
    446419// --------------------------------------------------------------------------
    447420//
    448 // Update the projections
     421// Create projection for all data, taking sum of bin contents of
     422// range (0, 90) - corresponding to the number of entries in this slice.
     423//
     424void MHFalseSource::ProjectAll(TH2D *h3)
     425{
     426    h3->SetTitle("Number of entries");
     427
     428    // Get projection for range
     429    TH2D *p = (TH2D*)fHist.Project3D("yx_all");
     430
     431    // Move contents from projection to h3
     432    h3->Reset();
     433    h3->Add(p);
     434    delete p;
     435
     436    // Set Minimum as minimum value Greater Than 0
     437    h3->SetMinimum(GetMinimumGT(*h3));
     438}
     439
     440// --------------------------------------------------------------------------
     441//
     442// Update the projections and paint them
    449443//
    450444void MHFalseSource::Paint(Option_t *opt)
    451445{
    452     // sigma = (s-b)/sqrt(s+k*k*b) mit k=s/b
    453 
     446    // Set pretty color palette
    454447    gStyle->SetPalette(1, 0);
    455448
    … …  
    457450
    458451    TH1D* h1;
     452    TH2D* h0;
    459453    TH2D* h2;
    460454    TH2D* h3;
    461455    TH2D* h4;
    462456
    463     padsave->cd(3);
     457    // Update projection of all-events
     458    padsave->GetPad(1)->cd(3);
     459    if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all")))
     460        ProjectAll(h0);
     461
     462    // Update projection of off-events
     463    padsave->GetPad(2)->cd(1);
     464    if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
     465        ProjectOff(h2);
     466
     467    // Update projection of on-events
     468    padsave->GetPad(2)->cd(2);
    464469    if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
    465470        ProjectOn(h3);
    466471
    467     padsave->cd(4);
    468     if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
    469         ProjectOff(h2);
    470 
    471     padsave->cd(2);
     472    // Update projection of significance
     473    padsave->GetPad(2)->cd(3);
    472474    if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig")))
    473475    {
    474476        const Int_t nx = h4->GetXaxis()->GetNbins();
    475477        const Int_t ny = h4->GetYaxis()->GetNbins();
    476         const Int_t nr = nx*nx + ny*ny;
     478        //const Int_t nr = nx*nx + ny*ny;
    477479
    478480        Int_t maxx=nx/2;
    … …  
    493495                h4->SetBinContent(n, sig);
    494496
    495                 if (sig>h4->GetBinContent(max) && sig>0 && ix*ix+iy*iy<nr*nr/9)
     497                if (sig>h4->GetBinContent(max) && sig>0/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
    496498                {
    497499                    max = n;
    … …  
    501503            }
    502504
    503         padsave->cd(1);
     505        // Update projection of 'the best alpha-plot'
     506        padsave->GetPad(1)->cd(1);
    504507        if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0)
    505508        {
    … …  
    520523// --------------------------------------------------------------------------
    521524//
     525// Get the MAstroCatalog corresponding to fRa, fDec. The limiting magnitude
     526// is set to 9, while the fov is equal to the current fov of the false
     527// source plot.
     528//
     529TObject *MHFalseSource::GetCatalog()
     530{
     531    const Double_t maxr = 0.98*TMath::Abs(fHist.GetBinCenter(1));
     532
     533    // Create catalog...
     534    MAstroCatalog stars;
     535    stars.SetLimMag(9);
     536    stars.SetGuiActive(kFALSE);
     537    stars.SetRadiusFOV(maxr);
     538    stars.SetRaDec(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
     539    stars.ReadBSC("bsc5.dat");
     540
     541    TObject *o = (MAstroCatalog*)stars.Clone();
     542    o->SetBit(kCanDelete);
     543
     544    return o;
     545}
     546
     547// --------------------------------------------------------------------------
     548//
    522549// Draw the histogram
    523550//
    … …  
    529556    AppendPad("");
    530557
    531     pad->Divide(2, 2);
     558    pad->Divide(1, 2, 0, 0.03);
     559
     560    TObject *catalog = GetCatalog();
    532561
    533562    // draw the 2D histogram Sigmabar versus Theta
    534563    pad->cd(1);
    535564    gPad->SetBorderMode(0);
     565    gPad->Divide(3, 1);
     566    delete pad->GetPad(1)->GetPad(2);
     567
     568    pad->GetPad(1)->cd(3);
     569    gPad->SetBorderMode(0);
     570    TH1 *h0 = fHist.Project3D("yx_all");
     571    h0->SetDirectory(NULL);
     572    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
     573    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
     574    h0->Draw("colz");
     575    h0->SetBit(kCanDelete);
     576    catalog->Draw("mirror same");
     577
     578    pad->GetPad(1)->cd(1);
     579    gPad->SetBorderMode(0);
     580
    536581    TH1 *h1 = fHist.ProjectionZ("Alpha");
    537582    h1->SetDirectory(NULL);
    … …  
    542587    h1->SetBit(kCanDelete);
    543588
    544     pad->cd(4);
     589    pad->cd(2);
     590    gPad->SetBorderMode(0);
     591    gPad->Divide(3, 1);
     592
     593    pad = gPad;
     594
     595    pad->cd(1);
    545596    gPad->SetBorderMode(0);
    546597    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
    … …  
    551602    h2->Draw("colz");
    552603    h2->SetBit(kCanDelete);
    553 
    554     pad->cd(3);
     604    catalog->Draw("mirror same");
     605
     606    pad->cd(2);
    555607    gPad->SetBorderMode(0);
    556608    fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
    … …  
    562614    h3->Draw("colz");
    563615    h3->SetBit(kCanDelete);
    564 
    565     pad->cd(2);
     616    catalog->Draw("mirror same");
     617
     618    pad->cd(3);
    566619    gPad->SetBorderMode(0);
    567620    fHist.GetZaxis()->SetRange(0,0);
    … …  
    574627    h4->Draw("colz");
    575628    h4->SetBit(kCanDelete);
     629    catalog->Draw("mirror same");
    576630}
    577631
    … …  
    596650    TVirtualPad *padsave = gPad;
    597651    padsave->Modified();
    598     padsave->cd(1);
     652    padsave->GetPad(1)->cd(1);
    599653    gPad->Modified();
    600     padsave->cd(2);
     654    padsave->GetPad(1)->cd(3);
    601655    gPad->Modified();
    602     padsave->cd(3);
     656    padsave->GetPad(2)->cd(1);
    603657    gPad->Modified();
    604     padsave->cd(4);
     658    padsave->GetPad(2)->cd(2);
     659    gPad->Modified();
     660    padsave->GetPad(2)->cd(3);
    605661    gPad->Modified();
    606662    gPad->cd();
    … …  
    639695void MHFalseSource::FitSignificance(Float_t sigint, Float_t sigmax, Float_t bgmin, Float_t bgmax, Byte_t polynom)
    640696{
     697    TObject *catalog = GetCatalog();
     698
    641699    TH1D h0a("A",          "", 50,   0, 4000);
    642700    TH1D h4a("chisq1",     "", 50,   0,   35);
    … …  
    686744    // Implementing the function yourself is only about 5% faster
    687745    TF1 func("", Form("gaus(0) + pol%d(3)", polynom), 0, 90);
     746    //TF1 func("", Form("[0]*(TMath::Gaus(x, [1], [2])+TMath::Gaus(x, -[1], [2]))+pol%d(3)", polynom), 0, 90);
    688747    TArrayD maxpar(func.GetNpar());
    689748
    690     /*
    691      func.SetParName(0, "A");
    692      func.SetParName(1, "mu");
    693      func.SetParName(2, "sigma");
    694      func.SetParName(3, "a");
    695      func.SetParName(4, "b");
    696      func.SetParName(5, "c");
     749    /*  func.SetParName(0, "A");
     750     *  func.SetParName(1, "mu");
     751     *  func.SetParName(2, "sigma");
    697752    */
    698753
    … …  
    704759    const Int_t nx = hist->GetXaxis()->GetNbins();
    705760    const Int_t ny = hist->GetYaxis()->GetNbins();
    706     const Int_t nr = nx*nx+ny*ny;
     761    //const Int_t nr = nx*nx+ny*ny;
    707762
    708763    Double_t maxalpha0=0;
    … …  
    748803
    749804            h->Fit(&func, "N0Q", "", bgmin, bgmax);
    750             //*fLog << dbg << ix << "/" << iy << ":  " << func.GetParameter(3) << "    " << func.GetParameter(4) << endl;
    751805
    752806            h4a.Fill(func.GetChisquare());
    … …  
    774828
    775829            h->Fit(&func, "N0Q", "", 0, sigmax);
    776             //*fLog << dbg << "     " << func.GetParameter(0) << "    " << func.GetParameter(1) << "    " << func.GetParameter(2) << endl;
    777830
    778831            TArrayD p(func.GetNpar(), func.GetParameters());
    … …  
    820873                h6.Fill(sig);
    821874
    822             if (sig>maxs && ix*ix+iy*iy<nr*nr/9)
     875            if (sig>maxs/* && (ix-nx/2)*(ix-nx/2)+(iy-ny/2)*(iy-ny/2)<nr*nr/9*/)
    823876            {
    824877                maxs = sig;
    … …  
    838891
    839892    clk.Stop();
    840     clk.Print();
     893    clk.Print("m");
    841894
    842895    TCanvas *c=new TCanvas;
    … …  
    849902    hists->Draw("colz");
    850903    hists->SetBit(kCanDelete);
     904    catalog->Draw("mirror same");
    851905    c->cd(2);
    852906    gPad->SetBorderMode(0);
    853907    hist->Draw("colz");
    854908    hist->SetBit(kCanDelete);
     909    catalog->Draw("mirror same");
    855910    c->cd(3);
    856911    gPad->SetBorderMode(0);
    857912    histb->Draw("colz");
    858913    histb->SetBit(kCanDelete);
     914    catalog->Draw("mirror same");
    859915    c->cd(4);
    860916    gPad->Divide(1,3, 0, 0);
    … …  
    942998            const Double_t b = f2.Integral(0, (float)i)/w;
    943999
    944             const Double_t sig = Significance(s, b);
     1000            const Double_t sig = SignificanceLiMa(s, b);
    9451001
    9461002            g->SetPoint(g->GetN(), i, sig);

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

    r3597 r3666  
    1212class TH2D;
    1313
    14 class MHillasSrc;
    15 class MEnergyEst;
    1614class MParList;
    1715class MTime;
    … …  
    2220{
    2321private:
    24     MTime        *fTime;        //! container to take the event time from
    25     MPointingPos *fPointPos;    //! container to take pointing position from
    26     MObservatory *fObservatory; //! conteiner to take observatory location from
     22    MTime         *fTime;        //! container to take the event time from
     23    MPointingPos  *fPointPos;    //! container to take pointing position from
     24    MObservatory  *fObservatory; //! conteiner to take observatory location from
    2725
    28     Float_t fMm2Deg;            // conversion factor for display in degrees
    29     Bool_t  fUseMmScale;        // which scale to use?
     26    Float_t fMm2Deg;             // conversion factor for display in degrees
    3027
    31     Float_t fAlphaCut;          // Alpha cut
    32     Float_t fBgMean;            // Background mean
     28    Float_t fAlphaCut;           // Alpha cut
     29    Float_t fBgMean;             // Background mean
    3330
    34     TH3D    fHist;              // Alpha vs. x and y
     31    Float_t fDistMin;            // Min dist
     32    Float_t fDistMax;            // Maximum distance in percent of dist
     33
     34    TH3D    fHist;               // Alpha vs. x and y
     35
     36    Double_t fRa;
     37    Double_t fDec;
    3538
    3639    Int_t DistancetoPrimitive(Int_t px, Int_t py);
    … …  
    3942    void ProjectOff(TH2D *h);
    4043    void ProjectOn(TH2D *h);
     44    void ProjectAll(TH2D *h);
     45
     46    TObject *GetCatalog();
    4147
    4248public:
    … …  
    4652    Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    4753
    48     void SetMmScale(Bool_t mmscale=kTRUE);
    49     void SetMm2Deg(Float_t mmdeg);
    50 
    5154    TH1 *GetHistByName(const TString name) { return &fHist; }
    5255
    5356    void FitSignificance(Float_t sigint=15, Float_t sigmax=70, Float_t bgmin=40, Float_t bgmax=70, Byte_t polynom=1); //*MENU*
    5457    void FitSignificanceStd() { FitSignificance(); } //*MENU*
     58
     59    void SetDistMin(Float_t dist)  { fDistMin = dist;  } // Absolute minimum distance
     60    void SetDistMax(Float_t ratio) { fDistMax = ratio; } // Maximum ratio between length/dist
    5561
    5662    void SetAlphaCut(Float_t alpha); //*MENU*

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

    r3594 r3666  
    3434// For a given a shower, a series of points along its main axis are filled
    3535// into the 2-dim histogram (x, y) of the camera plane.
    36 // 
     36//
    3737// Before filling a point (x, y) into the histogram it is
    3838//        - shifted by (xSrc, ySrc)   (the expected source position)
    3939//        - and rotated in order to compensate the rotation of the
    4040//          sky image in the camera
    41 //         
    42 //         
     41//
     42// To be done:
     43//  - simplification of the algorithm, by doing all translation before
     44//    calculation of the line
     45//  - the center of rotation need not to be the camera center
     46//
    4347/////////////////////////////////////////////////////////////////////////////
    4448#include "MHStarMap.h"

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

    r2624 r3666  
    4747// fConc1    removed
    4848//
     49// Version 3:
     50// ----------
     51// fMaxDist  added. Distance between center and most distant used pixel
     52//
    4953//
    5054// WARNING: Before you can use fAsym, fM3Long and fM3Trans you must
    … …  
    96100    fM3Long  =  0;
    97101    fM3Trans =  0;
     102
     103    fMaxDist = -1;
    98104}
    99105
    … …  
    109115    *fLog << " - 3.Moment Long  [mm]  = " << fM3Long  << endl;
    110116    *fLog << " - 3.Moment Trans [mm]  = " << fM3Trans << endl;
     117    *fLog << " - Max.Dist       [mm]  = " << fMaxDist << endl;
    111118}
    112119
    … …  
    123130    *fLog << " - 3.Moment Long  [deg] = " << fM3Long *geom.GetConvMm2Deg() << endl;
    124131    *fLog << " - 3.Moment Trans [deg] = " << fM3Trans*geom.GetConvMm2Deg() << endl;
     132    *fLog << " - Max.Dist       [deg] = " << fMaxDist*geom.GetConvMm2Deg() << endl;
    125133}
    126134
    … …  
    146154    const UInt_t npixevt = evt.GetNumPixels();
    147155
    148     Int_t maxpixid = 0;
    149     Float_t maxpix = 0;
     156    Int_t maxpixid  = 0;
     157    Float_t maxpix  = 0;
     158    Float_t maxdist = 0;
    150159
    151160    for (UInt_t i=0; i<npixevt; i++)
    … …  
    160169        const Double_t dx = gpix.GetX() - hil.GetMeanX();      // [mm]
    161170        const Double_t dy = gpix.GetY() - hil.GetMeanY();      // [mm]
     171
     172        const Double_t dist = dx*dx+dy*dy;
     173        if (dist>maxdist)
     174            maxdist=dist;                                      // [mm^2]
    162175
    163176        Double_t nphot = pix.GetNumPhotons();                  // [1]
    … …  
    197210    fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3);      // [mm]
    198211
     212    fMaxDist = TMath::Sqrt(maxdist);                           // [mm]
     213
    199214    SetReadyToSave();
    200215

  • trunk/MagicSoft/Mars/mimage/MHillasExt.h

    r2624 r3666  
    2020    Float_t fM3Trans; // [mm] 3rd moment (e-weighted) along minor axis
    2121
     22    Float_t fMaxDist; // Distance between center and most distant used pixel
     23
    2224public:
    2325    MHillasExt(const char *name=NULL, const char *title=NULL);
    … …  
    2931    Float_t GetM3Trans() const { return fM3Trans; }
    3032
     33    Float_t GetMaxDist() const { return fMaxDist; }
     34
    3135    Int_t Calc(const MGeomCam &geom, const MCerPhotEvt &pix, const MHillas &hil);
    3236
    … …  
    3640    void Set(const TArrayF &arr);
    3741
    38     ClassDef(MHillasExt, 2) // Storage Container for extended Hillas Parameter
     42    ClassDef(MHillasExt, 3) // Storage Container for extended Hillas Parameter
    3943};
    4044#endif

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

    r2470 r3666  
    114114Int_t MHillasSrcCalc::Process()
    115115{
    116 
    117116    if (!fHillasSrc->Calc(fHillas))
    118117    {

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

    r3574 r3666  
    4444//  - added fInnerLeakage2
    4545//  - added fInnerSize
     46//  - added fUsedArea
     47//  - added fCoreArea
     48//
    4649//
    4750/////////////////////////////////////////////////////////////////////////////
    … …  
    9295    fNumCorePixels = -1;
    9396
     97    fUsedArea = -1;
     98    fCoreArea = -1;
     99
    94100    fNumSaturatedPixels = -1;
    95101    fNumHGSaturatedPixels = -1;
    … …  
    105111    fNumUsedPixels = 0;
    106112    fNumCorePixels = 0;
     113
     114    fUsedArea = 0;
     115    fCoreArea = 0;
    107116
    108117    fNumSaturatedPixels = 0;
    … …  
    136145            continue;
    137146
     147        // Get geometry of pixel
     148        const Int_t pixid = pix.GetPixId();
     149        const MGeomPix &gpix = geom[pixid];
     150
    138151        // count used and core pixels
    139152        if (pix.IsPixelCore())
     153        {
    140154            fNumCorePixels++;
     155            fCoreArea += gpix.GetA();
     156        }
    141157
    142158        // count used pixels
    143159        fNumUsedPixels++;
    144 
    145         const Int_t pixid = pix.GetPixId();
    146 
    147         const MGeomPix &gpix = geom[pixid];
     160        fUsedArea += gpix.GetA();
    148161
    149162        Double_t nphot = pix.GetNumPhotons();
    … …  
    165178        // count inner pixels: To dependent on MAGIC Camera --> FIXME
    166179
    167         if (pixid<397){
    168           fInnerSize += nphot;
    169           if(pixid>270){
    170             edgepixin2 += nphot;
    171             if(pixid>330)
    172               edgepixin1 += nphot;
    173           }
    174         }
     180        if (pixid<397){
     181            fInnerSize += nphot;
     182            if(pixid>270){
     183                edgepixin2 += nphot;
     184                if(pixid>330)
     185                    edgepixin1 += nphot;
     186            }
     187        }
    175188
    176189        // Compute Concetration 1 -2
    177190
    178191        if (nphot>maxpix1)
    179           {
     192        {
    180193            maxpix2  = maxpix1;
    181194            maxpix1  = nphot;                            // [1]
    182195            continue;                                    // [1]
    183           }
    184        
     196        }
     197
    185198        if (nphot>maxpix2)
    186           maxpix2 = nphot;                             // [1]
     199            maxpix2 = nphot;                             // [1]
    187200    }
    188201   

  • trunk/MagicSoft/Mars/mimage/MNewImagePar.h

    r3574 r3666  
    1313{
    1414private:
    15     Float_t fLeakage1;           // (photons in most outer ring of pixels) over fSize
    16     Float_t fLeakage2;           // (photons in the 2 outer rings of pixels) over fSize
    17     Float_t fInnerLeakage1;      // (photons in most outer rings of inner pixels) over fInnerSize
    18     Float_t fInnerLeakage2;      // (photons in the 2 outer rings of inner pixels) over fInnerSize
    19     Float_t fInnerSize;          //
     15    Float_t fLeakage1;             // (photons in most outer ring of pixels) over fSize
     16    Float_t fLeakage2;             // (photons in the 2 outer rings of pixels) over fSize
     17    Float_t fInnerLeakage1;        // (photons in most outer rings of inner pixels) over fInnerSize
     18    Float_t fInnerLeakage2;        // (photons in the 2 outer rings of inner pixels) over fInnerSize
     19    Float_t fInnerSize;            //
    2020
    21     Float_t fConc;               // [ratio] concentration ratio: sum of the two highest pixels / fSize
    22     Float_t fConc1;              // [ratio] concentration ratio: sum of the highest pixel / fSize
     21    Float_t fConc;                 // [ratio] concentration ratio: sum of the two highest pixels / fSize
     22    Float_t fConc1;                // [ratio] concentration ratio: sum of the highest pixel / fSize
    2323
    24     Short_t fNumUsedPixels;      // Number of pixels which survived the image cleaning
    25     Short_t fNumCorePixels;      // number of core pixels
     24    Float_t fUsedArea;             // Area of pixels which survived the image cleaning
     25    Float_t fCoreArea;             // Area of core pixels
     26    Short_t fNumUsedPixels;        // Number of pixels which survived the image cleaning
     27    Short_t fNumCorePixels;        // number of core pixels
    2628    Short_t fNumHGSaturatedPixels; // number of pixels with saturating hi-gains
    2729    Short_t fNumSaturatedPixels;   // number of pixels with saturating lo-gains
    … …  
    4143    Short_t GetNumCorePixels() const { return fNumCorePixels; }
    4244
     45    Float_t GetNumUsedArea() const { return fUsedArea; }
     46    Float_t GetNumCoreArea() const { return fCoreArea; }
     47
    4348    Short_t GetNumSaturatedPixels() const { return fNumSaturatedPixels; }
    4449    Short_t GetNumHGSaturatedPixels() const { return fNumHGSaturatedPixels; }

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

    r3568 r3666  
    7373Double_t MPointingPos::RotationAngle(const MObservatory &o, const MTime &t) const
    7474{
    75     return MAstroSky2Local(t, o).RotationAngle(fRa*TMath::DegToRad()*15, fDec*TMath::DegToRad());
     75    return MAstroSky2Local(t, o).RotationAngle(GetRaRad(), GetDecRad());
    7676}

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

    r3568 r3666  
    1515    Double_t fAz;  // [deg] Azimuth
    1616
    17     Double_t fRa;  // [rad] Right ascension
    18     Double_t fHa;  // [rad] Hour angle
    19     Double_t fDec; // [rad] Declination
     17    Double_t fRa;  // [h]  Right ascension
     18    Double_t fHa;  // [h]  Hour angle
     19    Double_t fDec; // [deg] Declination
    2020
    2121public:
    22     MPointingPos()
     22    MPointingPos(const char *name=0, const char *title=0)
    2323    {
    24         fName  = "MPointingPos";
    25         fTitle = "Container storing the (corrected) telescope pointing position";
     24        fName  = name ? name   : "MPointingPos";
     25        fTitle = title ? title : "Container storing the (corrected) telescope pointing position";
    2626    }
    2727
    … …  
    3535    Double_t GetDec() const { return fDec; }
    3636
     37    Double_t GetRaRad() const  { return fRa*TMath::DegToRad()*15; }
     38    Double_t GetDecRad() const { return fDec*TMath::DegToRad(); }
     39
    3740    Double_t RotationAngle(const MObservatory &o) const;
    3841    Double_t RotationAngle(const MObservatory &o, const MTime &t) const;
     42    Double_t RotationAngle(const MObservatory &o, const MTime *t) const
     43    {
     44        return t ? RotationAngle(o, *t) : RotationAngle(o);
     45    }
    3946
    4047    ClassDef(MPointingPos, 1) //Container storing the (corrected) telescope pointing position

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

    r3592 r3666  
    4646// To be done:
    4747//   - wobble mode missing
     48//   - a switch between using sky-coordinates and time or local-coordinates
     49//     from MPointingPos for determin the rotation angle
     50//   - the center of rotation need not to be the camera center
    4851//
    4952//////////////////////////////////////////////////////////////////////////////
    … …  
    101104    {
    102105        *fLog << err << "MObservatory not found... aborting." << endl;
     106        return kFALSE;
     107    }
     108
     109    fTime = (MTime*)pList->FindObject("MTime");
     110    if (!fTime)
     111    {
     112        *fLog << err << "MTime not found... aborting." << endl;
    103113        return kFALSE;
    104114    }
    … …  
    166176        return kTRUE;
    167177
     178    // rotate the source position by the current rotation angle
     179    const Double_t rho = fPointPos->RotationAngle(*fObservatory, *fTime);
     180
    168181    // Define source position in the camera plain
    169182    TVector2 v(fX, fY);
    170 
    171     // rotate the source position by the current rotation angle
    172     const Double_t rho = fPointPos->RotationAngle(*fObservatory);
    173     v=v.Rotate(-rho);
     183    v=v.Rotate(rho);
    174184
    175185    // Convert coordinates into camera plain (mm)

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

    r3568 r3666  
    1010class MSrcPosCam;
    1111class MGeomCam;
     12class MTime;
    1213
    1314class MSrcPosCalc : public MTask
    … …  
    1819    MSrcPosCam   *fSrcPos;
    1920    MGeomCam     *fGeom;
     21    MTime        *fTime;
    2022
    2123    Double_t fR;    // Distance of source to a fitted star
    … …  
    3234    MSrcPosCalc(const char *name=NULL, const char *title=NULL);
    3335
     36    // Use is deprecated!
    3437    void SetOffset(Double_t r, Double_t drho)
    3538    {
Note: See TracChangeset for help on using the changeset viewer.