Changeset 4716 for trunk/MagicSoft


Ignore:
Timestamp:
08/23/04 17:49:53 (20 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
11 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r4714 r4716  
    7979   * mhist/MHEvent.[h,cc]:
    8080     - added new option for island index
     81     - added kEvtCleaningData
    8182
    8283   * mimage/MImgCleanStd.[h,cc]:
     
    8586     - added ReadEnv
    8687     - updated StreamPrimitive
     88     - added new cleaning method (probability cleaning)
    8789
    8890   * mimage/Makefile:
     
    142144   * msignal/MArrivalTime.h:
    143145     - removed operator()
     146     - added operator[] const
     147
     148   * manalysis/MCameraData.[h,cc]:
     149     - added algorithm for 'Probability cleaning'
     150
     151   * mbase/MMath.[h,cc]:
     152     - added GaussProb
     153
     154   * mjobs/MSequence.h:
     155     - added IsValid
    144156
    145157
  • trunk/MagicSoft/Mars/NEWS

    r4710 r4716  
    2525   - implemented new image parameters displaying the number of islands,
    2626     saturated hi-gain and lo-gain pixels (MImagePar, MHImagePar)
     27
     28   - event display in executable changed to support also calibrated files
     29     (done with MJCalibrateSignal)
     30
     31   - added a cleaning which takes signal height _and_ arrival time into
     32     account: probability cleaning (for more details see MImgCleanStd)
    2733
    2834
  • trunk/MagicSoft/Mars/manalysis/MCameraData.cc

    r4454 r4716  
    3434#include "MCameraData.h"
    3535
     36#include "MMath.h"
     37
     38#include "MLog.h"
     39#include "MLogManip.h"
     40
     41#include "MGeomCam.h"
     42#include "MGeomPix.h"
     43
     44#include "MPedPhotCam.h"
     45#include "MPedPhotPix.h"
     46
    3647#include "MCerPhotEvt.h"
    3748#include "MCerPhotPix.h"
    3849
    39 #include "MGeomCam.h"
    40 #include "MGeomPix.h"
    41 
    42 #include "MLog.h"
    43 #include "MLogManip.h"
    44 
    45 #include "MPedPhotCam.h"
    46 #include "MPedPhotPix.h"
    47 
    4850#include "MSigmabar.h"
     51#include "MArrivalTime.h"
    4952
    5053ClassImp(MCameraData);
     
    101104    const Int_t n = geom.GetNumPixels();
    102105
    103     fData.Set(n);
    104     fData.Reset();
    105 
    106     fValidity.Set(n);
    107     fValidity.Reset();
    108 
    109     const Int_t entries = evt.GetNumPixels();
    110 
    111     //
    112     // check the number of all pixels against the noise level and
    113     // set them to 'unused' state if necessary
    114     //
     106    // Reset arrays
     107    fData.Set(n);
     108    fData.Reset();
     109
     110    fValidity.Set(n);
     111    fValidity.Reset();
     112
     113    const Int_t entries = evt.GetNumPixels();
     114
     115    // calculate cleaning levels
    115116    for (Int_t i=0; i<entries; i++)
    116117    {
     
    149150    const Int_t n = geom.GetNumPixels();
    150151
    151     fData.Set(n);
    152     fData.Reset();
    153 
    154     fValidity.Set(n);
    155     fValidity.Reset();
    156 
    157     const Int_t entries   = evt.GetNumPixels();
     152    // reset arrays
     153    fData.Set(n);
     154    fData.Reset();
     155
     156    fValidity.Set(n);
     157    fValidity.Reset();
     158
     159    // check validity of rms with area index 0
    158160    const Float_t anoise0 = cam.GetArea(0).GetRms();
    159161    if (anoise0<=0)
    160162        return;
    161163
    162     //
    163     // check the number of all pixels against the noise level and
    164     // set them to 'unused' state if necessary
    165     //
     164    // calculate cleaning levels
     165    const Int_t entries = evt.GetNumPixels();
    166166    for (Int_t i=0; i<entries; i++)
    167167    {
     
    200200    const Int_t n = geom.GetNumPixels();
    201201
    202     fData.Set(n);
    203     fData.Reset();
    204 
    205     fValidity.Set(n);
    206     fValidity.Reset();
    207 
     202    // reset arrays
     203    fData.Set(n);
     204    fData.Reset();
     205
     206    fValidity.Set(n);
     207    fValidity.Reset();
     208
     209    // check validity of noise
    208210    if (noise<=0)
    209211        return;
    210212
    211     const Int_t entries = evt.GetNumPixels();
    212 
    213     //
    214     // check the number of all pixels against the noise level and
    215     // set them to 'unused' state if necessary
    216     //
     213    // calculate cleaning levels
     214    const Int_t entries = evt.GetNumPixels();
    217215    for (Int_t i=0; i<entries; i++)
    218216    {
     
    244242    const Int_t n = geom.GetNumPixels();
    245243
    246     fData.Set(n);
    247     fData.Reset();
    248 
    249     fValidity.Set(n);
    250     fValidity.Reset();
    251 
    252     const Int_t   entries = evt.GetNumPixels();
    253     const Float_t noise0  = cam.GetArea(0).GetRms();
     244    // reset arrays
     245    fData.Set(n);
     246    fData.Reset();
     247
     248    fValidity.Set(n);
     249    fValidity.Reset();
     250
     251    // check validity of rms with area index 0
     252    const Float_t noise0 = cam.GetArea(0).GetRms();
    254253    if (noise0<=0)
    255254        return;
    256255
    257     //
    258     // check the number of all pixels against the noise level and
    259     // set them to 'unused' state if necessary
    260     //
     256    // calculate cleaning levels
     257    const Int_t entries = evt.GetNumPixels();
    261258    for (Int_t i=0; i<entries; i++)
    262259    {
     
    280277// --------------------------------------------------------------------------
    281278//
     279// Function to calculate the cleaning level for all pixels in a given event.
     280// The level is the probability that the signal is a real signal (taking
     281// the signal height and the fluctuation of the background into account)
     282// times one minus the probability that the signal is a background
     283// fluctuation (calculated from the time spread of arrival times
     284// around the pixel with the highest signal)
     285//
     286// FIXME: Should the check noise<=0 be replaced by MBadPixels?
     287//
     288void MCameraData::CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam,
     289                                          const MGeomCam &geom,   const MArrivalTime *tcam)
     290{
     291    const Int_t n = geom.GetNumPixels();
     292
     293    // Reset arrays
     294    fData.Set(n);
     295    fData.Reset();
     296
     297    fValidity.Set(n);
     298    fValidity.Reset();
     299
     300    // check validity of noise
     301    const Float_t anoise0 = pcam.GetArea(0).GetRms();
     302    if (anoise0<=0)
     303        return;
     304
     305    const Int_t entries = evt.GetNumPixels();
     306
     307    // Find pixel with max signal
     308    Int_t maxidx = 0;
     309    if (tcam)
     310    {
     311        // Find pixel enty with maximum signal
     312        for (Int_t i=0; i<entries; i++)
     313        {
     314            const Double_t s0 = evt[i].GetNumPhotons()      * geom.GetPixRatio(i);
     315            const Double_t s1 = evt[maxidx].GetNumPhotons() * geom.GetPixRatio(maxidx);
     316            if (s0>s1)
     317                maxidx = i;
     318        }
     319
     320        // Get software index of pixel
     321        maxidx = evt[maxidx].GetPixId();
     322    }
     323
     324    const Double_t timemean = tcam ? (*tcam)[maxidx] : 0;
     325    const Double_t timerms  = 0.75; //[slices] rms time spread around highest pixel
     326
     327    // calculate cleaning levels
     328    for (Int_t i=0; i<entries; i++)
     329    {
     330        const MCerPhotPix &spix = evt[i];
     331
     332        const Int_t   idx = spix.GetPixId();
     333        const Float_t rms = pcam[idx].GetRms();
     334        if (rms<=0) // fData[idx]=0, fValidity[idx]=0
     335            continue;
     336
     337        fValidity[idx]=1;
     338
     339        // get signal and arrival time
     340        const UInt_t  aidx     = geom[idx].GetAidx();
     341        const Float_t ratio    = pcam.GetArea(aidx).GetRms()/anoise0;
     342
     343        const Double_t signal  = spix.GetNumPhotons() * geom.GetPixRatio(idx) * ratio / rms;
     344        const Double_t time    = tcam ? (*tcam)[idx] : 1;
     345
     346        // if signal<0 the probability is equal 0
     347        if (signal<0)
     348            continue;
     349
     350        // Just for convinience for easy readybility
     351        const Double_t means   = 0;
     352        const Double_t meant   = timemean;
     353
     354        const Double_t sigmas  = rms;
     355        const Double_t sigmat  = timerms;
     356
     357        // Get probabilities
     358        const Double_t psignal = MMath::GaussProb(signal, sigmas, means);
     359        const Double_t pbckgnd = MMath::GaussProb(time,   sigmat, meant);
     360
     361        // Set probability
     362        fData[idx]     = psignal*(1-pbckgnd);
     363        fValidity[idx] = 1;
     364
     365        // Make sure, that we don't run into trouble because of
     366        // a little numerical uncertanty
     367        if (fData[idx]>1)
     368            fData[idx]=1;
     369        if (fData[idx]<0)
     370            fData[idx]=0;
     371    }
     372}
     373
     374// --------------------------------------------------------------------------
     375//
    282376// Returns the contents of the pixel.
    283377//
  • trunk/MagicSoft/Mars/manalysis/MCameraData.h

    r4452 r4716  
    1919class MCerPhotEvt;
    2020class MPedPhotCam;
     21class MArrivalTime;
    2122
    2223class MCameraData : public MParContainer, public MCamEvent
     
    3839    void CalcCleaningLevel(const MCerPhotEvt &evt, Double_t noise,
    3940                           const MGeomCam &geom);
    40 
    4141    void CalcCleaningLevel2(const MCerPhotEvt &evt, const MPedPhotCam &fCam,
    42                            const MGeomCam &geom);
    43 
     42                            const MGeomCam &geom);
    4443    void CalcCleaningLevelDemocratic(const MCerPhotEvt &evt, const MPedPhotCam &cam,
    4544                                     const MGeomCam &geom);
     45    void CalcCleaningProbability(const MCerPhotEvt &evt, const MPedPhotCam &pcam,
     46                                 const MGeomCam &geom,   const MArrivalTime *tcam);
    4647
    4748    const TArrayD &GetData() const { return fData; }
  • trunk/MagicSoft/Mars/mbase/MMath.cc

    r3666 r4716  
    4747    const Double_t f = s+k*k*b;
    4848
    49     return f==0 ? 0 : (s-b)/TMath::Sqrt(f);
     49    return f==0 ? 0 : (s-b)/Sqrt(f);
    5050}
    5151
     
    8686        return -1;
    8787
    88     const Double_t l = s*TMath::Log(s/sum*(alpha+1)/alpha);
    89     const Double_t m = b*TMath::Log(b/sum*(alpha+1)      );
     88    const Double_t l = s*Log(s/sum*(alpha+1)/alpha);
     89    const Double_t m = b*Log(b/sum*(alpha+1)      );
    9090
    91     return l+m<0 ? -1 : TMath::Sqrt((l+m)*2);
     91    return l+m<0 ? -1 : Sqrt((l+m)*2);
    9292}
    9393
     
    104104        return 0;
    105105
    106     return TMath::Sign(sig, s-b);
     106    return Sign(sig, s-b);
    107107}
     108
     109// --------------------------------------------------------------------------
     110//
     111// Returns: 2/(sigma*sqrt(2))*integral[0,x](exp(-(x-mu)^2/(2*sigma^2)))
     112//
     113Double_t MMath::GaussProb(Double_t x, Double_t sigma, Double_t mean)
     114{
     115    static const Double_t sqrt2 = Sqrt(2.);
     116    return Erf((x-mean)/(sigma*sqrt2));
     117}
     118
  • trunk/MagicSoft/Mars/mbase/MMath.h

    r3666 r4716  
    99{
    1010public:
     11    static Double_t GaussProb(Double_t x, Double_t sigma, Double_t mean=0);
     12
    1113    static Double_t Significance(Double_t s, Double_t b);
    1214    static Double_t SignificanceSym(Double_t s, Double_t b);
  • trunk/MagicSoft/Mars/mhist/MHEvent.cc

    r4702 r4716  
    152152        fHist->SetYTitle("L");
    153153        break;
     154    case kEvtCleaningData:
     155        fHist->SetName("CleanData");
     156        fHist->SetYTitle("L");
     157        break;
    154158     case kEvtIdxMax:
    155159        fHist->SetName("Max Slice Idx");
     
    220224            fHist->SetLevels(lvl);
    221225        }
     226        break;
     227    case kEvtCleaningData:
     228        fHist->SetCamContent(*event, 0);
    222229        break;
    223230    case kEvtIdxMax:
  • trunk/MagicSoft/Mars/mhist/MHEvent.h

    r4702 r4716  
    2222        kEvtSignalRaw, kEvtSignalDensity, kEvtPedestal,
    2323        kEvtPedestalRMS, kEvtRelativeSignal, kEvtCleaningLevels,
     24        kEvtCleaningData,
    2425        kEvtIdxMax, kEvtArrTime, kEvtTrigPix, kEvtIslandIndex
    2526    };
  • trunk/MagicSoft/Mars/mimage/MImgCleanStd.cc

    r4702 r4716  
    549549        fData->CalcCleaningLevelDemocratic(*fEvt, *fPed, *fCam);
    550550        break;
    551         /*
    552551    case kProbability:
    553552        fData->CalcCleaningProbability(*fEvt, *fPed, *fCam, fTime);
    554         break;*/
     553        break;
    555554    default:
    556555        break;
  • trunk/MagicSoft/Mars/mjobs/MSequence.h

    r4695 r4716  
    4343    void Print(Option_t *o="") const;
    4444
     45    Bool_t IsValid() const { return fSequence!=(UInt_t)-1; }
     46
    4547    void SetupPedRuns(MDirIter &iter, const char *path=0) const;
    4648    void SetupDatRuns(MDirIter &iter, const char *path=0) const;
  • trunk/MagicSoft/Mars/msignal/MArrivalTime.h

    r4714 r4716  
    4040    const TArrayF &GetDataErr() const { return fDataErr; }
    4141
    42     Double_t operator[](int i) { return fData[i]; }
     42    Double_t operator[](int i)       { return fData[i]; }
     43    Double_t operator[](int i) const { return fData[i]; }
    4344
    4445    // Do not do such things! It is highly dangerous to use two very similar operators
Note: See TracChangeset for help on using the changeset viewer.