Ignore:
Timestamp:
09/09/04 17:23:23 (20 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/mcalib
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h

    r4858 r4900  
    66
    77#pragma link C++ class MCalibColorSet+;
     8#pragma link C++ class MCalibColorSteer+;
    89
    910#pragma link C++ class MCalibrate+;
    … …  
    1112#pragma link C++ class MCalibrateRelTimes+;
    1213
     14#pragma link C++ class MCalibrationIntensityCam+;
     15#pragma link C++ class MCalibrationIntensityChargeCam+;
     16#pragma link C++ class MCalibrationIntensityQECam+;
     17#pragma link C++ class MCalibrationIntensityRelTimeCam+;
    1318#pragma link C++ class MCalibrationCam+;
    1419#pragma link C++ class MCalibrationPix+;
    … …  
    3338
    3439#pragma link C++ class MHCalibrationCam+;
     40#pragma link C++ class MHCalibrationPix+;
    3541#pragma link C++ class MHCalibrationChargeCam+;
    3642#pragma link C++ class MHCalibrationChargePix+;

  • trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc

    r4873 r4900  
    101101using namespace std;
    102102
    103 const Int_t    MHGausEvents::fgBinsAfterStripping   = 40;
    104 const Float_t  MHGausEvents::fgBlackoutLimit        = 5.;
    105103const Int_t    MHGausEvents::fgNDFLimit             = 2;
    106 const Float_t  MHGausEvents::fgPickupLimit          = 5.;
    107104const Float_t  MHGausEvents::fgProbLimit            = 0.001;
    108 const Int_t    MHGausEvents::fgPowerProbabilityBins = 20;
     105const Int_t    MHGausEvents::fgPowerProbabilityBins = 30;
    109106// --------------------------------------------------------------------------
    110107//
    … …  
    114111// - the default Probability Limit for fProbLimit  (fgProbLimit)
    115112// - the default NDF Limit for fNDFLimit           (fgNDFLimit)
    116 // - the default number for fPickupLimit           (fgPickupLimit)
    117 // - the default number for fBlackoutLimit         (fgBlackoutLimit)
    118113// - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping)
    119114// - the default name of the fHGausHist            ("HGausHist")
    … …  
    132127//
    133128MHGausEvents::MHGausEvents(const char *name, const char *title)
    134     : fEventFrequency(0), fHPowerProbability(NULL),
     129    : fHPowerProbability(NULL),
    135130      fPowerSpectrum(NULL),
     131      fFGausFit(NULL), fFExpFit(NULL),
     132      fFirst(0.),
    136133      fGraphEvents(NULL), fGraphPowerSpectrum(NULL),
    137       fFGausFit(NULL), fFExpFit(NULL),
    138       fFirst(0.), fLast(100.),
    139       fNbins(100), fPixId(-1)
     134      fLast(100.), fNbins(100)
    140135{
    141136
    … …  
    146141 
    147142  SetBinsAfterStripping();
    148   SetBlackoutLimit();
    149143  SetNDFLimit();
    150   SetPickupLimit();
    151144  SetPowerProbabilityBins();
    152145  SetProbLimit();
    … …  
    206199// Sets:
    207200// - all other pointers to NULL
    208 // - all variables to 0., except fPixId to -1 and keep fEventFrequency
     201// - all variables to 0.
    209202// - all flags to kFALSE
    210203//
    … …  
    226219  fProb              = 0.;
    227220
    228   fSaturated         = 0;
    229221  fCurrentSize       = 0;
    230222
    … …  
    271263
    272264
    273 // -----------------------------------------------------------------------------
    274 //
    275 // Bypasses the Gauss fit by taking mean and RMS from the histogram
    276 //
    277 // Errors are determined in the following way:
    278 // MeanErr  = RMS / Sqrt(entries)
    279 // SigmaErr = RMS / (2.*Sqrt(entries) )
    280 //
    281 void MHGausEvents::BypassFit()
    282 {
    283 
    284   const Stat_t entries = fHGausHist.GetEntries();
    285  
    286   if (entries <= 0.)
    287     {
    288       *fLog << warn << GetDescriptor()
    289             << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
    290       return;
    291     }
    292  
    293   fMean     = fHGausHist.GetMean();
    294   fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
    295   fSigma    = fHGausHist.GetRMS() ;
    296   fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
    297 }
    298 
    299 
    300 
    301 // --------------------------------------------------------------------------
    302 //
    303 // - Set fPixId to id
    304 //
    305 // Add id to names and titles of:
    306 // - fHGausHist
    307 //
    308 void MHGausEvents::ChangeHistId(const Int_t id)
    309 {
    310 
    311   fPixId = id;
    312 
    313   fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
    314   fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
    315 
    316   fName  = Form("%s%d", fName.Data(),  id);
    317   fTitle = Form("%s%d", fTitle.Data(), id);
    318 
    319 }
    320265
    321266// -----------------------------------------------------------------------------
    … …  
    328273  Float_t *xaxis = new Float_t[n]; 
    329274
    330   if (fEventFrequency)
    331     for (Int_t i=0;i<n;i++)
    332       xaxis[i] = (Float_t)i/fEventFrequency;
    333   else
    334     for (Int_t i=0;i<n;i++)
    335       xaxis[i] = (Float_t)i;
    336 
     275  for (Int_t i=0;i<n;i++)
     276    xaxis[i] = (Float_t)i;
     277 
    337278  return xaxis;
    338279                 
    … …  
    344285// Create the fourier spectrum using the class MFFT.
    345286// The result is projected into a histogram and fitted by an exponential
    346 //   fHPowerProbability->SetDirectory(NULL);
     287//
    347288void MHGausEvents::CreateFourierSpectrum()
    348289{
    … …  
    353294  if (fEvents.GetSize() < 8)
    354295    {
    355       *fLog << warn << "Cannot create Fourier spectrum in pixel: " << fPixId
     296      *fLog << warn << "Cannot create Fourier spectrum in: " << fName
    356297            << ". Number of events smaller than 8 " << endl;
    357298      return;
    358299    }
    359300 
    360 
    361301  //
    362302  // The number of entries HAS to be a potence of 2,
    … …  
    436376  const Int_t n = fEvents.GetSize();
    437377
     378  if (n==0)
     379    return;
     380
    438381  fGraphEvents = new TGraph(n,CreateEventXaxis(n),fEvents.GetArray()); 
    439382  fGraphEvents->SetTitle("Evolution of Events with time");
    440   fGraphEvents->GetXaxis()->SetTitle((fEventFrequency) ? "Time [s]" : "Event Nr.");
    441383  fGraphEvents->GetYaxis()->SetTitle(fHGausHist.GetXaxis()->GetTitle());
    442384  fGraphEvents->GetYaxis()->CenterTitle();
    … …  
    457399  fGraphPowerSpectrum = new TGraph(n,CreatePSDXaxis(n),fPowerSpectrum->GetArray()); 
    458400  fGraphPowerSpectrum->SetTitle("Power Spectrum Density");
    459   fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency");
    460401  fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)");
    461402  fGraphPowerSpectrum->GetYaxis()->CenterTitle();
    … …  
    473414  Float_t *xaxis = new Float_t[n];
    474415
    475   if (fEventFrequency)
    476     for (Int_t i=0;i<n;i++)
    477       xaxis[i] = 0.5*(Float_t)i*fEventFrequency/n;
    478   else
    479     for (Int_t i=0;i<n;i++)
    480       xaxis[i] = 0.5*(Float_t)i/n;
    481 
     416  for (Int_t i=0;i<n;i++)
     417    xaxis[i] = 0.5*(Float_t)i/n;
     418 
    482419  return xaxis;
    483420                 
    … …  
    534471  pad->cd(1);
    535472
    536   if (!IsEmpty())
     473  if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow())
    537474    gPad->SetLogy();
    538475
    … …  
    571508    CreateGraphEvents();
    572509
     510  if (!fGraphEvents)
     511    return;
     512
    573513  fGraphEvents->SetBit(kCanDelete);
    574514  fGraphEvents->SetTitle("Events with time");
    … …  
    666606  //
    667607  // First, strip the zeros from the edges which contain only zeros and rebin
    668   // to about fBinsAfterStripping bins.
     608  // to fBinsAfterStripping bins. If fBinsAfterStripping is 0, reduce bins only by
     609  // a factor 10.
    669610  //
    670611  // (ATTENTION: The Chisquare method is more sensitive,
    671612  // the _less_ bins, you have!)
    672613  //
    673   StripZeros(&fHGausHist,fBinsAfterStripping);
     614  StripZeros(&fHGausHist,
     615             fBinsAfterStripping ? fBinsAfterStripping
     616             : (fNbins > 1000 ? fNbins/10 : 0));
    674617 
    675618  TAxis *axe = fHGausHist.GetXaxis();
    … …  
    693636    {
    694637    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
    695           << "in pixel: " << fPixId << endl;
     638          << "in: " << fName << endl;
    696639    return kFALSE;
    697640    }
    … …  
    731674}
    732675
    733 // -------------------------------------------------------------------------------
    734 //
    735 // Return the number of "blackout" events, which are events with values higher
    736 // than fBlackoutLimit sigmas from the mean
    737 //
    738 //
    739 const Double_t MHGausEvents::GetBlackout() const
    740 {
    741  
    742   if ((fMean == 0.) && (fSigma == 0.))
    743     return -1.;
    744 
    745   const Int_t first = fHGausHist.GetXaxis()->GetFirst();
    746   const Int_t last  = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma);
    747 
    748   if (first >= last)
    749     return 0.;
    750  
    751   return fHGausHist.Integral(first, last, "width");
    752 }
    753676
    754677const Double_t MHGausEvents::GetChiSquare()  const
    … …  
    786709}
    787710
    788 
    789 // -------------------------------------------------------------------------------
    790 //
    791 // Return the number of "pickup" events, which are events with values higher
    792 // than fPickupLimit sigmas from the mean
    793 //
    794 //
    795 const Double_t MHGausEvents::GetPickup() const
    796 {
    797  
    798   if ((fMean == 0.) && (fSigma == 0.))
    799     return -1.;
    800 
    801   const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
    802   const Int_t last  = fHGausHist.GetXaxis()->GetLast();
    803 
    804   if (first >= last)
    805     return 0.;
    806  
    807   return fHGausHist.Integral(first, last, "width");
    808 }
    809711
    810712
    … …  
    890792 
    891793  *fLog << all                                                        << endl;
    892   *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId     << endl;
     794  *fLog << all << "Results of the Gauss Fit in: " << fName            << endl;
    893795  *fLog << all << "Mean: "             << GetMean()                   << endl;
    894796  *fLog << all << "Sigma: "            << GetSigma()                  << endl;
    … …  
    900802}
    901803
    902 // --------------------------------------------------------------------------
    903 //
    904 // Re-normalize the results, has to be overloaded
    905 //
    906 void  MHGausEvents::Renorm()
    907 {
    908 }
    909 
    910 // -----------------------------------------------------------------------------
    911 //
    912 // If flag IsGausFitOK() is set (histogram already successfully fitted),
    913 // returns kTRUE
    914 //
    915 // If both fMean and fSigma are still zero, call FitGaus()
    916 //
    917 // Repeats the Gauss fit in a smaller range, defined by:
    918 //
    919 // min = GetMean() - fBlackoutLimit * GetSigma();
    920 // max = GetMean() + fPickupLimit   * GetSigma();
    921 //
    922 // The fit results are retrieved and stored in class-own variables. 
    923 //
    924 // A flag IsGausFitOK() is set according to whether the fit probability
    925 // is smaller or bigger than fProbLimit, whether the NDF is bigger than
    926 // fNDFLimit and whether results are NaNs.
    927 //
    928 Bool_t MHGausEvents::RepeatFit(const Option_t *option)
    929 {
    930 
    931   if (IsGausFitOK())
    932     return kTRUE;
    933 
    934   if ((fMean == 0.) && (fSigma == 0.))
    935     return FitGaus();
    936 
    937   //
    938   // Get new fitting ranges
    939   //
    940   Axis_t rmin = fMean - fBlackoutLimit * fSigma;
    941   Axis_t rmax = fMean + fPickupLimit   * fSigma;
    942 
    943   Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
    944   Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
    945 
    946   fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
    947 
    948   fHGausHist.Fit(fFGausFit,option);
    949 
    950   fMean     = fFGausFit->GetParameter(1);
    951   fSigma    = fFGausFit->GetParameter(2);
    952   fMeanErr  = fFGausFit->GetParError(1) ;
    953   fSigmaErr = fFGausFit->GetParError(2) ;
    954   fProb     = fFGausFit->GetProb()      ;     
    955 
    956   //
    957   // The fit result is accepted under condition:
    958   // 1) The results are not nan's
    959   // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
    960   // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
    961   //
    962   if (   TMath::IsNaN ( fMean     )
    963       || TMath::IsNaN ( fMeanErr  )
    964       || TMath::IsNaN ( fProb     )   
    965       || TMath::IsNaN ( fSigma    )
    966       || TMath::IsNaN ( fSigmaErr )
    967       || fFGausFit->GetNDF() < fNDFLimit
    968       || fProb < fProbLimit )
    969     return kFALSE;
    970  
    971   SetGausFitOK(kTRUE);
    972   return kTRUE;
    973 
    974 }
    975804
    976805// --------------------------------------------------------------------------

  • trunk/MagicSoft/Mars/mcalib/MHGausEvents.h

    r4882 r4900  
    2121private:
    2222
    23   const static Int_t    fgBinsAfterStripping;   //! Default for fBinsAfterStripping   (now set to: 40)
    24   const static Float_t  fgBlackoutLimit;        //! Default for fBlackoutLimit        (now set to: 5. )
    2523  const static Int_t    fgNDFLimit;             //! Default for fNDFLimit             (now set to: 2)
    2624  const static Float_t  fgProbLimit;            //! Default for fProbLimit            (now set to: 0.001)
    27   const static Float_t  fgPickupLimit;          //! Default for fPickupLimit          (now set to: 5. )
    2825  const static Int_t    fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20)
    2926 
    3027  Int_t    fBinsAfterStripping;        // Bins for the Gauss Histogram after stripping off the zeros at both ends
    3128  Int_t    fCurrentSize;               // Current size of the array fEvents
    32   Float_t  fEventFrequency;            // Event frequency in Hertz (to be set)
    3329  Byte_t   fFlags;                     // Bit field for the fit result bits
    3430  Int_t    fPowerProbabilityBins;      // Bins for the projected power spectrum
    … …  
    3632  TH1I    *fHPowerProbability;         // Fourier transform of fEvents projected on y-axis
    3733  TArrayF *fPowerSpectrum;             // Fourier transform of fEvents
    38   TGraph  *fGraphEvents;               //! TGraph to display the event array (will not be cloned!!)
    39   TGraph  *fGraphPowerSpectrum;        //! TGraph to display the power spectrum array (will not be cloned!!)
    4034
    41   enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK, kExcluded }; // Bits for information about fit results
     35  enum  { kGausFitOK,
     36          kExpFitOK,
     37          kFourierSpectrumOK,
     38          kExcluded };                 // Bits for information about fit results
    4239 
    4340protected:
    4441
    45   Float_t  fBlackoutLimit;             // Lower number sigmas from mean until event is considered blackout
    4642  TArrayF  fEvents;                    // Array which holds the entries of GausHist
    4743  TF1     *fFGausFit;                  // Gauss fit for fHGausHist
    4844  TF1     *fFExpFit;                   // Exponential fit for FHPowerProbability
    4945  Axis_t   fFirst;                     // Lower histogram edge  for fHGausHist (used by InitBins())
     46  TGraph  *fGraphEvents;               //! TGraph to display the event array (will not be cloned!!)
     47  TGraph  *fGraphPowerSpectrum;        //! TGraph to display the power spectrum array (will not be cloned!!)
    5048  TH1F     fHGausHist;                 // Histogram to hold the Gaussian distribution
    5149  Axis_t   fLast;                      // Upper histogram edge  for fHGausHist (used by InitBins())
    … …  
    5452  Int_t    fNbins;                     // Number histogram bins for fHGausHist (used by InitBins())
    5553  Int_t    fNDFLimit;                  // NDF limit for judgement if fit is OK
    56   Int_t    fSaturated;                 // Number of events classified as saturated
    5754  Double_t fSigma;                     // Sigma of the Gauss fit
    5855  Double_t fSigmaErr;                  // Error of the sigma of the Gauss fit
    59   Float_t  fPickupLimit;               // Upper number sigmas from mean until event is considered pickup
    60   Int_t    fPixId;                     // Pixel ID
    6156  Double_t fProb;                      // Probability of the Gauss fit
    6257  Float_t  fProbLimit;                 // Probability limit for judgement if fit is OK
    6358
    64   Float_t *CreateEventXaxis(Int_t n);  // Create an x-axis for the Event TGraphs
    65   Float_t *CreatePSDXaxis(Int_t n);    // Create an x-axis for the PSD TGraphs
     59  virtual Float_t *CreateEventXaxis(Int_t n);  // Create an x-axis for the Event TGraphs
     60  virtual Float_t *CreatePSDXaxis(Int_t n);    // Create an x-axis for the PSD TGraphs
     61  virtual void     CreateGraphEvents();        // Create the TGraph fGraphEvents of fEvents
     62  virtual void     CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum
    6663
    6764  void DrawPowerSpectrum(TVirtualPad &pad, Int_t i);  // Draw graph of fPowerSpectrum and fHPowerProbability
    6865
    6966  // Setters
    70   void  SetBinsAfterStripping   ( const Int_t nbins=fgBinsAfterStripping   ) { fBinsAfterStripping  =nbins; }
     67  void  SetBinsAfterStripping   ( const Int_t nbins=0   )                    { fBinsAfterStripping  =nbins; }
    7168  void  SetPowerProbabilityBins ( const Int_t nbins=fgPowerProbabilityBins ) { fPowerProbabilityBins=nbins; }
    7269
    73  public:
     70public:
    7471
    7572  MHGausEvents(const char* name=NULL, const char* title=NULL);
    … …  
    8582 
    8683  // Getters
    87   const Double_t GetBlackout()           const; 
    8884  const Double_t GetChiSquare()          const;
    8985  const Double_t GetExpChiSquare()       const;
    … …  
    110106  const Int_t    GetNdf()                const;
    111107  const Double_t GetOffset()             const;
    112   const Double_t GetPickup()             const;
    113   const Int_t    GetPixId()              const { return fPixId;              }
    114108        TArrayF *GetPowerSpectrum()            { return fPowerSpectrum;      } 
    115109  const TArrayF *GetPowerSpectrum()      const { return fPowerSpectrum;      }
    116110  const Double_t GetProb()               const { return fProb;               }
    117   const Float_t  GetSaturated()          const { return fSaturated;          }
    118111  const Double_t GetSigma()              const { return fSigma;              }
    119112  const Double_t GetSigmaErr()           const { return fSigmaErr;           }
    … …  
    137130                   const Double_t xmin=0.,
    138131                   const Double_t xmax=0.);       // Fit the histogram HGausHist with a Gaussian
    139   Bool_t RepeatFit(const Option_t *option="RQ0"); // Repeat fit within limits defined by fPickupLimit
    140   void   BypassFit();                             // Take mean and RMS from the histogram
    141132 
    142133  // Prints
    … …  
    144135 
    145136  // Setters
    146   void  SetBlackoutLimit    ( const Float_t  lim=fgBlackoutLimit ) { fBlackoutLimit  = lim; }
    147   void  SetEventFrequency   ( const Float_t  f                   ) { fEventFrequency = f;   }
    148137  void  SetExcluded         ( const Bool_t   b=kTRUE             ); 
    149138  void  SetExpFitOK         ( const Bool_t   b=kTRUE             );
    … …  
    156145  void  SetNbins            ( const Int_t    i                   ) { fNbins          = i;   } 
    157146  void  SetNDFLimit         ( const Int_t    lim=fgNDFLimit      ) { fNDFLimit       = lim; } 
    158   void  SetPickupLimit      ( const Float_t  lim=fgPickupLimit   ) { fPickupLimit    = lim; }
    159   void  SetPixId            ( const Int_t    i                   ) { fPixId          = i;   }
    160147  void  SetProb             ( const Double_t d                   ) { fProb           = d;   }
    161148  void  SetProbLimit        ( const Float_t  lim=fgProbLimit     ) { fProbLimit      = lim; }
    162   void  SetSaturated        ( const Int_t    i                   ) { fSaturated     += i;   }
    163149  void  SetSigma            ( const Double_t d                   ) { fSigma          = d;   }
    164150  void  SetSigmaErr         ( const Double_t d                   ) { fSigmaErr       = d;   }
    165151
    166   // Miscelleaneous
    167   virtual void ChangeHistId(const Int_t id);      // Changes names and titles of the histogram
    168   virtual void Renorm();                          // Re-normalize the results
    169  
    170152  void CreateFourierSpectrum();                   // Create the fourier spectrum out of fEvents
    171   void CreateGraphEvents();                       // Create the TGraph fGraphEvents of fEvents
    172   void CreateGraphPowerSpectrum();                // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum
    173153 
    174154  ClassDef(MHGausEvents, 1) // Base class for events with Gaussian distributed values

  • trunk/MagicSoft/Mars/mcalib/Makefile

    r4858 r4900  
    3333SRCFILES = MCalibrate.cc \
    3434           MCalibColorSet.cc \
     35           MCalibColorSteer.cc \
    3536           MCalibrateData.cc \
    3637           MCalibrateRelTimes.cc \
     38           MCalibrationIntensityCam.cc \
     39           MCalibrationIntensityChargeCam.cc \
     40           MCalibrationIntensityQECam.cc \
     41           MCalibrationIntensityRelTimeCam.cc \
    3742           MCalibrationCam.cc \
    3843           MCalibrationPix.cc  \
    … …  
    5964           MHCalibrationChargePix.cc \
    6065           MHCalibrationCam.cc \
     66           MHCalibrationPix.cc \
    6167           MHCalibrationChargeCam.cc \
    6268           MHCalibrationChargeHiGainPix.cc \
Note: See TracChangeset for help on using the changeset viewer.