Changeset 9226 for trunk


Ignore:
Timestamp:
01/17/09 14:52:43 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
9 edited

Legend:

Unmodified
Added
Removed

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

    r9219 r9226  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.43 2009-01-14 12:31:36 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.44 2009-01-17 14:52:37 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    2020!   Author(s): Thomas Bretz  3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
    2121!
    22 !   Copyright: MAGIC Software Development, 2000-2005
     22!   Copyright: MAGIC Software Development, 2000-2009
    2323!
    2424!
    … …  
    5050#ifndef ROOT_TComplex
    5151#include <TComplex.h>
     52#endif
     53
     54#ifndef ROOT_TRandom
     55#include <TRandom.h>  // gRandom in RndmExp
    5256#endif
    5357
    … …  
    875879    e = error.Atof();
    876880}
     881
     882Double_t MMath::RndmExp(Double_t tau)
     883{
     884    // returns an exponential deviate.
     885    //
     886    //          exp( -t/tau )
     887
     888    const Double_t x = gRandom->Rndm(); // uniform on ] 0, 1 ]
     889
     890    return -tau * TMath::Log(x);       // convert to exponential distribution
     891}

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

    r9219 r9226  
    8787
    8888    void Format(Double_t &v, Double_t &e);
     89
     90    Double_t RndmExp(Double_t tau);
    8991}
    9092

  • trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h

    r7942 r9226  
    11#ifndef __CINT__
    22
     3#include "../mbase/MArrayF.h"
     4
    35#endif // __CINT__

  • trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc

    r9212 r9226  
    1818!   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
    1919!
    20 !   Copyright: MAGIC Software Development, 2002-2008
     20!   Copyright: MAGIC Software Development, 2002-2009
    2121!
    2222!
    … …  
    120120// --------------------------------------------------------------------------
    121121//
    122 // Returns the highest x value in [min;max[ at which the spline in
    123 // the bin i is equal to y
    124 //
    125 // min and max are defined to be [0;1]
    126 //
    127 // The default for min is 0, the default for max is 1
    128 // The defaule for y is 0
    129 //
    130 Double_t MExtralgoSpline::FindY(Int_t i, Bool_t downwards, Double_t y, Double_t min, Double_t max) const
     122// Solve the polynomial
     123//
     124//    y = a*x^3 + b*x^2 + c*x + d'
     125//    0 = a*x^3 + b*x^2 + c*x + d' - y
     126//
     127// to find y in the i-th bin. Return the result as x1, x2, x3 and the return
     128// code from MMath::SolvPol3.
     129//
     130Int_t MExtralgoSpline::SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const
    131131{
    132132    // y = a*x^3 + b*x^2 + c*x + d'
    … …  
    145145    //     return -2;
    146146
     147    return MMath::SolvePol3(a, b, c, d, x1, x2, x3);
     148}
     149
     150// --------------------------------------------------------------------------
     151//
     152// Returns the highest x value in [min;max[ at which the spline in
     153// the bin i is equal to y
     154//
     155// min and max must be in the range [0;1]
     156//
     157// The default for min is 0, the default for max is 1
     158// The default for y is 0
     159//
     160Double_t MExtralgoSpline::FindYdn(Int_t i, Double_t y, Double_t min, Double_t max) const
     161{
    147162    Double_t x1, x2, x3;
    148     const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3);
    149 
    150     if (downwards==kTRUE)
    151     {
    152         Double_t x = -1;
    153 
    154         if (rc>0 && x1>=min && x1<max && x1>x)
    155             x = x1;
    156         if (rc>1 && x2>=min && x2<max && x2>x)
    157             x = x2;
    158         if (rc>2 && x3>=min && x3<max && x3>x)
    159             x = x3;
    160 
    161         return x<0 ? -2 : x+i;
    162     }
    163     else
    164     {
    165         Double_t x = 2;
    166 
    167         if (rc>0 && x1>min && x1<=max && x1<x)
    168             x = x1;
    169         if (rc>1 && x2>min && x2<=max && x2<x)
    170             x = x2;
    171         if (rc>2 && x3>min && x3<=max && x3<x)
    172             x = x3;
    173 
    174         return x>1 ? -2 : x+i;
    175     }
    176 
    177     return -2;
     163    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
     164
     165    Double_t x = -1;
     166
     167    if (rc>0 && x1>=min && x1<max && x1>x)
     168        x = x1;
     169    if (rc>1 && x2>=min && x2<max && x2>x)
     170        x = x2;
     171    if (rc>2 && x3>=min && x3<max && x3>x)
     172        x = x3;
     173
     174    return x<0 ? -2 : x+i;
     175}
     176
     177// --------------------------------------------------------------------------
     178//
     179// Returns the lowest x value in [min;max[ at which the spline in
     180// the bin i is equal to y
     181//
     182// min and max must be in the range [0;1]
     183//
     184// The default for min is 0, the default for max is 1
     185// The default for y is 0
     186//
     187Double_t MExtralgoSpline::FindYup(Int_t i, Double_t y, Double_t min, Double_t max) const
     188{
     189    Double_t x1, x2, x3;
     190    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
     191
     192    Double_t x = 2;
     193
     194    if (rc>0 && x1>min && x1<=max && x1<x)
     195        x = x1;
     196    if (rc>1 && x2>min && x2<=max && x2<x)
     197        x = x2;
     198    if (rc>2 && x3>min && x3<=max && x3<x)
     199        x = x3;
     200
     201    return x>1 ? -2 : x+i;
    178202}
    179203
    … …  
    181205//
    182206// Search analytically downward for the value y of the spline, starting
    183 // at x, until x==0. If y is not found -2 is returned.
     207// at x, until x==0. If y is not found or out of range -2 is returned.
    184208//
    185209Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
    … …  
    192216        return -2;
    193217
    194     Double_t rc = FindY(i, kTRUE, y, 0, x-i);
     218    Double_t rc = FindYdn(i, y, 0, x-i);
    195219    while (--i>=0 && rc<0)
    196         rc = FindY(i, kTRUE, y);
     220        rc = FindYdn(i, y);
    197221
    198222    return rc;
    199223}
    200224
     225// --------------------------------------------------------------------------
     226//
     227// Search analytically upwards for the value y of the spline, starting
     228// at x, until x==fNum-1. If y is not found or out of range -2 is returned.
     229//
    201230Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const
    202231{
    … …  
    208237        return -2;
    209238
    210     Double_t rc = FindY(i, kFALSE, y, x-i, 1.);
     239    Double_t rc = FindYup(i, y, x-i, 1.);
    211240    while (++i<fNum-1 && rc<0)
    212         rc = FindY(i, kFALSE, y);
     241        rc = FindYup(i, y);
    213242
    214243    return rc;
    … …  
    222251Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
    223252{
    224     // In the future we will calculate the intgeral analytically.
    225     // It has been tested that it gives identical results within
    226     // acceptable differences.
    227 
    228253    // We allow extrapolation of 1/2 slice.
    229254    const Float_t min = fRiseTime;        //-0.5+fRiseTime;
    … …  
    236261
    237262    return EvalInteg(pos-fRiseTime, pos+fFallTime);
     263}
     264
     265MArrayF MExtralgoSpline::GetIntegral(bool norm) const
     266{
     267    MArrayF val(fNum);
     268
     269    //val[0] = 0;
     270
     271    Double_t integ = 0;
     272    for (int i=0; i<fNum-1; i++)
     273    {
     274        integ += EvalInteg(i);
     275
     276        val[i+1] = integ;
     277    }
     278
     279    if (norm)
     280        for (int i=0; i<fNum-1; i++)
     281            val[i+1] /= val[fNum-1];
     282
     283    return val;
    238284}
    239285

  • trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h

    r9212 r9226  
    66#endif
    77
     8class MArrayF;
    89class TComplex;
    910
    … …  
    7576    }
    7677
    77     Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
     78    Int_t SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const;
     79    Double_t FindYdn(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
     80    Double_t FindYup(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
     81    //Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
    7882
    7983    Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
    … …  
    153157    }
    154158
    155     // Identical to sum EvalInteg(i, 0, 1) for i=0 to i<b but much faster
     159    // Identical to sum of EvalInteg(i, 0, 1) for i=a to i=b-1,
     160    // but much faster
     161    // It is identical to EvalInteg(fVal[a], fVal[b])
    156162    // Be carefull: NO RANGECHECK!
    157163    inline Double_t EvalInteg(Int_t a, Int_t b) const
    … …  
    329335    Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); }
    330336    Double_t SearchYup(Float_t y) const { return SearchYup(0,    y); }
     337
     338    MArrayF GetIntegral(bool norm=true) const;
    331339};
    332340

  • trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc

    r8565 r9226  
    11/* ======================================================================== *\
    2 ! $Name: not supported by cvs2svn $:$Id: MPedestalSubtractedEvt.cc,v 1.6 2007-06-16 22:08:00 tbretz Exp $
     2! $Name: not supported by cvs2svn $:$Id: MPedestalSubtractedEvt.cc,v 1.7 2009-01-17 14:52:41 tbretz Exp $
    33! --------------------------------------------------------------------------
    44!
    … …  
    3232//
    3333/////////////////////////////////////////////////////////////////////////////
    34 
    3534#include "MPedestalSubtractedEvt.h"
     35
     36#include "MLogManip.h"
    3637
    3738ClassImp(MPedestalSubtractedEvt);
    … …  
    185186}
    186187
     188void MPedestalSubtractedEvt::Print(Option_t *o) const
     189{
     190    *fLog << all << GetDescriptor() << endl;
     191    *fLog << " Num Pixels:  " << fNumPixels << " (" << fNumSamples << " samples)" << endl;
     192    *fLog << " Samples raw:" << hex << endl;;
     193    for (UInt_t idx=0; idx<fNumPixels; idx++)
     194    {
     195        *fLog << setw(4) << dec << idx << hex << ":";
     196        for (UInt_t i=0; i<fNumSamples; i++)
     197            *fLog << " " << fSamplesRaw[idx*fNumSamples+i];
     198        *fLog << endl;
     199    }
     200    *fLog << dec << endl;
     201    *fLog << " Samples:" << endl;;
     202    for (UInt_t idx=0; idx<fNumPixels; idx++)
     203    {
     204        *fLog << setw(4) << idx << ":";
     205        for (UInt_t i=0; i<fNumSamples; i++)
     206            *fLog << " " << fSamples[idx*fNumSamples+i];
     207        *fLog << endl;
     208    }
     209    *fLog << endl;
     210}
     211
    187212/*
    188213#include <TSpline.h>

  • trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h

    r8571 r9226  
    2020{
    2121private:
    22     MArrayF fSamples;        // list of all samples with pedestal subtracted
    23     MArrayS fSamplesRaw;     // list of all samples (raw)
     22    MArrayF fSamples;         // list of all samples with pedestal subtracted
     23    MArrayS fSamplesRaw;      // list of all samples (raw)
    2424
    25     UInt_t fNumSamples;      // number of samples per pixel
    26     UInt_t fNumPixels;       // number of pixels
     25    UInt_t fNumSamples;       // number of samples per pixel
     26    UInt_t fNumPixels;        // number of pixels
    2727
    2828public:
    … …  
    114114    }
    115115
     116    void Print(Option_t *o="") const;
     117
    116118    Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const { return kTRUE; }
    117119    void   DrawPixelContent(Int_t num) const { }

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

    r9212 r9226  
    284284    name += pix.GetPixelId();
    285285
     286
    286287    Bool_t same = str.Contains("same");
    287288
    … …  
    392393}
    393394
     395void MRawEvtData::ResetPixels()
     396{
     397    if (!fRunHeader)
     398        return;
     399
     400    // FIXME: Better give NumNormalPixels if numconnected==0 ?
     401
     402    ResetPixels(fRunHeader->GetNumNormalPixels(), fRunHeader->GetNumNormalPixels()-1);
     403}
     404
    394405// --------------------------------------------------------------------------
    395406//
    … …  
    435446void MRawEvtData::Set(const MArrayI &data)
    436447{
    437     const UInt_t n = fRunHeader->GetNumSamples()*fConnectedPixels;
     448    fConnectedPixels = fHiGainPixId->GetSize();
     449
     450    UInt_t n = fConnectedPixels*fRunHeader->GetNumSamplesHiGain();
     451    if (n!=data.GetSize())
     452    {
     453        *fLog << err << "MRawEvtData::Set: Size mismatch." << endl;
     454        *fLog << " fConnectedPixels="   << fConnectedPixels << endl;
     455        *fLog << " NumHiGainSamples="   << fRunHeader->GetNumSamplesHiGain() << endl;
     456        *fLog << " data.GetSize()="     << data.GetSize() << endl;
     457        return;
     458    }
    438459
    439460    Byte_t *dest = fHiGainFadcSamples->GetArray();
    … …  
    448469            Byte_t *ptr = reinterpret_cast<Byte_t*>(dest);
    449470            while (src<end)
    450                 *ptr++ = *src++;
     471                *ptr++ = Byte_t(*src++);
    451472        }
    452473        return;
    … …  
    456477            UShort_t *ptr = reinterpret_cast<UShort_t*>(dest);
    457478            while (src<end)
    458                 *ptr++ = *src++;
     479                *ptr++ = UShort_t(*src++);
    459480        }
    460481        return;
    … …  
    464485        return;
    465486    }
     487}
     488
     489void MRawEvtData::SetIndices(const MArrayS &idx)
     490{
     491    if (idx.GetSize()!=fHiGainPixId->GetSize())
     492        return;
     493
     494    memcpy(fHiGainPixId->GetArray(), idx.GetArray(), idx.GetSize()*sizeof(UShort_t));
     495}
     496
     497void MRawEvtData::SetIndices()
     498{
     499    for (UInt_t i=0; i<fHiGainPixId->GetSize(); i++)
     500        (*fHiGainPixId)[i] = i;
    466501}
    467502

  • trunk/MagicSoft/Mars/mraw/MRawEvtData.h

    r9212 r9226  
    5757    }
    5858
     59
    5960public:
    6061    MRawEvtData(const char *name=NULL, const char *title=NULL);
    … …  
    7475    void Draw (Option_t * = NULL);
    7576
     77    void ResetPixels();
    7678    void ResetPixels(UShort_t npix, UShort_t maxid);
    7779    void AddPixel(UShort_t nOfPixel, const TArrayC &data);
    7880    void Set(const MArrayI &data);
     81    void SetIndices(const MArrayS &idx);
     82    void SetIndices();
    7983
    8084    UShort_t GetNumHiGainSamples() const;
Note: See TracChangeset for help on using the changeset viewer.