Changeset 3636


Ignore:
Timestamp:
04/03/04 17:27:50 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
17 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r3635 r3636  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2004/04/02: Markus Gaug
     22   * mcalib/MHGausEvents.[h,cc]
     23   * mcalib/MHCalibrationChargePix.[h,cc]
     24   * mcalib/MHCalibrationChargeHiGainPix.[h,cc]
     25   * mcalib/MHCalibrationChargeLoGainPix.[h,cc]
     26   * mcalib/MHCalibrationRelTimePix.[h,cc]
     27     - put fPixId, fPickup, fPickupLimt, CountPickup(), RepeatFit() and
     28       ChangeHistId() into MHGausEvents (before in the derived classes)
     29     - put fChargeNbins, fChargeFirst, fChargeLast,
     30     - put fRelTimeNbins, fRelTimeFirst, fRelTimeLast together
     31       into MHGausEvents as fNbins, fFirst and fLast
     32
     33   * mcalib/MHCalibrationRelTimePix.[h,cc]
     34     - remove Renormalization to time slices. Need to think about
     35       more direct way to implement
     36
     37   * mcalib/MHCalibrationCam.[h,cc]
     38   * mcalib/MHCalibrationChargeCam.[h,cc]
     39   * mcalib/MHCalibrationRelTimeCam.[h,cc]
     40     - put most of the functionality into the base class MHCalibrationCam
     41     - derived classes overload the functions SetupHists, ReInitHists,
     42       FillHists, FinalizeHists and FinalizeBadPixels.
     43     - functions FitHiGainArrays, FitLoGainArrays, FitHiGainHists,
     44       FitLoGainHists and InitHists can be used from base class.
    2045
    2146
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc

    r3631 r3636  
    5353#include "MLogManip.h"
    5454
     55#include "MCalibrationPix.h"
     56#include "MCalibrationCam.h"
     57
    5558#include "MHGausEvents.h"
    5659
     60#include "MBadPixelsPix.h"
     61#include "MBadPixelsCam.h"
     62
    5763#include "MGeomCam.h"
     64#include "MGeomPix.h"
     65
     66#include "MParList.h"
    5867
    5968ClassImp(MHCalibrationCam);
     
    6170using namespace std;
    6271
     72const Int_t   MHCalibrationCam::fgAverageNbins    = 2000;
    6373const Int_t   MHCalibrationCam::fgPulserFrequency = 500;
    6474// --------------------------------------------------------------------------
    6575//
    6676// Default Constructor.
     77//
     78// Sets:
     79// - all pointers to NULL
    6780//
    6881// Initializes and sets owner of:
     
    7588//
    7689MHCalibrationCam::MHCalibrationCam(const char *name, const char *title)
     90    :  fGeom(NULL), fBadPixels(NULL), fCam(NULL)
    7791{
    7892    fName  = name  ? name  : "MHCalibrationCam";
     
    97111    fAverageLoGainSectors->SetOwner();
    98112
     113    SetAverageNbins();
    99114    SetPulserFrequency();
    100115}
     
    306321// --------------------------------------------------------------------------
    307322//
     323// Gets the pointers to:
     324// - MGeomCam
     325//
     326// Calls SetupHists(const MParList *pList)
     327//
     328// Calls Delete-Function of:
     329// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
     330// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     331// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     332//
     333Bool_t MHCalibrationCam::SetupFill(const MParList *pList)
     334{
     335 
     336  fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
     337  if (!fGeom)
     338  {
     339      *fLog << err << "MGeomCam not found... aborting." << endl;
     340      return kFALSE;
     341  }
     342
     343  fHiGainArray->Delete();
     344  fLoGainArray->Delete();
     345
     346  fAverageHiGainAreas->Delete();
     347  fAverageLoGainAreas->Delete();
     348
     349  fAverageHiGainSectors->Delete();
     350  fAverageLoGainSectors->Delete();
     351
     352  return SetupHists(pList);
     353}
     354
     355
     356Bool_t MHCalibrationCam::SetupHists(const MParList *pList)
     357{
     358  return kTRUE;
     359}
     360
     361// --------------------------------------------------------------------------
     362//
     363// Gets or creates the pointers to:
     364// - MBadPixelsCam
     365//
     366// Searches pointer to:
     367// - MArrivalTimeCam
     368//
     369// Initializes, if empty to MArrivalTimeCam::GetSize() for:
     370// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
     371//
     372// Initializes, if empty to MGeomCam::GetNumAreas() for:
     373// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     374//
     375// Initializes, if empty to MGeomCam::GetNumSectors() for:
     376// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     377//
     378// Initializes TArray's to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively
     379// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()):
     380// - MHCalibrationCam::fAverageAreaNum[area index]
     381// - MHCalibrationCam::fAverageSectorNum[area index]
     382//
     383// Calls InitializeHists() for every entry in:
     384// - MHCalibrationCam::fHiGainArray
     385// - MHCalibrationCam::fAverageHiGainAreas
     386// - MHCalibrationCam::fAverageHiGainSectors
     387//
     388// Sets Titles and Names for the Histograms
     389// - MHCalibrationCam::fAverageHiGainAreas
     390// - MHCalibrationCam::fAverageHiGainSectors
     391//
     392Bool_t MHCalibrationCam::ReInit(MParList *pList)
     393{
     394
     395  fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
     396  if (!fBadPixels)
     397      return kFALSE;
     398
     399
     400  const Int_t npixels  = fGeom->GetNumPixels();
     401  const Int_t nsectors = fGeom->GetNumSectors();
     402  const Int_t nareas   = fGeom->GetNumAreas();
     403
     404  //
     405  // The function TArrayF::Set() already sets all entries to 0.
     406  //
     407  fAverageAreaNum.        Set(nareas);
     408  fAverageAreaSat.        Set(nareas);           
     409  fAverageAreaSigma.      Set(nareas);     
     410  fAverageAreaSigmaErr.   Set(nareas);   
     411  fAverageAreaRelSigma.   Set(nareas);   
     412  fAverageAreaRelSigmaErr.Set(nareas);
     413  fAverageSectorNum.      Set(nsectors);
     414
     415  for (Int_t i=0; i<npixels; i++)
     416    {
     417
     418      if ((*fBadPixels)[i].IsBad())
     419        continue;
     420
     421      fAverageAreaNum  [(*fGeom)[i].GetAidx()  ]++;
     422      fAverageSectorNum[(*fGeom)[i].GetSector()]++;
     423    }
     424
     425  return ReInitHists(pList);
     426}
     427
     428
     429Bool_t MHCalibrationCam::ReInitHists(MParList *pList)
     430{
     431  return kTRUE;
     432}
     433
     434
     435
     436//--------------------------------------------------------------------------------
     437//
     438// Retrieves from MGeomCam:
     439// - number of pixels
     440// - number of pixel areas
     441// - number of sectors
     442//
     443// For all TObjArray's (including the averaged ones), the following steps are performed:
     444//
     445// 1) Test size and return kFALSE if not matching
     446// 2)
     447//
     448Bool_t MHCalibrationCam::Fill(const MParContainer *par, const Stat_t w)
     449{
     450
     451  const Int_t npixels  = fGeom->GetNumPixels();
     452  const Int_t nareas   = fGeom->GetNumAreas();
     453  const Int_t nsectors = fGeom->GetNumSectors();
     454 
     455  if (fHiGainArray->GetEntries() != npixels)
     456    {
     457      gLog << err << "ERROR - Size mismatch... abort." << endl;
     458      return kFALSE;
     459    }
     460 
     461  if (fLoGainArray->GetEntries() != npixels)
     462    {
     463      gLog << err << "ERROR - Size mismatch... abort." << endl;
     464      return kFALSE;
     465    }
     466 
     467  if (fAverageHiGainAreas->GetEntries() != nareas)
     468    {
     469      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
     470      return kFALSE;
     471    }
     472
     473  if (fAverageLoGainAreas->GetEntries() != nareas)
     474    {
     475      *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
     476      return kFALSE;
     477    }
     478
     479  if (fAverageHiGainSectors->GetEntries() != nsectors)
     480    {
     481      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
     482      return kFALSE;
     483    }
     484
     485  if (fAverageLoGainSectors->GetEntries() != nsectors)
     486    {
     487      *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
     488      return kFALSE;
     489    }
     490
     491  return FillHists(par, w);
     492}
     493
     494Bool_t MHCalibrationCam::FillHists(const MParContainer *par, const Stat_t w)
     495{
     496  return kTRUE;
     497}
     498
     499// --------------------------------------------------------------------------
     500//
     501// 1) FinalizeHists()
     502// 2) FinalizeBadPixels()
     503// 3) CalcAverageSigma()
     504//
     505Bool_t MHCalibrationCam::Finalize()
     506{
     507  if (!FinalizeHists())
     508    return kFALSE;
     509
     510  FinalizeBadPixels();
     511  CalcAverageSigma();
     512
     513  return kTRUE;
     514}
     515
     516Bool_t MHCalibrationCam::FinalizeHists()
     517{
     518  return kTRUE;
     519}
     520
     521void MHCalibrationCam::FinalizeBadPixels()
     522{
     523}
     524
     525
     526// -------------------------------------------------------------
     527//
     528// If MBadPixelsPix::IsBad():
     529// - calls MHGausEvents::SetExcluded()
     530//
     531// Calls:
     532// - MHGausEvents::InitBins()
     533// - MHGausEvents::ChangeHistId(i)
     534// - MHGausEvents::SetEventFrequency(fPulserFrequency)
     535//
     536void MHCalibrationCam::InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i)
     537{
     538
     539  if (bad.IsBad())
     540    hist.SetExcluded();
     541
     542  hist.InitBins();
     543  hist.ChangeHistId(i);
     544  hist.SetEventFrequency(fPulserFrequency);
     545         
     546}
     547
     548void MHCalibrationCam::FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
     549                                       MBadPixelsPix::UncalibratedType_t fittyp,
     550                                       MBadPixelsPix::UncalibratedType_t osctyp)
     551{
     552 
     553  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
     554    {
     555     
     556      MHGausEvents &hist = (*this)[i];
     557     
     558      if (hist.IsExcluded())
     559        continue;
     560     
     561      MCalibrationPix &pix    = calcam[i];
     562      MBadPixelsPix   &bad    = badcam[i];
     563     
     564      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
     565     
     566    }
     567
     568  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
     569    {
     570     
     571      MHGausEvents     &hist = GetAverageHiGainArea(j);     
     572      MCalibrationPix  &pix  = calcam.GetAverageArea(j);
     573      MBadPixelsPix    &bad  = calcam.GetAverageBadArea(j);       
     574     
     575      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
     576  }
     577 
     578
     579  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
     580    {
     581     
     582      MHGausEvents     &hist = GetAverageHiGainSector(j);     
     583      MCalibrationPix  &pix  = calcam.GetAverageSector(j);
     584      MBadPixelsPix    &bad  = calcam.GetAverageBadSector(j);       
     585     
     586      FitHiGainHists(hist,pix,bad,fittyp,osctyp);
     587    }
     588
     589}
     590
     591void MHCalibrationCam::FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
     592                                            MBadPixelsPix::UncalibratedType_t fittyp,
     593                                            MBadPixelsPix::UncalibratedType_t osctyp)
     594{
     595 
     596  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
     597    {
     598     
     599      MHGausEvents &hist = (*this)(i);
     600     
     601      if (hist.IsExcluded())
     602        continue;
     603     
     604      MCalibrationPix &pix    = calcam[i];
     605      MBadPixelsPix   &bad    = badcam[i];
     606     
     607      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
     608     
     609    }
     610
     611  for (Int_t j=0; j<fAverageLoGainAreas->GetSize(); j++)
     612    {
     613     
     614      MHGausEvents     &hist = GetAverageLoGainArea(j);     
     615      MCalibrationPix  &pix  = calcam.GetAverageArea(j);
     616      MBadPixelsPix    &bad  = calcam.GetAverageBadArea(j);       
     617     
     618      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
     619  }
     620 
     621
     622  for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
     623    {
     624     
     625      MHGausEvents     &hist = GetAverageLoGainSector(j);     
     626      MCalibrationPix  &pix  = calcam.GetAverageSector(j);
     627      MBadPixelsPix    &bad  = calcam.GetAverageBadSector(j);       
     628     
     629      FitLoGainHists(hist,pix,bad,fittyp,osctyp);
     630    }
     631}
     632
     633//------------------------------------------------------------
     634//
     635// For all averaged areas, the fitted sigma is multiplied with the square root of
     636// the number involved pixels
     637//
     638void MHCalibrationCam::CalcAverageSigma()
     639{
     640 
     641  for (Int_t j=0; j<fGeom->GetNumAreas(); j++)
     642    {
     643 
     644      MCalibrationPix &pix    = (*fCam).GetAverageArea(j);
     645
     646      fAverageAreaSigma[j]    = pix.GetSigma    () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);
     647      fAverageAreaSigmaErr[j] = pix.GetSigmaErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);
     648
     649      pix.SetSigma   (fAverageAreaSigma[j]);
     650      pix.SetSigmaErr(fAverageAreaSigmaErr[j]);
     651
     652      fAverageAreaRelSigma[j]   = fAverageAreaSigma[j] / pix.GetMean();
     653     
     654      Float_t relsigmaerr       =  fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j]
     655                                / (fAverageAreaSigma[j]   *fAverageAreaSigma[j]   );
     656      relsigmaerr               += pix.GetMeanErr()*pix.GetMeanErr()
     657                                / (pix.GetMean()   *pix.GetMean()   );
     658      relsigmaerr               *= fAverageAreaRelSigma[j];
     659      fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr);
     660    }
     661}
     662
     663// ---------------------------------------------------------------------------
     664//
     665// Returns if the histogram is empty and sets the following flag:
     666// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
     667//
     668// Fits the histograms with a Gaussian, in case of failure
     669// calls MHGausEvents::RepeatFit(), in case of repeated failure
     670// calls MHGausEvents::BypassFit() and sets the following flags:
     671// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp )
     672// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
     673//
     674// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK().
     675// In case no, sets the following flags:
     676// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp )
     677// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
     678//
     679// Retrieves the results and store them in MCalibrationPix
     680//
     681void MHCalibrationCam::FitHiGainHists(MHGausEvents &hist,
     682                                      MCalibrationPix &pix,
     683                                      MBadPixelsPix &bad,
     684                                      MBadPixelsPix::UncalibratedType_t fittyp,
     685                                      MBadPixelsPix::UncalibratedType_t osctyp)
     686{
     687
     688  if (hist.IsEmpty())
     689    return;
     690 
     691  //
     692  // 2) Fit the Hi Gain histograms with a Gaussian
     693  //
     694  if (!hist.FitGaus())
     695    //
     696    // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
     697    //
     698    if (!hist.RepeatFit())
     699      {
     700        hist.BypassFit();
     701        bad.SetUncalibrated( fittyp );
     702      }
     703 
     704  //
     705  // 4) Check for oscillations
     706  //
     707  hist.CreateFourierSpectrum();
     708 
     709  if (!hist.IsFourierSpectrumOK())
     710    bad.SetUncalibrated( osctyp );
     711 
     712  //
     713  // 5) Retrieve the results and store them in this class
     714  //
     715  pix.SetHiGainMean      ( hist.GetMean()      );
     716  pix.SetHiGainMeanErr   ( hist.GetMeanErr()   );
     717  pix.SetHiGainSigma     ( hist.GetSigma()     );
     718  pix.SetHiGainSigmaErr  ( hist.GetSigmaErr()  );
     719  pix.SetHiGainProb      ( hist.GetProb()      );
     720  pix.SetHiGainNumPickup ( hist.GetPickup()    );
     721 
     722}
     723
     724
     725// ---------------------------------------------------------------------------
     726//
     727// Returns if the histogram is empty and sets the following flag:
     728// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
     729//
     730// Fits the histograms with a Gaussian, in case of failure
     731// calls MHGausEvents::RepeatFit(), in case of repeated failure
     732// calls MHGausEvents::BypassFit() and sets the following flags:
     733// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp )
     734// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
     735//
     736// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK().
     737// In case no, sets the following flags:
     738// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp )
     739// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
     740//
     741// Retrieves the results and store them in MCalibrationPix
     742//
     743void MHCalibrationCam::FitLoGainHists(MHGausEvents &hist,
     744                                      MCalibrationPix &pix,
     745                                      MBadPixelsPix &bad,
     746                                      MBadPixelsPix::UncalibratedType_t fittyp,
     747                                      MBadPixelsPix::UncalibratedType_t osctyp)
     748{
     749
     750  if (hist.IsEmpty())
     751      return;
     752 
     753
     754  //
     755  // 2) Fit the Hi Gain histograms with a Gaussian
     756  //
     757  if (!hist.FitGaus())
     758    //
     759    // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
     760    //
     761    if (!hist.RepeatFit())
     762      {
     763        hist.BypassFit();
     764        bad.SetUncalibrated( fittyp );
     765      }
     766 
     767  //
     768  // 4) Check for oscillations
     769  //
     770  hist.CreateFourierSpectrum();
     771 
     772  if (!hist.IsFourierSpectrumOK())
     773    bad.SetUncalibrated( osctyp );
     774 
     775  //
     776  // 5) Retrieve the results and store them in this class
     777  //
     778  pix.SetLoGainMean      ( hist.GetMean()      );
     779  pix.SetLoGainMeanErr   ( hist.GetMeanErr()   );
     780  pix.SetLoGainSigma     ( hist.GetSigma()     );
     781  pix.SetLoGainSigmaErr  ( hist.GetSigmaErr()  );
     782  pix.SetLoGainProb      ( hist.GetProb()      );
     783  pix.SetLoGainNumPickup ( hist.GetPickup()    );
     784 
     785}
     786
     787
     788
     789// --------------------------------------------------------------------------
     790//
    308791// Dummy, needed by MCamEvent
    309792//
     
    343826  for (Int_t i=0; i<nareas;i++)
    344827    {
     828
    345829      pad->cd(2*(i+1)-1);
    346830      GetAverageHiGainArea(i).Draw(opt);
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h

    r3631 r3636  
    2121#endif
    2222
     23#ifndef MARS_MBadPixelsPix
     24#include "MBadPixelsPix.h"
     25#endif
     26
    2327class TText;
    2428class TArrayI;
    2529class TArrayF;
    2630class MHGausEvents;
     31class MGeomCam;
     32class MCalibrationCam;
     33class MCalibrationPix;
     34class MBadPixelsCam;
     35class MBadPixelsPix;
    2736class MHCalibrationCam : public MH, public MCamEvent
    2837{
    2938 
     39private:
     40
     41  static const Int_t fgAverageNbins;     //! The default for fAverageNbins    (now set to: 2000)
     42  static const Int_t fgPulserFrequency;  //! The default for fPulserFrequency (now set to: 500)
     43
    3044protected:
    3145
    32   static const Int_t fgPulserFrequency;  // The default for fPulserFrequency (now set to: 500)
     46  Int_t   fAverageNbins;                // Number of bins for the average histograms
    3347  Int_t   fPulserFrequency;             // Light pulser frequency
    3448 
     
    4054  TObjArray *fAverageLoGainSectors;     //-> Array of calibration pixels, one per camera sector
    4155
     56  MGeomCam         *fGeom;              //!  Camera geometry
     57  MBadPixelsCam    *fBadPixels;         //!  Bad Pixels storage container
     58  MCalibrationCam  *fCam;               //!  Calibration Cam with the results
     59
    4260  TArrayI fAverageAreaNum;              // Number of pixels in average pixels per area
    4361  TArrayI fAverageAreaSat;              // Number of saturated slices in average pixels per area
     
    4866  TArrayI fAverageSectorNum;            // Number of pixels in average pixels per sector
    4967
     68  virtual Bool_t SetupHists(const MParList *pList);
     69  virtual Bool_t ReInitHists(MParList *pList); 
     70  virtual Bool_t FillHists(const MParContainer *par, const Stat_t w=1);
     71  virtual Bool_t FinalizeHists();
     72  virtual void   FinalizeBadPixels();
     73 
     74  void InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i);
     75
     76  void FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
     77                       MBadPixelsPix::UncalibratedType_t fittyp,
     78                       MBadPixelsPix::UncalibratedType_t osctyp);
     79 
     80  void FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam,
     81                       MBadPixelsPix::UncalibratedType_t fittyp,
     82                       MBadPixelsPix::UncalibratedType_t osctyp);
     83 
     84  void FitHiGainHists(MHGausEvents &hist,
     85                      MCalibrationPix &pix,
     86                      MBadPixelsPix &bad,
     87                      MBadPixelsPix::UncalibratedType_t fittyp,
     88                      MBadPixelsPix::UncalibratedType_t osctyp);
     89 
     90  void FitLoGainHists(MHGausEvents &hist,
     91                      MCalibrationPix &pix,
     92                      MBadPixelsPix &bad,
     93                      MBadPixelsPix::UncalibratedType_t fittyp,
     94                      MBadPixelsPix::UncalibratedType_t osctyp);
     95
     96  void CalcAverageSigma();
     97 
    5098  void DrawAverageSigma(Bool_t sat, Bool_t inner,
    5199                        Float_t sigma, Float_t sigmaerr,
     
    57105  ~MHCalibrationCam();
    58106
     107  void SetAverageNbins(   const Int_t bins=fgAverageNbins ) { fAverageNbins = bins; }
    59108  void SetPulserFrequency(const Int_t f=fgPulserFrequency) { fPulserFrequency = f; }
    60109 
     
    76125  MHGausEvents  &GetAverageLoGainSector(UInt_t i);
    77126  const MHGausEvents   &GetAverageLoGainSector(UInt_t i)  const;
     127
     128  virtual Bool_t SetupFill(const MParList *pList);
     129  virtual Bool_t ReInit   (      MParList *pList);
     130  virtual Bool_t Fill     (const MParContainer *par, const Stat_t w=1);
     131  virtual Bool_t Finalize ( );
    78132
    79133  // Clone
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc

    r3631 r3636  
    127127#include "MGeomPix.h"
    128128
     129#include "MHGausEvents.h"
     130
    129131#include "MBadPixelsCam.h"
    130132#include "MBadPixelsPix.h"
     
    140142using namespace std;
    141143
    142 const Int_t   MHCalibrationChargeCam::fgAverageNbins             = 4000;
    143144const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01;
    144145const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005;
     
    150151// - all pointers to NULL
    151152//
    152 // Initializes and sets owner of (done in MHCalibrationCam):
    153 // - fHiGainArray, fLoGainArray
    154 // - fAverageHiGainAreas, fAverageLoGainAreas
    155 // - fAverageHiGainSectors, fAverageLoGainSectors
    156 //
    157153// Initializes:
    158154// - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit
    159155// - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit
    160 // - fPulserFrequency to fgPulserFrequency (in MHCalibrationCam)
    161156//
    162157MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title)
    163     : fCam(NULL), fRawEvt(NULL), fGeom(NULL), fBadPixels(NULL)
     158    : fRawEvt(NULL)
    164159{
    165160    fName  = name  ? name  : "MHCalibrationChargeCam";
     
    174169// Gets the pointers to:
    175170// - MRawEvtData
    176 // - MGeomCam
    177 //
    178 // Calls Delete-Function of:
    179 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
    180 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
    181 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
    182 //
    183 Bool_t MHCalibrationChargeCam::SetupFill(const MParList *pList)
     171//
     172Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList)
    184173{
    185174 
     
    191180  }
    192181
    193   fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
    194   if (!fGeom)
    195   {
    196       *fLog << err << "MGeomCam not found... aborting." << endl;
    197       return kFALSE;
    198   }
    199 
    200   fHiGainArray->Delete();
    201   fLoGainArray->Delete();
    202 
    203   fAverageHiGainAreas->Delete();
    204   fAverageLoGainAreas->Delete();
    205 
    206   fAverageHiGainSectors->Delete();
    207   fAverageLoGainSectors->Delete();
    208 
    209182  return kTRUE;
    210183}
     
    213186//
    214187// Gets or creates the pointers to:
    215 // - MBadPixelsCam
    216188// - MCalibrationChargeCam
    217189//
     
    219191// - MExtractedSignalCam
    220192//
    221 // Initializes, if empty to MExtractedSignalCam::GetSize() for:
    222 // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
    223 //
    224 // Initializes, if empty to MGeomCam::GetNumAreas() for:
    225 // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
    226 //
    227 // Initializes, if empty to MGeomCam::GetNumSectors() for:
    228 // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
    229 //
    230 // Initializes TArrays to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively
    231 // Fills with number of valid pixels (if !MBadPixelsPix::IsBad()):
    232 // - MHCalibrationCam::fAverageAreaNum[area index]
    233 // - MHCalibrationCam::fAverageSectorNum[area index]
    234 //
    235 // Calls InitializeHiGainHists() and InitializeLoGainHists() for every entry in:
     193// Calls InitializeHists() for every entry in:
    236194// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
    237195// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     
    239197//
    240198// Sets Titles and Names for the Charge Histograms and
    241 // Sets number of bins to fAverageNbins for:
     199// Sets number of bins to MHCalibrationCam::fAverageNbins for:
    242200// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
    243201// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
    244202//
    245 Bool_t MHCalibrationChargeCam::ReInit(MParList *pList)
     203Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList)
    246204{
    247205
    248   fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam");
    249   if (!fBadPixels)
    250       return kFALSE;
    251 
    252 
    253   fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam");
     206  fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationChargeCam");
    254207  if (!fCam)
    255208      return kFALSE;
     
    262215  }
    263216
    264   const Int_t n        = signal->GetSize();
     217  const Int_t npixels  = fGeom->GetNumPixels();
    265218  const Int_t nsectors = fGeom->GetNumSectors();
    266219  const Int_t nareas   = fGeom->GetNumAreas();
    267220
    268   //
    269   // The function TArrayF::Set() already sets all entries to 0.
    270   //
    271   fAverageAreaNum.        Set(nareas);
    272   fAverageAreaSat.        Set(nareas);           
    273   fAverageAreaSigma.      Set(nareas);     
    274   fAverageAreaSigmaErr.   Set(nareas);   
    275   fAverageAreaRelSigma.   Set(nareas);   
    276   fAverageAreaRelSigmaErr.Set(nareas);
    277   fAverageSectorNum.      Set(nsectors);
    278 
    279   for (Int_t i=0; i<n; i++)
    280     {
    281       if ((*fBadPixels)[i].IsBad())
    282         continue;
    283       fAverageAreaNum  [(*fGeom)[i].GetAidx()  ]++;
    284       fAverageSectorNum[(*fGeom)[i].GetSector()]++;
    285     }
    286 
    287221  if (fHiGainArray->GetEntries()==0)
    288222  {
    289       fHiGainArray->Expand(n);
    290       for (Int_t i=0; i<n; i++)
     223      fHiGainArray->Expand(npixels);
     224      for (Int_t i=0; i<npixels; i++)
    291225      {
    292226          (*fHiGainArray)[i] = new MHCalibrationChargeHiGainPix;
    293           InitializeHiGainHists((MHCalibrationChargePix&)(*this)[i],(*fBadPixels)[i],i);
     227          InitHists((*this)[i],(*fBadPixels)[i],i);
    294228      }
    295229  }
    296230
    297  
    298231  if (fLoGainArray->GetEntries()==0)
    299232  {
    300       fLoGainArray->Expand(n);
    301      
    302       for (Int_t i=0; i<n; i++)
     233      fLoGainArray->Expand(npixels);
     234     
     235      for (Int_t i=0; i<npixels; i++)
    303236      {
    304237          (*fLoGainArray)[i] = new MHCalibrationChargeLoGainPix;
    305           InitializeLoGainHists((MHCalibrationChargePix&)(*this)(i),(*fBadPixels)[i],i);
     238          InitHists((*this)(i),(*fBadPixels)[i],i);
    306239      }
    307240     
     
    318251                                           "Average HiGain FADC sums area idx ");
    319252
     253        InitHists(GetAverageHiGainArea(j),fCam->GetAverageBadArea(j),j);
     254
    320255        MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainArea(j);
    321256
    322         InitializeHiGainHists(hist,fCam->GetAverageBadArea(j),j);
    323 
    324257        hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Area Idx ");
    325         hist.SetChargeNbins(fAverageNbins);
     258        hist.SetNbins(fAverageNbins);
    326259        hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx ");
    327260      }
     
    340273        MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainArea(j);
    341274
    342         InitializeLoGainHists(hist,fCam->GetAverageBadArea(j),j);
     275        InitHists(hist,fCam->GetAverageBadArea(j),j);
    343276
    344277        hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Area Idx ");
    345         hist.SetChargeNbins(fAverageNbins);
     278        hist.SetNbins(fAverageNbins);
    346279        hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx ");
    347280
     
    361294          MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainSector(j);
    362295
    363           InitializeHiGainHists(hist,fCam->GetAverageBadSector(j),j);
     296          InitHists(hist,fCam->GetAverageBadSector(j),j);
    364297         
    365298          hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector ");
    366           hist.SetChargeNbins(fAverageNbins);
     299          hist.SetNbins(fAverageNbins);
    367300          hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector ");
    368301
     
    382315          MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainSector(j);
    383316
    384           InitializeLoGainHists(hist,fCam->GetAverageBadSector(j),j);
     317          InitHists(hist,fCam->GetAverageBadSector(j),j);
    385318         
    386319          hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector ");
    387           hist.SetChargeNbins(fAverageNbins);
     320          hist.SetNbins(fAverageNbins);
    388321          hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector ");
    389322
     
    394327}
    395328
    396 // -------------------------------------------------------------
    397 //
    398 // If MBadPixelsPix::IsBad():
    399 // - calls MHCalibrationChargePix::SetExcluded()
    400 //
    401 // Calls:
    402 // - MHCalibrationChargePix::Init()
    403 // - MHCalibrationChargePix::ChangeHistId(i)
    404 // - MHGausEvents::SetEventFrequency(fPulserFrequency)
    405 //
    406 void MHCalibrationChargeCam::InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i)
    407 {
    408  
    409   if (bad.IsBad())
    410     {
    411       *fLog << warn << "Excluded pixel: " << i << " from calibration " << endl;
    412       hist.SetExcluded();
    413     }
    414  
    415   hist.Init();
    416   hist.ChangeHistId(i);
    417   hist.SetEventFrequency(fPulserFrequency);
    418          
    419 }
    420 
    421 // -------------------------------------------------------------
    422 //
    423 // If MBadPixelsPix::IsBad():
    424 // - calls MHCalibrationChargePix::SetExcluded()
    425 //
    426 // Calls:
    427 // - MHCalibrationChargePix::Init()
    428 // - MHCalibrationChargePix::ChangeHistId(i)
    429 // - MHGausEvents::SetEventFrequency(fPulserFrequency)
    430 //
    431 void MHCalibrationChargeCam::InitializeLoGainHists(MHCalibrationChargePix &hist,MBadPixelsPix &bad, const Int_t i)
    432 {
    433 
    434   if (bad.IsBad())
    435     hist.SetExcluded();
    436  
    437   hist.Init();
    438   hist.ChangeHistId(i);
    439   hist.SetEventFrequency(fPulserFrequency);
    440 }
    441 
    442329 
    443330// --------------------------------------------------------------------------
    444331//
    445332// Retrieves from MExtractedSignalCam:
     333// - first used LoGain FADC slice
     334//
     335// Retrieves from MGeomCam:
    446336// - number of pixels
    447 // - first used LoGain FADC slice
    448 //
    449 // Retrieves from MGeomCam:
    450337// - number of pixel areas
    451338// - number of sectors
     
    453340// For all TObjArray's (including the averaged ones), the following steps are performed:
    454341//
    455 // 1) Test size and return kFALSE if not matching
    456 //
    457 // 2) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with:
     342// 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with:
    458343// - MExtractedSignalPix::GetExtractedSignalHiGain();
    459344// - MExtractedSignalPix::GetExtractedSignalLoGain();
    460345//
    461 // 3) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:
     346// 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with:
    462347// - MExtractedSignalPix::GetNumHiGainSaturated();
    463348// - MExtractedSignalPix::GetNumLoGainSaturated();
    464349//
    465 // 4) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:
     350// 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with:
    466351// - MRawEvtPixelIter::GetIdxMaxHiGainSample();       
    467352// - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice);
    468353//
    469 Bool_t MHCalibrationChargeCam::Fill(const MParContainer *par, const Stat_t w)
     354Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w)
    470355{
    471356
     
    477362    }
    478363 
    479   const Int_t npixels  = signal->GetSize();
    480   const Int_t lofirst  = signal->GetFirstUsedSliceLoGain();
     364  const Int_t npixels  = fGeom->GetNumPixels();
    481365  const Int_t nareas   = fGeom->GetNumAreas();
    482366  const Int_t nsectors = fGeom->GetNumSectors();
    483 
    484   if (fHiGainArray->GetEntries() != npixels)
    485     {
    486       *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl;
    487       return kFALSE;
    488     }
    489 
    490   if (fLoGainArray->GetEntries() != npixels)
    491     {
    492       *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl;
    493       return kFALSE;
    494     }
    495 
    496   if (fAverageHiGainAreas->GetEntries() != nareas)
    497     {
    498       *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
    499       return kFALSE;
    500     }
    501 
    502   if (fAverageLoGainAreas->GetEntries() != nareas)
    503     {
    504       *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl;
    505       return kFALSE;
    506     }
    507 
    508   if (fAverageHiGainSectors->GetEntries() != nsectors)
    509     {
    510       *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
    511       return kFALSE;
    512     }
    513 
    514   if (fAverageLoGainSectors->GetEntries() != nsectors)
    515     {
    516       *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl;
    517       return kFALSE;
    518     }
     367  const Int_t lofirst  = signal->GetFirstUsedSliceLoGain();
    519368
    520369  Float_t sumhiarea  [nareas],   sumloarea  [nareas],   timehiarea  [nareas],   timeloarea  [nareas];
     
    526375    {
    527376      sumhiarea  [j] = sumloarea  [j] = timehiarea  [j] =  timeloarea  [j] = 0.;
     377      sathiarea  [j] = satloarea  [j] = 0;
     378    }
     379 
     380  for (UInt_t j=0; j<nsectors; j++)
     381    {
    528382      sumhisector[j] = sumlosector[j] = timehisector[j] =  timelosector[j] = 0.;
    529       sathiarea  [j] = satloarea  [j] = 0;
    530383      sathisector[j] = satlosector[j] = 0;
    531384    }
    532385 
     386
    533387  for (Int_t i=0; i<npixels; i++)
    534388    {
     
    596450    }
    597451 
    598 
    599452  for (UInt_t j=0; j<nareas; j++)
    600453    {
     
    642495// For all TObjArray's (including the averaged ones), the following steps are performed:
    643496//
    644 // 1) Return if the pixel is excluded
    645 //
    646 // 2) Call to FinalizeHiGainHists() and FinalizeLoGainHists()
    647 //
    648 // For all averaged areas, the fitted sigma is multiplied with the square root of
    649 // the number involved pixels
    650 //
    651 Bool_t MHCalibrationChargeCam::Finalize()
     497// 1) Returns if the pixel is excluded.
     498// 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation()
     499//    or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated )
     500// 3) Store the absolute arrival times in the MCalibrationChargePix's. If flag
     501//    MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored,
     502//    otherwise the Hi-Gain ones.
     503// 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays()
     504//    and sets the flags (if  occurring):
     505//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
     506//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
     507//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
     508//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
     509// 5) Retrieves the results and store them in MCalibrationChargePix
     510//
     511Bool_t MHCalibrationChargeCam::FinalizeHists()
    652512{
    653513
     
    656516     
    657517      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)(*this)[i];
    658 
     518      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)(*fCam)[i];
     519     
    659520      if (histhi.IsExcluded())
    660521        continue;
    661522     
    662       MCalibrationChargePix  &pix    = (*fCam)[i];
     523      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
     524        {
     525          pix.SetHiGainSaturation();
     526          histhi.CreateFourierSpectrum();
     527          continue;
     528        }
     529
     530      pix.SetAbsTimeMean ( histhi.GetAbsTimeMean());
     531      pix.SetAbsTimeRms  ( histhi.GetAbsTimeRms() );
     532    }
     533
     534  for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
     535    {
     536     
     537      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
    663538      MBadPixelsPix          &bad    = (*fBadPixels)[i];
    664      
    665       FinalizeHiGainHists(histhi,pix,bad);
    666      
    667     }
    668  
    669   for (Int_t i=0; i<fLoGainArray->GetSize(); i++)
    670     {
    671      
    672       MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i);
    673539
    674540      if (histlo.IsExcluded())
    675541        continue;
    676542     
    677       MCalibrationChargePix        &pix    = (*fCam)[i];
    678       MBadPixelsPix                &bad    = (*fBadPixels)[i];
    679      
    680       FinalizeLoGainHists(histlo,pix,bad);
    681      
    682     }
    683  
     543      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
     544        {
     545          *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl;
     546          bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
     547          histlo.CreateFourierSpectrum();
     548          continue;
     549        }
     550 
     551      MCalibrationChargePix &pix    = (MCalibrationChargePix&)(*fCam)[i];
     552     
     553      if (pix.IsHiGainSaturation())
     554        {
     555          pix.SetAbsTimeMean     ( histlo.GetAbsTimeMean());
     556          pix.SetAbsTimeRms      ( histlo.GetAbsTimeRms() );
     557        }           
     558    }
     559
    684560  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
    685561    {
    686562     
    687563      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j);     
    688       MCalibrationChargePix  &pix    = fCam->GetAverageArea(j);
    689       MBadPixelsPix          &bad    = fCam->GetAverageBadArea(j);       
    690      
    691       FinalizeHiGainHists(histhi,pix,bad);
     564      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageArea(j);
     565     
     566      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
     567        {
     568          pix.SetHiGainSaturation();
     569          histhi.CreateFourierSpectrum();
     570          continue;
     571        }
     572
     573      pix.SetAbsTimeMean ( histhi.GetAbsTimeMean());
     574      pix.SetAbsTimeRms  ( histhi.GetAbsTimeRms() );
    692575    }
    693576 
     
    696579     
    697580      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j);     
    698       MCalibrationChargePix  &pix    = fCam->GetAverageArea(j);
    699       MBadPixelsPix          &bad    = fCam->GetAverageBadArea(j);       
    700      
    701       FinalizeLoGainHists(histlo,pix,bad);
     581      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageArea(j);
     582     
     583      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
     584        {
     585          *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl;
     586          histlo.CreateFourierSpectrum();
     587          continue;
     588        }
     589
     590      if (pix.IsHiGainSaturation())
     591        {
     592          pix.SetAbsTimeMean  ( histlo.GetAbsTimeMean());
     593          pix.SetAbsTimeRms   ( histlo.GetAbsTimeRms() );
     594        }           
    702595    }
    703596
     
    706599     
    707600      MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j);     
    708       MCalibrationChargePix  &pix    = fCam->GetAverageSector(j);
     601      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageSector(j);
     602     
     603      if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries())
     604        {
     605          pix.SetHiGainSaturation();
     606          histhi.CreateFourierSpectrum();
     607          continue;
     608        }
     609
     610      pix.SetAbsTimeMean ( histhi.GetAbsTimeMean());
     611      pix.SetAbsTimeRms  ( histhi.GetAbsTimeRms() );
     612    }
     613 
     614  for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
     615    {
     616     
     617      MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j);     
     618      MCalibrationChargePix  &pix    = (MCalibrationChargePix&)fCam->GetAverageSector(j);
    709619      MBadPixelsPix          &bad    = fCam->GetAverageBadSector(j);       
    710620     
    711       FinalizeHiGainHists(histhi,pix,bad);
    712     }
    713  
    714   for (Int_t j=0; j<fAverageLoGainSectors->GetSize(); j++)
    715     {
    716      
    717       MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j);     
    718       MCalibrationChargePix  &pix    = fCam->GetAverageSector(j);
    719       MBadPixelsPix          &bad    = fCam->GetAverageBadSector(j);       
    720      
    721       FinalizeLoGainHists(histlo,pix,bad);
    722     }
    723  
    724   for (Int_t j=0; j<fGeom->GetNumAreas();j++)
    725     {
    726      
    727       MCalibrationChargePix &pix = fCam->GetAverageArea(j);
     621      if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries())
     622        {
     623          *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl;
     624          bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
     625          histlo.CreateFourierSpectrum();
     626          continue;
     627        }
    728628
    729629      if (pix.IsHiGainSaturation())
    730         fAverageAreaSat[j]++;
    731 
    732       fAverageAreaSigma[j]    = pix.GetSigmaCharge    () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);
    733       fAverageAreaSigmaErr[j] = pix.GetSigmaChargeErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]);
    734 
    735       pix.SetSigmaCharge   (fAverageAreaSigma[j]);
    736       pix.SetSigmaChargeErr(fAverageAreaSigmaErr[j]);
    737 
    738       fAverageAreaRelSigma[j]   = fAverageAreaSigma[j] / pix.GetMeanCharge();
    739      
    740       Float_t relsigmaerr       =  fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j]
    741                                 / (fAverageAreaSigma[j]   *fAverageAreaSigma[j]   );
    742       relsigmaerr               += pix.GetMeanChargeErr()*pix.GetMeanChargeErr()
    743                                 / (pix.GetMeanCharge()   *pix.GetMeanCharge()   );
    744       relsigmaerr               *= fAverageAreaRelSigma[j];
    745       fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr);
    746 
    747     }
    748 
     630        {
     631          pix.SetAbsTimeMean  ( histlo.GetAbsTimeMean());
     632          pix.SetAbsTimeRms   ( histlo.GetAbsTimeRms() );
     633        }           
     634    }
     635 
     636  //
     637  // Perform the fitting for the High Gain (done in MHCalibrationCam)
     638  //
     639  FitHiGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
     640                  MBadPixelsPix::kHiGainNotFitted,
     641                  MBadPixelsPix::kHiGainOscillating);
     642  //
     643  // Perform the fitting for the Low Gain (done in MHCalibrationCam)
     644  //
     645  FitLoGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels),
     646                  MBadPixelsPix::kLoGainNotFitted,
     647                  MBadPixelsPix::kLoGainOscillating);
     648     
    749649  return kTRUE;
    750650}
    751651
    752 // ---------------------------------------------------------------------------
    753 //
    754 // Returns if the histogram is empty and sets the following flag:
    755 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
    756 //
    757 // Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is
    758 // higher than the allowed limit (fNumHiGainSaturationLimit*histogram-entries), creates
    759 // the fourier spectrum and sets the following flags:
    760 // - MCalibrationChargePix::SetHiGainSaturation()
    761 //
    762 // Fits the histograms with a Gaussian, in case of failure
    763 // calls MHCalibrationChargePix::RepeatFit(), in case of repeated failure
    764 // calls MHGausEvents::BypassFit() and sets the following flags:
    765 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
    766 // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
    767 //
    768 // Counts the number of pickup events
    769 //
    770 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK().
    771 // In case no, sets the following flags:
    772 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
    773 // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
    774 //
    775 // Retrieves the results and store them in MCalibrationChargePix
    776 //
    777 void MHCalibrationChargeCam::FinalizeHiGainHists(MHCalibrationChargePix &hist,
    778                                                  MCalibrationChargePix &pix,
    779                                                  MBadPixelsPix &bad)
     652// --------------------------------------------------------------------------
     653//
     654// Takes the decisions under which conditions a pixel is declared:
     655// MBadPixelsPix::kUnsuitableRun  or MBadPixelsPix::kUnreliableRun, namely:
     656// * if MBadPixelsPix::kHiGainNotFitted   and !MCalibrationPix::IsHiGainSaturation()
     657//   sets MBadPixelsPix::kUnreliableRun
     658// * if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation()
     659//   sets MBadPixelsPix::kUnreliableRun
     660// * if MBadPixelsPix::kLoGainNotFitted   and  MCalibrationPix::IsLoGainSaturation()
     661//   sets MBadPixelsPix::kUnreliableRun
     662// * if MBadPixelsPix::kLoGainOscillating and  MCalibrationPix::IsLoGainSaturation()
     663//   sets MBadPixelsPix::kUnreliableRun
     664// * if MBadPixelsPix::kLoGainSaturation
     665//   sets MBadPixelsPix::kUnsuitableRun
     666//
     667void MHCalibrationChargeCam::FinalizeBadPixels()
    780668{
    781 
    782     if (hist.IsEmpty())
    783       {
    784         *fLog << warn << "Empty Hi Gain histogram in pixel: " << pix.GetPixId() << endl;
    785         bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);       
    786         return;
    787       }
    788    
    789     if (hist.GetSaturated() > fNumHiGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    790     {
    791       *fLog << warn << "Saturated Hi Gain histogram in pixel: " << pix.GetPixId() << endl;
    792       pix.SetHiGainSaturation();
    793       hist.CreateFourierSpectrum();
    794       return;
    795     }
    796 
    797     //
    798     // 2) Fit the Hi Gain histograms with a Gaussian
    799     //
    800     if (!hist.FitGaus())
    801     //
    802     // 3) In case of failure set the bit Fitted to false and take histogram means and RMS
    803     //
    804       if (!hist.RepeatFit())
    805         {
    806           hist.BypassFit();
    807           *fLog << warn << "Hi Gain could not be fitted in pixel: " << pix.GetPixId() << endl;
    808           bad.SetUncalibrated( MBadPixelsPix::kHiGainNotFitted );
    809           bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun   );
    810         }
    811 
    812     //
    813     // 4) Check for pickup
    814     //
    815     hist.CountPickup();
    816    
    817     //
    818     // 5) Check for oscillations
    819     //
    820     hist.CreateFourierSpectrum();
    821    
    822     if (!hist.IsFourierSpectrumOK())
    823       {
    824         *fLog << warn << "Oscillating Hi Gain in pixel: " << pix.GetPixId() << endl;
    825         bad.SetUncalibrated( MBadPixelsPix::kHiGainOscillating );
    826         bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun     );
    827       }
    828 
    829     //
    830     // 6) Retrieve the results and store them in this class
    831     //
    832     pix.SetHiGainMeanCharge(     hist.GetMean()      );
    833     pix.SetHiGainMeanChargeErr(  hist.GetMeanErr()   );
    834     pix.SetHiGainSigmaCharge(    hist.GetSigma()     );
    835     pix.SetHiGainSigmaChargeErr( hist.GetSigmaErr()  );
    836     pix.SetHiGainChargeProb    ( hist.GetProb()      );
    837    
    838     pix.SetAbsTimeMean         ( hist.GetAbsTimeMean());
    839     pix.SetAbsTimeRms          ( hist.GetAbsTimeRms() );
    840    
    841     pix.SetHiGainNumPickup     ( hist.GetPickup()     );
    842    
     669     
     670  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
     671    {
     672     
     673      MBadPixelsPix    &bad    = (*fBadPixels)[i];
     674      MCalibrationPix  &pix    = (*fCam)[i];
     675
     676      if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted ))
     677        if (!pix.IsHiGainSaturation())
     678          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
     679 
     680      if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating ))
     681        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
     682
     683      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted ))
     684        if (pix.IsHiGainSaturation())
     685          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
     686 
     687      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating ))
     688        if (pix.IsHiGainSaturation())
     689          bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
     690
     691      if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation ))
     692          bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
     693    }
    843694}
    844 
    845 // ---------------------------------------------------------------------------
    846 //
    847 // Returns if the histogram is empty. If yes and the flag MCalibrationChargePix::IsHiGainSaturation()
    848 // is set, sets the following flag:
    849 // - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun)
    850 //
    851 // Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is
    852 // higher than the allowed limit (fNumLoGainSaturationLimit*histogram-entries), creates
    853 // the fourier spectrum and sets the following flags:
    854 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
    855 // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
    856 //
    857 // Fits the histograms with a Gaussian, in case of failure calls MHCalibrationChargePix::RepeatFit(),
    858 // in case of repeated failure calls MHGausEvents::BypassFit() and sets the following flags:
    859 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
    860 //   If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:
    861 // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   )
    862 //
    863 // Counts the number of pickup events
    864 //
    865 // Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK().
    866 // In case no, sets the following flags:
    867 // - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
    868 //   If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally:
    869 // - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
    870 //
    871 // Retrieves the results and store them in MCalibrationChargePix
    872 //
    873 void MHCalibrationChargeCam::FinalizeLoGainHists(MHCalibrationChargePix &hist,
    874                                                  MCalibrationChargePix &pix,
    875                                                  MBadPixelsPix &bad)
    876 {
    877 
    878     if (hist.IsEmpty())
    879       {
    880         if (pix.IsHiGainSaturation())
    881           bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    882         return;
    883       }
    884 
    885     if (hist.GetSaturated() > fNumLoGainSaturationLimit*hist.GetHGausHist()->GetEntries())
    886     {
    887       *fLog << warn << "Saturated Lo Gain histogram in pixel: " << pix.GetPixId() << endl;
    888       bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation );
    889       bad.SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    );
    890       hist.CreateFourierSpectrum();
    891       return;
    892     }
    893 
    894     //
    895     // 2) Fit the Lo Gain histograms with a Gaussian
    896     //
    897     if (!hist.FitGaus())
    898     //
    899     // 3) In case of failure set the bit kFitted to false and take histogram means and RMS
    900     //
    901       if (!hist.RepeatFit())
    902         {
    903           hist.BypassFit();
    904           bad.SetUncalibrated( MBadPixelsPix::kLoGainNotFitted );
    905           if (pix.IsHiGainSaturation())
    906             bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
    907         }
    908 
    909     //
    910     // 4) Check for pickup
    911     //
    912     hist.CountPickup();
    913    
    914     //
    915     // 5) Check for oscillations
    916     //
    917     hist.CreateFourierSpectrum();
    918    
    919     if (!hist.IsFourierSpectrumOK())
    920       {
    921         bad.SetUncalibrated( MBadPixelsPix::kLoGainOscillating );
    922         if (pix.IsHiGainSaturation())
    923           bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun);
    924       }
    925     //
    926     // 6) Retrieve the results and store them in this class
    927     //
    928     pix.SetLoGainMeanCharge(     hist.GetMean()     );
    929     pix.SetLoGainMeanChargeErr(  hist.GetMeanErr()  );
    930     pix.SetLoGainSigmaCharge(    hist.GetSigma()    );
    931     pix.SetLoGainSigmaChargeErr( hist.GetSigmaErr() );
    932     pix.SetLoGainChargeProb    ( hist.GetProb()     );
    933    
    934     if (pix.IsHiGainSaturation())
    935     {
    936         pix.SetAbsTimeMean     ( hist.GetAbsTimeMean());
    937         pix.SetAbsTimeRms      ( hist.GetAbsTimeRms() );
    938     }       
    939    
    940     pix.SetLoGainNumPickup     (  hist.GetPickup()    );
    941 
    942 }   
    943695
    944696// --------------------------------------------------------------------------
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h

    r3631 r3636  
    77
    88class MRawEvtData;
    9 class MGeomCam;
    10 class MBadPixelsCam;
    119class MBadPixelsPix;
    1210class MCalibrationChargeCam;
     
    1715private:
    1816
    19   static const Int_t   fgAverageNbins;               // The default for fAverageNbins
    2017  static const Float_t fgNumHiGainSaturationLimit;   // The default for fNumHiGainSaturationLimit
    2118  static const Float_t fgNumLoGainSaturationLimit;   // The default for fNumLoGainSaturationLimit
    2219 
    23   Int_t   fAverageNbins;              // Number of bins for the average histograms
    2420  Float_t fNumHiGainSaturationLimit;  // Rel. amount sat. higain FADC slices until pixel is called saturated
    2521  Float_t fNumLoGainSaturationLimit;  // Rel. amount sat. logain FADC slices until pixel is called saturated
    2622 
    27   MCalibrationChargeCam  *fCam;       //!  Calibration Cam with the results
    28   MRawEvtData            *fRawEvt;    //!  Raw event data
    29   MGeomCam               *fGeom;      //!  Camera geometry
    30   MBadPixelsCam          *fBadPixels; //!  Bad Pixels storage container
     23  MRawEvtData *fRawEvt;               //!  Raw event data
    3124
    32   void InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i);
    33   void InitializeLoGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i); 
    34  
     25  Bool_t SetupHists(const MParList *pList);
     26  Bool_t ReInitHists(MParList *pList);
     27  Bool_t FillHists(const MParContainer *par, const Stat_t w=1);
     28  Bool_t FinalizeHists();
     29  void   FinalizeBadPixels();
     30
    3531  void FinalizeHiGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad);
    3632  void FinalizeLoGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad);
     
    4137  ~MHCalibrationChargeCam() {}
    4238 
    43 
    44   void SetAverageNbins(    const Int_t bins=fgAverageNbins  )          { fAverageNbins = bins; }
    4539  void SetNumLoGainSaturationLimit( const Float_t lim=fgNumLoGainSaturationLimit) { fNumLoGainSaturationLimit = lim; }
    4640  void SetNumHiGainSaturationLimit( const Float_t lim=fgNumHiGainSaturationLimit) { fNumHiGainSaturationLimit = lim; }
     
    4943  Float_t GetNumLoGainSaturationLimit()      const  { return fNumLoGainSaturationLimit; }
    5044
    51   Bool_t SetupFill(const MParList *pList);
    52   Bool_t ReInit   (      MParList *pList);
    53   Bool_t Fill     (const MParContainer *par, const Stat_t w=1);
    54   Bool_t Finalize ( );
    55 
    5645  Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
    5746  void   DrawPixelContent( Int_t num )  const;   
    5847
    59   ClassDef(MHCalibrationChargeCam, 1)   // Histogram class for camera calibration
     48  ClassDef(MHCalibrationChargeCam, 1)   // Histogram class for camera charge calibration
    6049};
    6150
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.cc

    r3613 r3636  
    6060//
    6161// Sets:
    62 // - the default number for fChargeNbins  (fgChargeNbins)
    63 // - the default number for fChargeFirst  (fgChargeFirst)
    64 // - the default number for fChargeLast   (fgChargeLast)
     62// - the default number for fNbins        (fgChargeNbins)
     63// - the default number for fFirst        (fgChargeFirst)
     64// - the default number for fLast         (fgChargeLast)
    6565// - the default number for fAbsTimeNbins (fgAbstTimeNbins)
    6666// - the default number for fAbsTimeFirst (fgAbsTimeFirst)
     
    7979  fTitle = title ? title : "Fill the FADC sums of the HiGainPix events and perform the fits Pixel ";
    8080 
    81   SetChargeNbins();
    82   SetChargeFirst();
    83   SetChargeLast();
     81  SetNbins ( fgChargeNbins );
     82  SetFirst ( fgChargeFirst );
     83  SetLast  ( fgChargeLast  );
    8484 
    8585  SetAbsTimeNbins();
     
    8787  SetAbsTimeLast();
    8888
    89   fHGausHist.SetName("HCalibrationChargeHiGainPix");
     89  fHGausHist.SetName ("HCalibrationChargeHiGainPix");
    9090  fHGausHist.SetTitle("Distribution of Summed Hi Gain FADC slices Pixel "); 
    9191
    92   fHAbsTime.SetName("HAbsTimeHiGainPix");
     92  fHAbsTime.SetName ("HAbsTimeHiGainPix");
    9393  fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Hi Gain Pixel "); 
    9494}
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.h

    r3613 r3636  
    1212private:
    1313
    14   static const Int_t   fgChargeNbins;     // Default for fChargeNBins (now set to: 2000  )
    15   static const Axis_t  fgChargeFirst;     // Default for fChargeFirst (now set to: -0.5  )
    16   static const Axis_t  fgChargeLast;      // Default for fChargeLast   (now set to: 1999.5)
    17   static const Int_t   fgAbsTimeNbins;    // Default for fAbsTimeNbins (now set to: 20    )
    18   static const Axis_t  fgAbsTimeFirst;    // Default for fAbsTimeFirst (now set to: -0.5  )
    19   static const Axis_t  fgAbsTimeLast;     // Default for fAbsTimeLast  (now set to: 19.5  )
     14  static const Int_t   fgChargeNbins;     // Default for MHGausEvents::fNbins (now set to: 2000  )
     15  static const Axis_t  fgChargeFirst;     // Default for MHGausEvents::fFirst (now set to: -0.5  )
     16  static const Axis_t  fgChargeLast;      // Default for MHGausEvents::fLast  (now set to: 1999.5)
     17  static const Int_t   fgAbsTimeNbins;    // Default for fAbsTimeNbins        (now set to: 20    )
     18  static const Axis_t  fgAbsTimeFirst;    // Default for fAbsTimeFirst        (now set to: -0.5  )
     19  static const Axis_t  fgAbsTimeLast;     // Default for fAbsTimeLast         (now set to: 19.5  )
    2020
    2121public:
     
    2525
    2626  // Setters
    27   void SetChargeNbins(const Int_t  bins =fgChargeNbins)          { fChargeNbins = bins;     }
    28   void SetChargeFirst(const Axis_t first=fgChargeFirst)          { fChargeFirst = first;    }
    29   void SetChargeLast (const Axis_t last =fgChargeLast)           { fChargeLast  = last;     }
    30  
    3127  void SetAbsTimeNbins(const Int_t  bins =fgAbsTimeNbins)        { fAbsTimeNbins = bins;    }
    3228  void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst)        { fAbsTimeFirst = first;   }
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.cc

    r3614 r3636  
    6060//
    6161// Sets:
    62 // - the default number for fChargeNbins  (fgChargeNbins)
    63 // - the default number for fChargeFirst  (fgChargeFirst)
    64 // - the default number for fChargeLast   (fgChargeLast)
     62// - the default number for fNbins        (fgChargeNbins)
     63// - the default number for fFirst        (fgChargeFirst)
     64// - the default number for fLast         (fgChargeLast)
    6565// - the default number for fAbsTimeNbins (fgAbstTimeNbins)
    6666// - the default number for fAbsTimeFirst (fgAbsTimeFirst)
     
    7979  fTitle = title ? title : "Fill the FADC sums of the Low Gain events and perform the fits Pixel ";
    8080 
    81   SetChargeNbins();
    82   SetChargeFirst();
    83   SetChargeLast();
     81  SetNbins ( fgChargeNbins );
     82  SetFirst ( fgChargeFirst );
     83  SetLast  ( fgChargeLast  );
    8484 
    8585  SetAbsTimeNbins();
     
    8787  SetAbsTimeLast();
    8888
    89   fHGausHist.SetName("HCalibrationChargeLoGainPix");
     89  fHGausHist.SetName ("HCalibrationChargeLoGainPix");
    9090  fHGausHist.SetTitle("Distribution of Summed Lo Gain FADC slices Pixel "); 
    9191
    92   fHAbsTime.SetName("HAbsTimeLoGainPix");
     92  fHAbsTime.SetName ("HAbsTimeLoGainPix");
    9393  fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Lo Gain Pixel "); 
    9494}
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.h

    r3614 r3636  
    1111private:
    1212
    13   static const Int_t   fgChargeNbins;      // Default for fChargeNBins (now set to: 200   )
    14   static const Axis_t  fgChargeFirst;      // Default for fChargeFirst (now set to: -0.5  )
    15   static const Axis_t  fgChargeLast;       // Default for fChargeLast   (now set to: 199.5 )
    16   static const Int_t   fgAbsTimeNbins;     // Default for fAbsTimeNbins (now set to: 15    )
    17   static const Axis_t  fgAbsTimeFirst;     // Default for fAbsTimeFirst (now set to: -0.5  )
    18   static const Axis_t  fgAbsTimeLast;      // Default for fAbsTimeLast  (now set to: 14.5  )
     13  static const Int_t   fgChargeNbins;      // Default for MHGausEvents::fNbins (now set to: 200   )
     14  static const Axis_t  fgChargeFirst;      // Default for MHGausEvents::fFirst (now set to: -0.5  )
     15  static const Axis_t  fgChargeLast;       // Default for MHGausEvents::fLast  (now set to: 199.5 )
     16  static const Int_t   fgAbsTimeNbins;     // Default for fAbsTimeNbins        (now set to: 15    )
     17  static const Axis_t  fgAbsTimeFirst;     // Default for fAbsTimeFirst        (now set to: -0.5  )
     18  static const Axis_t  fgAbsTimeLast;      // Default for fAbsTimeLast         (now set to: 14.5  )
    1919
    2020public:
     
    2424
    2525  // Setters
    26   void SetChargeNbins(const Int_t  bins =fgChargeNbins)          { fChargeNbins = bins;     }
    27   void SetChargeFirst(const Axis_t first=fgChargeFirst)          { fChargeFirst = first;    }
    28   void SetChargeLast (const Axis_t last =fgChargeLast)           { fChargeLast  = last;     }
    29  
    3026  void SetAbsTimeNbins(const Int_t  bins =fgAbsTimeNbins)        { fAbsTimeNbins = bins;    }
    3127  void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst)        { fAbsTimeFirst = first;   }
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc

    r3625 r3636  
    5656const Axis_t  MHCalibrationChargePix::fgAbsTimeFirst    = -0.5;
    5757const Axis_t  MHCalibrationChargePix::fgAbsTimeLast     = 14.5;
    58 const Float_t MHCalibrationChargePix::fgPickupLimit     = 5.;
    5958// --------------------------------------------------------------------------
    6059//
     
    6261//
    6362// Sets:
    64 // - the default number for fChargeNbins  (fgChargeNbins)
    65 // - the default number for fChargeFirst  (fgChargeFirst)
    66 // - the default number for fChargeLast   (fgChargeLast)
     63// - the default number for fNbins        (fgChargeNbins)
     64// - the default number for fFirst        (fgChargeFirst)
     65// - the default number for fLast         (fgChargeLast)
    6766// - the default number for fAbsTimeNbins (fgAbstTimeNbins)
    6867// - the default number for fAbsTimeFirst (fgAbsTimeFirst)
    6968// - the default number for fAbsTimeLast  (fgAbsTimeLast)
    70 // - the default number for fPickupLimit  (fgPickupLimit)
    7169//
    7270// - the default name of the  fHGausHist ("HCalibrationCharge")
     
    8886//
    8987MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title)
    90     : fHAbsTime(), fFlags(0)
     88    : fHAbsTime()
    9189{
    9290 
     
    9492  fTitle = title ? title : "Statistics of the FADC sums of calibration events";
    9593
    96   SetChargeNbins();
    97   SetChargeFirst();
    98   SetChargeLast();
     94  SetNbins ( fgChargeNbins );
     95  SetFirst ( fgChargeFirst );
     96  SetLast  ( fgChargeLast  );
    9997
    10098  SetAbsTimeNbins();
    10199  SetAbsTimeFirst();
    102100  SetAbsTimeLast();
    103 
    104   SetPickupLimit();
    105101
    106102  fHGausHist.SetName("HCalibrationCharge");
     
    124120//
    125121// Sets:
    126 // - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
     122// - fHGausHist.SetBins(fNbins,fFirst,fLast);
    127123// - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);
    128124//
    129 void MHCalibrationChargePix::Init()
    130 {
    131 
    132   fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
     125void MHCalibrationChargePix::InitBins()
     126{
     127  MHGausEvents::InitBins();
    133128  fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast);
    134129}
     
    137132//
    138133// Sets:
    139 // - all variables to 0., except fPixId to -1
    140 // - all flags to kFALSE
     134// - all variables to 0.
    141135//
    142136// Executes:
     
    146140{
    147141
    148   fPixId     = -1;
    149142  fSaturated = 0.;
    150   fPickup    = 0.;
    151 
    152   SetExcluded ( kFALSE );
    153  
     143
    154144  MHGausEvents::Clear();
    155145  return;
     
    166156// --------------------------------------------------------------------------
    167157//
    168 // Set fPixId to id
     158// Calls MHGausEvents::ChangeHistId()
    169159//
    170160// Add id to names and titles of:
    171 // - this
    172 // - fHGausHist
    173161// - fHAbsTime
    174162//
     
    176164{
    177165
    178   fPixId = id;
    179 
    180   fHGausHist.SetName(Form("%s%d", fHGausHist.GetName(), id));
    181   fHGausHist.SetTitle(Form("%s%d", fHGausHist.GetTitle(), id));
    182 
    183   fHAbsTime.SetName(Form("%s%d", fHAbsTime.GetName(), id));
     166  MHGausEvents::ChangeHistId(id);
     167
     168  fHAbsTime.SetName (Form("%s%d", fHAbsTime.GetName(),  id));
    184169  fHAbsTime.SetTitle(Form("%s%d", fHAbsTime.GetTitle(), id));
    185170
    186   fName  = Form("%s%d", fName.Data(), id);
    187   fTitle = Form("%s%d", fTitle.Data(), id);
    188 }
    189 
    190 // --------------------------------------------------------------------------
    191 //
    192 // Set Excluded bit from outside
    193 //
    194 void MHCalibrationChargePix::SetExcluded(const Bool_t b)
    195 {
    196     b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
    197 }
    198 
    199 const Bool_t MHCalibrationChargePix::IsExcluded()      const
    200 {
    201   return TESTBIT(fFlags,kExcluded);
    202171}
    203172
     
    317286
    318287}
    319 
    320 // -----------------------------------------------------------------------------
    321 //
    322 // Repeats the Gauss fit in a smaller range, defined by:
    323 //
    324 // min = GetMean() - fPickupLimit * GetSigma();
    325 // max = GetMean() + fPickupLimit * GetSigma();
    326 //
    327 Bool_t MHCalibrationChargePix::RepeatFit(const Option_t *option)
    328 {
    329 
    330   //
    331   // Get new fitting ranges
    332   //
    333   Axis_t rmin = GetMean() - fPickupLimit * GetSigma();
    334   Axis_t rmax = GetMean() + fPickupLimit * GetSigma();
    335 
    336   GetFGausFit()->SetRange(rmin,rmax);
    337 
    338   GetHGausHist()->Fit(GetFGausFit(),option);
    339 
    340   SetMean     ( GetFGausFit()->GetParameter(1) );
    341   SetMeanErr  ( GetFGausFit()->GetParameter(2) );
    342   SetSigma    ( GetFGausFit()->GetParError(1)  );
    343   SetSigmaErr ( GetFGausFit()->GetParError(2)  );
    344   SetProb     ( GetFGausFit()->GetProb()       );     
    345 
    346   //
    347   // The fit result is accepted under condition:
    348   // 1) The results are not nan's
    349   // 2) The NDF is not smaller than fNDFLimit (5)
    350   // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
    351   //
    352   if (   TMath::IsNaN ( GetMean()     )
    353       || TMath::IsNaN ( GetMeanErr()  )
    354       || TMath::IsNaN ( GetProb()     )   
    355       || TMath::IsNaN ( GetSigma()    )
    356       || TMath::IsNaN ( GetSigmaErr() )
    357       || GetFGausFit()->GetNDF() < fNDFLimit
    358       || GetProb() < fProbLimit )
    359     return kFALSE;
    360  
    361   SetGausFitOK(kTRUE);
    362   return kTRUE;
    363 
    364 }
    365 
    366 
    367 
    368 void MHCalibrationChargePix::CountPickup()
    369 {
    370    fPickup  = (Int_t)GetHGausHist()->Integral(GetHGausHist()->GetXaxis()->FindBin(GetMean()+fPickupLimit*GetSigma()),
    371                                                GetHGausHist()->GetXaxis()->GetLast(),
    372                                                "width");
    373 }
    374 
    375 
    376 
    377 
    378 
    379 
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h

    r3631 r3636  
    1212private:
    1313
    14   static const Int_t   fgChargeNbins;        // Default for fChargeNBins  (now set to: 2000  )
    15   static const Axis_t  fgChargeFirst;        // Default for fChargeFirst  (now set to: -0.5  )
    16   static const Axis_t  fgChargeLast;         // Default for fChargeLast   (now set to: 1999.5)
    17   static const Int_t   fgAbsTimeNbins;       // Default for fAbsTimeNbins (now set to: 15    )
    18   static const Axis_t  fgAbsTimeFirst;       // Default for fAbsTimeFirst (now set to: -0.5  )
    19   static const Axis_t  fgAbsTimeLast;        // Default for fAbsTimeLast  (now set to: 14.5  )
    20   static const Float_t fgPickupLimit;        // Default for fPickupLimit  (now set to:  5.   )
     14  static const Int_t   fgChargeNbins;        // Default for MHGausEvents::fNBins (now set to: 2000  )
     15  static const Axis_t  fgChargeFirst;        // Default for MHGausEvents::fFirst (now set to: -0.5  )
     16  static const Axis_t  fgChargeLast;         // Default for MHGausEvents::fLast  (now set to: 1999.5)
     17  static const Int_t   fgAbsTimeNbins;       // Default for fAbsTimeNbins        (now set to: 15    )
     18  static const Axis_t  fgAbsTimeFirst;       // Default for fAbsTimeFirst        (now set to: -0.5  )
     19  static const Axis_t  fgAbsTimeLast;        // Default for fAbsTimeLast         (now set to: 14.5  )
    2120
    2221protected:
    2322
    24   Int_t    fPixId;         // The pixel ID
    25 
    2623  TH1F     fHAbsTime;      // Histogram containing the absolute arrival times
    2724                         
    28   Int_t    fChargeNbins;   // Number of  bins used for the fHGausHist
    29   Axis_t   fChargeFirst;   // Lower bound bin used for the fHGausHist
    30   Axis_t   fChargeLast;    // Upper bound bin used for the fHGausHist
    31 
    3225  Int_t    fAbsTimeNbins;  // Number of  bins used for the fHAbsTime
    3326  Axis_t   fAbsTimeFirst;  // Lower bound bin used for the fHAbsTime
    3427  Axis_t   fAbsTimeLast;   // Upper bound bin used for the fHAbsTime
    3528
    36   Float_t  fPickupLimit;   // Upper number of sigmas from the fitted mean above which events are considered as pickup
    37 
    3829  Float_t  fSaturated;     // Number of events classified as saturated
    39   Float_t  fPickup;        // Number of events classified as pick-up
    40 
    41   Byte_t   fFlags;         // Bit-field for the flags
    42   enum     { kExcluded };  // Possible bits to be set
    4330
    4431public:
     
    4936  virtual void Clear(Option_t *o="");
    5037  virtual void Reset(); 
    51   virtual void Init();
    52   virtual void ChangeHistId(Int_t i);
     38  virtual void InitBins();
    5339 
    5440  // Setters
    55   virtual void SetChargeNbins(const Int_t  bins =fgChargeNbins)    { fChargeNbins = bins;   }
    56   virtual void SetChargeFirst(const Axis_t first=fgChargeFirst)    { fChargeFirst = first;  }
    57   virtual void SetChargeLast( const Axis_t last =fgChargeLast)     { fChargeLast  = last;   }
    58  
    5941  virtual void SetAbsTimeNbins(const Int_t  bins =fgAbsTimeNbins)  { fAbsTimeNbins = bins;  }
    6042  virtual void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst)  { fAbsTimeFirst = first; }
    6143  virtual void SetAbsTimeLast( const Axis_t last =fgAbsTimeLast)   { fAbsTimeLast  = last;  }
    6244
    63   virtual void SetPickupLimit( const Float_t  lim =fgPickupLimit)  { fPickupLimit  = lim;   }
    64 
    6545  void SetSaturated          ( const Float_t f    )                { fSaturated += f;      }
    66   void SetExcluded           ( const Bool_t b=kTRUE );
    6746
    6847  // Getters
     
    7049  const TH1F *GetHAbsTime()             const { return &fHAbsTime;  }
    7150
    72   const Float_t  GetIntegral()          const;
    73  
    7451  const Float_t  GetAbsTimeMean(  )     const;
    7552  const Float_t  GetAbsTimeRms()        const;
    76 
     53  const Float_t  GetIntegral()          const;
    7754  const Float_t  GetSaturated()         const { return fSaturated;   }
    78   const Float_t  GetPickup()            const { return fPickup;      }
    79 
    80   const Bool_t   IsExcluded()          const;
    8155
    8256  // Fill histos
    8357  Bool_t FillAbsTime(const Float_t t);
    8458
    85   // Fits
    86   Bool_t RepeatFit(const Option_t *option="RQ0");
    87  
    8859  // Draws
    8960  virtual void Draw(Option_t *opt="");
    9061
    9162  // Miscelleaneous
    92   void CountPickup();
    93 
     63  void ChangeHistId(Int_t id);
     64 
    9465  ClassDef(MHCalibrationChargePix, 1)     // Base Histogram class for a Charge Calibration Pixel
    9566};
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc

    r3625 r3636  
    2222!
    2323\* ======================================================================== */
    24 
    2524/////////////////////////////////////////////////////////////////////////////
    2625//                                                                         //
    27 // MHCalibrationRelTimeCam                                                            //
     26// MHCalibrationRelTimeCam                                                 //
    2827//                                                                         //
    29 // Hold the CalibrationRelTime information for all pixels in the camera              //
    30 //                                                                         //
     28// Fills the extracted relative arrival times of MArrivalTimeCam into
     29// the MHGausEvents-classes MHCalibrationRelTimePix for every:
     30//
     31// - pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
     32//
     33// - average pixel per area index (e.g. inner and outer for the MAGIC camera),
     34//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas
     35//
     36// - average pixel per camera sector (e.g. sectors 1-6 for the MAGIC camera),
     37//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
     38//
     39// Every time is filled into a histogram and an array, in order to perform
     40// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
     41// event-by-event basis and written into the corresponding average pixels.
     42//
     43// The histograms are fitted to a Gaussian, mean and sigma with its errors
     44// and the fit probability are extracted. If none of these values are NaN's and
     45// if the probability is bigger than fProbLimit (default: 0.5%), the fit is valid.
     46// Otherwise, the fit is repeated within ranges of the previous mean +- 5 sigma.
     47// In case this does not make the fit valid, the histogram means and RMS's are
     48// taken directly and the following flags are set:
     49// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeNotFitted ) and
     50// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun    )
     51//
     52// Outliers of more than MHCalibrationChargePix::fPickupLimit (default: 5) sigmas
     53// from the mean are counted as PickUp events (stored in MHCalibrationRelTimePix::fPickup)
     54//
     55// The class also fills arrays with the signal vs. event number, creates a fourier
     56// spectrum and investigates if the projected fourier components follow an exponential
     57// distribution. In case that the probability of the exponential fit is less than
     58// fProbLimit (default: 0.5%), the following flags are set:
     59// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeOscillating ) and
     60// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun      )
     61//
     62// This same procedure is performed for the average pixels.
     63//
     64// The following results are written into MCalibrationRelTimeCam:
     65//
     66// - MCalibrationPix::SetMean()
     67// - MCalibrationPix::SetMeanErr()
     68// - MCalibrationPix::SetSigma()
     69// - MCalibrationPix::SetSigmaErr()
     70// - MCalibrationPix::SetProb()
     71// - MCalibrationPix::SetNumPickup()
     72//
     73// For all averaged areas, the fitted sigma is multiplied with the square root of
     74// the number involved pixels in order to be able to compare it to the average of
     75// sigmas in the camera.
     76//                                                                         
    3177/////////////////////////////////////////////////////////////////////////////
    3278#include "MHCalibrationRelTimeCam.h"
     79#include "MHCalibrationRelTimePix.h"
    3380
    3481#include "MLog.h"
     
    3784#include "MParList.h"
    3885
    39 #include "MHCalibrationRelTimePix.h"
     86#include "MCalibrationRelTimeCam.h"
     87#include "MCalibrationPix.h"
    4088
    4189#include "MArrivalTimeCam.h"
    4290#include "MArrivalTimePix.h"
    4391
     92#include "MGeomCam.h"
     93#include "MGeomPix.h"
     94
     95#include "MBadPixelsCam.h"
     96#include "MBadPixelsPix.h"
     97
    4498ClassImp(MHCalibrationRelTimeCam);
    4599
    46100using namespace std;
    47 const Float_t MHCalibrationRelTimeCam::fgTimeSliceWidth  = 3.3;
    48 const Int_t   MHCalibrationRelTimeCam::fgPulserFrequency = 500;
    49101// --------------------------------------------------------------------------
    50102//
    51 // Default constructor. Creates a MHCalibrationRelTimePix object for each pixel
     103// Default Constructor.
    52104//
    53105MHCalibrationRelTimeCam::MHCalibrationRelTimeCam(const char *name, const char *title)
     
    55107
    56108  fName  = name  ? name  : "MHCalibrationRelTimeCam";
    57   fTitle = title ? title : "Histogram container for the relative time calibration of the camera";
    58  
    59   //
    60   // loop over all Pixels and create two histograms
    61   // one for the Low and one for the High gain
    62   // connect all the histogram with the container fHist
    63   //
    64   fArray = new TObjArray;
    65   fArray->SetOwner();
    66  
    67 
    68   SetTimeSliceWidth();
    69   SetPulserFrequency();
     109  fTitle = title ? title : "Histogram class for the relative time calibration of the camera";
     110 
    70111}
    71112
    72113// --------------------------------------------------------------------------
    73114//
    74 // Delete the array conatining the pixel pedest information
    75 //
    76 MHCalibrationRelTimeCam::~MHCalibrationRelTimeCam()
    77 {
    78   delete fArray;
    79 }
    80 
    81 // --------------------------------------
    82 //
    83 void MHCalibrationRelTimeCam::Clear(Option_t *o)
    84 {
    85 
    86   fArray->ForEach(TObject, Clear)();
    87   return;
    88 }
    89 
    90 // --------------------------------------------------------------------------
    91 //
    92 // Get i-th pixel (pixel number)
    93 //
    94 MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i)
    95 {
    96   return *(MHCalibrationRelTimePix*)(fArray->At(i));
    97 }
    98 
    99 // --------------------------------------------------------------------------
    100 //
    101 // Get i-th pixel (pixel number)
    102 //
    103 const MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i) const
    104 {
    105   return *(MHCalibrationRelTimePix*)(fArray->At(i));
    106 }
    107 
    108 
    109 // --------------------------------------------------------------------------
    110 //
    111 // Our own clone function is necessary since root 3.01/06 or Mars 0.4
    112 // I don't know the reason
    113 //
    114 TObject *MHCalibrationRelTimeCam::Clone(const char *) const
    115 {
    116 
    117   const Int_t n = fArray->GetSize();
    118  
    119   //
    120   // FIXME, this might be done faster and more elegant, by direct copy.
    121   //
    122   MHCalibrationRelTimeCam *cam = new MHCalibrationRelTimeCam;
    123  
    124   cam->fArray->Expand(n);
    125  
    126   for (int i=0; i<n; i++)
    127     {
    128       delete (*cam->fArray)[i];
    129       (*cam->fArray)[i] = (*fArray)[i]->Clone();
    130     }
    131 
    132   return cam;
    133 }
    134 
    135 
    136 
    137 // --------------------------------------------------------------------------
    138 //
    139 //
    140 Bool_t MHCalibrationRelTimeCam::SetupFill(const MParList *pList)
    141 {
     115// Gets or creates the pointers to:
     116// - MCalibrationRelTimeCam
     117//
     118// Searches pointer to:
     119// - MArrivalTimeCam
     120//
     121// Initializes, if empty to MArrivalTimeCam::GetSize() for:
     122// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
     123//
     124// Initializes, if empty to MGeomCam::GetNumAreas() for:
     125// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     126//
     127// Initializes, if empty to MGeomCam::GetNumSectors() for:
     128// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     129//
     130// Calls InitializeHists() for every entry in:
     131// - MHCalibrationCam::fHiGainArray
     132// - MHCalibrationCam::fAverageHiGainAreas
     133// - MHCalibrationCam::fAverageHiGainSectors
     134//
     135// Sets Titles and Names for the Histograms
     136// - MHCalibrationCam::fAverageHiGainAreas
     137// - MHCalibrationCam::fAverageHiGainSectors
     138//
     139// Sets number of bins to MHCalibrationCam::fAverageNbins for:
     140// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     141// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     142//
     143Bool_t MHCalibrationRelTimeCam::ReInitHists(MParList *pList)
     144{
     145
     146  fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationRelTimeCam");
     147  if (!fCam)
     148      return kFALSE;
     149
     150  MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam");
     151  if (!signal)
     152  {
     153      *fLog << err << "MArrivalTimeCam not found... abort." << endl;
     154      return kFALSE;
     155  }
     156
     157  const Int_t npixels  = fGeom->GetNumPixels();
     158  const Int_t nsectors = fGeom->GetNumSectors();
     159  const Int_t nareas   = fGeom->GetNumAreas();
     160
     161  if (fHiGainArray->GetEntries()==0)
     162  {
     163      fHiGainArray->Expand(npixels);
     164      for (Int_t i=0; i<npixels; i++)
     165      {
     166          (*fHiGainArray)[i] = new MHCalibrationRelTimePix;
     167          InitHists((*this)[i],(*fBadPixels)[i],i);
     168      }
     169  }
     170
     171  if (fLoGainArray->GetEntries()==0)
     172  {
     173      fLoGainArray->Expand(npixels);
     174      for (Int_t i=0; i<npixels; i++)
     175      {
     176          (*fLoGainArray)[i] = new MHCalibrationRelTimePix;
     177          InitHists((*this)(i),(*fBadPixels)[i],i);
     178      }
     179  }
     180
     181  if (fAverageHiGainAreas->GetEntries()==0)
     182  {
     183    fAverageHiGainAreas->Expand(nareas);
     184   
     185    for (Int_t j=0; j<nareas; j++)
     186      {
     187        (*fAverageHiGainAreas)[j] =
     188          new MHCalibrationRelTimePix("AverageHiGainArea",
     189                                      "Average Relative Arrival Times area idx ");
     190
     191        InitHists(GetAverageHiGainArea(j),fCam->GetAverageBadArea(j),j);
     192
     193        GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx ");
     194        GetAverageHiGainArea(j).SetNbins(fAverageNbins);
     195      }
     196  }
     197
     198  if (fAverageLoGainAreas->GetEntries()==0)
     199  {
     200    fAverageLoGainAreas->Expand(nareas);
     201   
     202    for (Int_t j=0; j<nareas; j++)
     203      {
     204        (*fAverageLoGainAreas)[j] =
     205          new MHCalibrationRelTimePix("AverageLoGainArea",
     206                                      "Relative Arrival Times average Area idx ");
     207
     208        InitHists(GetAverageLoGainArea(j),fCam->GetAverageBadArea(j),j);
     209
     210        GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx ");
     211        GetAverageLoGainArea(j).SetNbins(fAverageNbins);
     212      }
     213  }
     214
     215  if (fAverageHiGainSectors->GetEntries()==0)
     216  {
     217      fAverageHiGainSectors->Expand(nsectors);
     218
     219      for (Int_t j=0; j<nsectors; j++)
     220      {
     221          (*fAverageHiGainSectors)[j] =
     222            new MHCalibrationRelTimePix("AverageHiGainSector",
     223                                        "Relative Arrival Times average sector ");
     224
     225          GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Relative Arrival Times average Sector ");
     226          GetAverageHiGainSector(j).SetNbins(fAverageNbins);
     227      }
     228  }
     229
     230  if (fAverageLoGainSectors->GetEntries()==0)
     231  {
     232      fAverageLoGainSectors->Expand(nsectors);
     233
     234      for (Int_t j=0; j<nsectors; j++)
     235      {
     236          (*fAverageLoGainSectors)[j] =
     237            new MHCalibrationRelTimePix("AverageLoGainSector",
     238                                        "Relative Arrival Times average sector ");
     239
     240          InitHists(GetAverageLoGainSector(j),fCam->GetAverageBadSector(j),j);
     241         
     242          GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Relative Arrival Times average Sector ");
     243          GetAverageLoGainSector(j).SetNbins(fAverageNbins);
     244      }
     245  }
    142246
    143247  return kTRUE;
    144 
    145 }
    146 
    147 // --------------------------------------------------------------------------
    148 //
    149 Bool_t MHCalibrationRelTimeCam::Fill(const MParContainer *par, const Stat_t w)
     248}
     249
     250
     251//--------------------------------------------------------------------------------
     252//
     253// Retrieves pointer to MArrivalTimeCam:
     254//
     255// Retrieves from MGeomCam:
     256// - number of pixels
     257// - number of pixel areas
     258// - number of sectors
     259//
     260// Fill RelTimes histograms (MHGausEvents::FillHistAndArray()) with:
     261// - MArrivalTimePix::GetArrivalTime(pixid) - MArrivalTimePix::GetArrivalTime(1);
     262//   (i.e. the time difference between pixel i and pixel 1 (hardware number: 2)
     263//
     264Bool_t MHCalibrationRelTimeCam::FillHists(const MParContainer *par, const Stat_t w)
    150265{
    151266
     
    157272    }
    158273 
    159   const Int_t n = arrtime->GetSize();
    160  
    161   if (fArray->GetEntries()==0)
    162     {
    163       fArray->Expand(n);
     274  const Int_t npixels  = fGeom->GetNumPixels();
     275  const Int_t nareas   = fGeom->GetNumAreas();
     276  const Int_t nsectors = fGeom->GetNumSectors();
     277
     278  Float_t sumareahi  [nareas],   sumarealo  [nareas];
     279  Float_t sumsectorhi[nsectors], sumsectorlo[nsectors];
     280  Int_t   numareahi  [nareas],   numarealo  [nareas];
     281  Int_t   numsectorhi[nsectors], numsectorlo[nsectors];
     282
     283  for (UInt_t j=0; j<nareas; j++)
     284    {
     285      sumareahi[j] = sumarealo[j] = 0.;
     286      numareahi[j] = numarealo[j] = 0;
     287    }
     288 
     289  for (UInt_t j=0; j<nsectors; j++)
     290    {
     291      sumsectorhi[j] = sumsectorlo[j] = 0.;
     292      numsectorhi[j] = numsectorlo[j] = 0;
     293    }
     294 
     295  const MArrivalTimePix &refpix = (*arrtime)[1];
     296  const Float_t reftime = refpix.IsLoGainUsed()
     297    ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain();
     298
     299  for (Int_t i=0; i<npixels; i++)
     300    {
     301
     302      MHGausEvents &histhi = (*this)[i];
     303      MHGausEvents &histlo = (*this)(i);
     304
     305      if (histhi.IsExcluded())
     306        continue;
     307
     308      const MArrivalTimePix &pix = (*arrtime)[i];
     309      const Float_t time    = pix.IsLoGainUsed()
     310        ? pix.GetArrivalTimeLoGain() : pix.GetArrivalTimeHiGain(); 
     311      const Float_t reltime = time - reftime;
    164312     
    165       for (Int_t i=0; i<n; i++)
     313      const Int_t aidx   = (*fGeom)[i].GetAidx();
     314      const Int_t sector = (*fGeom)[i].GetSector();
     315
     316      if (pix.IsLoGainUsed())
     317        {
     318          histlo.FillHistAndArray(reltime);
     319          sumarealo  [aidx]   += reltime;
     320          numarealo  [aidx]   ++;
     321          sumsectorlo[sector] += reltime;
     322          numsectorlo[sector] ++;
     323        }
     324      else
    166325        {
    167           (*fArray)[i] = new MHCalibrationRelTimePix;
    168           MHCalibrationRelTimePix &hist = (*this)[i];
    169           hist.ChangeHistId(i);
    170           hist.Init();
    171           hist.SetEventFrequency(fPulserFrequency);
     326          histhi.FillHistAndArray(reltime) ;
     327          sumareahi  [aidx]   += reltime;
     328          numareahi  [aidx]   ++;
     329          sumsectorhi[sector] += reltime;
     330          numsectorhi[sector] ++;
    172331        }
    173332    }
    174333 
    175   if (fArray->GetEntries() != n)
    176     {
    177       gLog << err << "ERROR - Size mismatch... abort." << endl;
    178       return kFALSE;
    179     }
    180  
    181   const MArrivalTimePix &refpix = (*arrtime)[1];
    182   const Float_t reftime = refpix.IsLoGainUsed() ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain();
    183 
    184   for (Int_t i=0; i<n; i++)
    185     {
    186         const MArrivalTimePix &pix = (*arrtime)[i];
    187         const Float_t time = pix.IsLoGainUsed() ? pix.GetArrivalTimeLoGain() : pix.GetArrivalTimeHiGain(); 
    188 
    189         MHCalibrationRelTimePix  &hist = (*this)[i];
    190         hist.FillHistAndArray(time - reftime);
    191     }
    192  
     334  for (UInt_t j=0; j<nareas; j++)
     335    {
     336      MHGausEvents &histhi = GetAverageHiGainArea(j);
     337      histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
     338
     339      MHGausEvents &histlo = GetAverageLoGainArea(j);
     340      histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]);
     341    }
     342 
     343  for (UInt_t j=0; j<nsectors; j++)
     344    {
     345      MHGausEvents &histhi = GetAverageHiGainSector(j);
     346      histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
     347
     348      MHGausEvents &histlo = GetAverageLoGainSector(j);
     349      histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]);
     350    }
     351
    193352  return kTRUE;
    194353}
    195354
    196 Bool_t MHCalibrationRelTimeCam::Finalize()
    197 {
    198 
    199   for (Int_t i=0; i<fArray->GetSize(); i++)
     355// --------------------------------------------------------------------------
     356//
     357// Calls:
     358// - MHCalibrationCam::FitHiGainArrays() with flags:
     359//   MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating
     360// - MHCalibrationCam::FitLoGainArrays() with flags:
     361//   MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating
     362//
     363Bool_t MHCalibrationRelTimeCam::FinalizeHists()
     364{
     365
     366  FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
     367                  MBadPixelsPix::kRelTimeNotFitted,
     368                  MBadPixelsPix::kRelTimeOscillating);
     369
     370  FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
     371                  MBadPixelsPix::kRelTimeNotFitted,
     372                  MBadPixelsPix::kRelTimeOscillating);
     373
     374  return kTRUE;
     375}
     376
     377// --------------------------------------------------------------------------
     378//
     379// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set:
     380// - MBadPixelsPix::kRelTimeNotFitted
     381// - MBadPixelsPix::kRelTimeOscillating
     382//
     383void MHCalibrationRelTimeCam::FinalizeBadPixels()
     384{
     385
     386  for (Int_t i=0; i<fBadPixels->GetSize(); i++)
    200387    {
    201388     
    202       MHCalibrationRelTimePix &hist = (*this)[i];
     389      MBadPixelsPix          &bad    = (*fBadPixels)[i];
     390
     391      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted ))
     392        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
     393
     394      if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating))
     395        bad.SetUnsuitable(   MBadPixelsPix::kUnreliableRun    );
    203396     
    204       //
    205       // 1) Return if the charge distribution is already succesfully fitted
    206       //    or if the histogram is empty
    207       //
    208       if (hist.IsGausFitOK() || hist.IsEmpty())
    209         continue;
    210      
    211       //
    212       // 2) Fit the Hi Gain histograms with a Gaussian
    213       //
    214       hist.FitGaus();
    215      
    216       //
    217       // 3) If fit does not succeed , bypass the fit and take the histogram means and sigmas
    218       //
    219       if (!hist.IsGausFitOK())
    220         hist.BypassFit();
    221      
    222       //
    223       // 4) Create the fourier transform of the arrays
    224       //
    225       hist.CreateFourierSpectrum();
    226      
    227       //
    228       // 5) Renormalize to the real time in ns.
    229       //
    230       hist.Renorm(fTimeSliceWidth);
    231      
    232     }
    233   return kTRUE;
     397    }
    234398}
    235399
     
    250414// =============================================
    251415//
    252 // 4: Returned probability of Gauss fit to Charge distribution
     416// 4: Returned probability of Gauss fit to RelTime distribution
    253417//
    254418// Localized defects:
     
    261425{
    262426
    263   if (fArray->GetSize() <= idx)
     427  if (fHiGainArray->GetSize() <= idx)
    264428    return kFALSE;
    265429
     430  MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx];     
     431
    266432  switch (type)
    267433    {
    268434    case 0:
    269       val = (*this)[idx].GetMean();
     435      val = hist.GetMean();
    270436      break;
    271437    case 1:
    272       val = (*this)[idx].GetMeanErr();
     438      val = hist.GetMeanErr();
    273439      break;
    274440    case 2:
    275       val = (*this)[idx].GetSigma();
     441      val = hist.GetSigma();
    276442      break;
    277443    case 3:
    278       val = (*this)[idx].GetSigmaErr();
     444      val = hist.GetSigmaErr();
    279445      break;
    280446    case 4:
    281       val = (*this)[idx].GetProb();
     447      val = hist.GetProb();
    282448      break;
    283449    case 5:
    284       if (!(*this)[idx].IsGausFitOK())
     450      if (!hist.IsGausFitOK())
    285451        val = 1.;
    286452      break;
    287453    case 6:
    288       if (!(*this)[idx].IsFourierSpectrumOK())
     454      if (!hist.IsFourierSpectrumOK())
    289455        val = 1.;
    290456      break;
     
    297463void MHCalibrationRelTimeCam::DrawPixelContent(Int_t idx) const
    298464{
    299   (*this)[idx].DrawClone();
    300 }
     465  const MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx];     
     466  hist.DrawClone();
     467}
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h

    r3625 r3636  
    22#define MARS_MHCalibrationRelTimeCam
    33
    4 #ifndef ROOT_TObjArray
    5 #include <TObjArray.h>
     4#ifndef MARS_MHCalibrationCam
     5#include "MHCalibrationCam.h"
    66#endif
    77
    8 #ifndef MARS_MH
    9 #include "MH.h"
    10 #endif
    11 #ifndef MARS_MCamEvent
    12 #include "MCamEvent.h"
    13 #endif
    14 
    15 class MHCalibrationRelTimePix;
    16 class MHCalibrationRelTimeCam : public MH, public MCamEvent
     8class MGeomCam;
     9class MHCalibrationRelTimeCam : public MHCalibrationCam
    1710{
    1811
    1912private:
    2013
    21   static const Float_t fgTimeSliceWidth;       // Default for fTimeSliceWidth
    22   static const Int_t   fgPulserFrequency;      // Default for fPulserFrequency
    23 
    24   Float_t fTimeSliceWidth;                    // FADC slice time width
    25   Int_t   fPulserFrequency;                   // The pulser frequency
     14  Bool_t ReInitHists(MParList *pList);
     15  Bool_t FillHists(const MParContainer *par, const Stat_t w=1);
     16  Bool_t FinalizeHists();
     17  void   FinalizeBadPixels();
    2618 
    27   TObjArray  *fArray;       //-> List of MHCalibrationRelTimePix's
    28 
    2919public:
    3020
    3121  MHCalibrationRelTimeCam(const char *name=NULL, const char *title=NULL);
    32   ~MHCalibrationRelTimeCam();
     22  ~MHCalibrationRelTimeCam() {}
    3323
    34   void Clear(Option_t *o="");
    35  
    36   MHCalibrationRelTimePix &operator[](UInt_t i);
    37   const MHCalibrationRelTimePix &operator[](UInt_t i) const;
    38  
    39   Bool_t SetupFill(const MParList *pList);
    40   Bool_t Fill(const MParContainer *par, const Stat_t w=1);
    41   Bool_t Finalize();
    42 
    43   // Setters
    44   void SetTimeSliceWidth( const Float_t width=fgTimeSliceWidth)  {  fTimeSliceWidth = width; }
    45   void SetPulserFrequency( const Int_t  f=fgPulserFrequency)     { fPulserFrequency = f;     }
    46  
    47   TObject *Clone(const char *) const;
    48  
    4924  Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
    5025  void DrawPixelContent(Int_t idx) const;
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc

    r3625 r3636  
    4040using namespace std;
    4141//
    42 const Int_t   MHCalibrationRelTimePix::fgChargeNbins    = 900;
    43 const Axis_t  MHCalibrationRelTimePix::fgChargeFirst    = -13.;
    44 const Axis_t  MHCalibrationRelTimePix::fgChargeLast     =  13.;
     42const Int_t   MHCalibrationRelTimePix::fgRelTimeNbins    = 900;
     43const Axis_t  MHCalibrationRelTimePix::fgRelTimeFirst    = -13.;
     44const Axis_t  MHCalibrationRelTimePix::fgRelTimeLast     =  13.;
    4545// --------------------------------------------------------------------------
    4646//
     
    4848//
    4949// Sets:
    50 // - the default number for fChargeNbins  (fgChargeNbins)
    51 // - the default number for fChargeFirst  (fgChargeFirst)
    52 // - the default number for fChargeLast   (fgChargeLast)
     50// - the default number for fNbins  (fgRelTimeNbins)
     51// - the default number for fFirst  (fgRelTimeFirst)
     52// - the default number for fLast   (fgRelTimeLast)
    5353//
    5454// - the default name of the  fHGausHist ("HCalibrationRelTime")
     
    6464//
    6565MHCalibrationRelTimePix::MHCalibrationRelTimePix(const char *name, const char *title)
    66     : fPixId(-1)
    6766{
    6867
     
    7069  fTitle = title ? title : "Histogrammed Calibration Relative Arrival Time events";
    7170
    72   SetChargeNbins();
    73   SetChargeFirst();
    74   SetChargeLast();
     71  SetNbins ( fgRelTimeNbins );
     72  SetFirst ( fgRelTimeFirst );
     73  SetLast  ( fgRelTimeLast  );
    7574
    7675  // Create a large number of bins, later we will rebin
     
    9695
    9796  fPixId = -1;
     97
    9898  MHGausEvents::Clear();
    9999  return;
     
    109109}
    110110
    111 // --------------------------------------------------------------------------
    112 //
    113 // Sets:
    114 // - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
    115 //
    116 void MHCalibrationRelTimePix::Init()
    117 {
    118111
    119   fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast);
    120 
    121 }
    122 
    123 // --------------------------------------------------------------------------
    124 //
    125 // - Set fPixId to id
    126 //
    127 // Add id to names and titles of:
    128 // - fHGausHist
    129 //
    130 void MHCalibrationRelTimePix::ChangeHistId(Int_t id)
    131 {
    132 
    133   fPixId = id;
    134 
    135   fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
    136   fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
    137 
    138 }
    139 
    140 
    141 // ----------------------------------------------------------------------
    142 //
    143 // Renorm the results from FADC slices to times.
    144 // The parameters slicewidth is the inverse of the FADC frequency and has the unit ns.
    145 //
    146 void MHCalibrationRelTimePix::Renorm(const Float_t slicewidth)
    147 {
    148 
    149   SetMean(     GetMean()    * slicewidth  );
    150   SetMeanErr(  GetMeanErr() * slicewidth  );
    151   SetSigma(    GetSigma()   * slicewidth  );
    152   SetSigmaErr( GetSigmaErr()* slicewidth  );
    153  
    154 }
    155 
  • trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h

    r3625 r3636  
    1111private:
    1212
    13   static const Int_t   fgChargeNbins;      // Default for fChargeNBins  (now set to: 900   )
    14   static const Axis_t  fgChargeFirst;      // Default for fChargeFirst  (now set to: -13.5 )
    15   static const Axis_t  fgChargeLast;       // Default for fChargeLast   (now set to:  13.5 )
    16 
    17   Int_t   fChargeNbins;                   // Number of  bins used for the fHGausHist
    18   Axis_t  fChargeFirst;                   // Lower bound bin used for the fHGausHist
    19   Axis_t  fChargeLast;                    // Upper bound bin used for the fHGausHist
    20 
    21   Int_t fPixId;                           // The pixel ID
     13  static const Int_t   fgRelTimeNbins;      //! Default for MHGausEvents::fNBins  (now set to: 900   )
     14  static const Axis_t  fgRelTimeFirst;      //! Default for MHGausEvents::fFirst  (now set to: -13.5 )
     15  static const Axis_t  fgRelTimeLast;       //! Default for MHGausEvents::fLast   (now set to:  13.5 )
    2216
    2317public:
     
    2620  ~MHCalibrationRelTimePix() {}
    2721 
    28 
    2922  void Clear(Option_t *o="");
    3023  void Reset();
    31   void Init();
    32  
    33   // Setters
    34   void SetChargeNbins(const Int_t  bins =fgChargeNbins)    { fChargeNbins = bins; }
    35   void SetChargeFirst(const Axis_t first=fgChargeFirst)    { fChargeFirst = first; }
    36   void SetChargeLast( const Axis_t last =fgChargeLast)     { fChargeLast  = last; }
    37 
    38   // Others
    39   void ChangeHistId(Int_t i);
    40   void Renorm(const Float_t slicewidth);
    4124 
    4225  ClassDef(MHCalibrationRelTimePix, 1)     // Histogram class for Relative Time Calibration pixel
  • trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc

    r3625 r3636  
    4040//     It is created with a default name and title and resides in directory NULL.
    4141//     - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist
    42 //       (e.g. with the function TH1F::SetBins(..) )
     42//       (e.g. by setting the variables fNbins, fFirst, fLast and calling the function
     43//       InitBins() or by directly calling fHGausHist.SetBins(..) )
    4344//     - The histogram is filled with the functions FillHist() or FillHistAndArray().
    4445//     - The histogram can be fitted with the function FitGaus(). This involves stripping
     
    4849//       The NDF is compared to fNDFLimit and a check is made whether results are NaNs.
    4950//       Accordingly, a flag IsGausFitOK() is set.
     51//     - One can repeat the fit within a given amount of sigmas from the previous mean
     52//       (specified by the variable fPickupLimit) with the function RepeatFit()
     53//     - One can completely skip the fitting to set mean, sigma and its errors directly
     54//       from the histograms with the function BypassFit()
    5055//
    5156//  2) The TArrayF fEvents:
     
    96101const Float_t  MHGausEvents::fgProbLimit            = 0.001;
    97102const Int_t    MHGausEvents::fgNDFLimit             = 2;
     103const Float_t  MHGausEvents::fgPickupLimit          = 5.;
    98104const Int_t    MHGausEvents::fgPowerProbabilityBins = 20;
    99105const Int_t    MHGausEvents::fgBinsAfterStripping   = 40;
     
    105111// - the default Probability Limit for fProbLimit  (fgProbLimit)
    106112// - the default NDF Limit for fNDFLimit  (fgNDFLimit)
     113// - the default number for fPickupLimit  (fgPickupLimit)
    107114// - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping)
    108115// - the default name of the fHGausHist   ("HGausHist")
    109116// - the default title of the fHGausHist  ("Histogram of Events with Gaussian Distribution")
    110117// - the default directory of the fHGausHist (NULL)
     118// - the default for fNbins (100)
     119// - the default for fFirst (0.)
     120// - the default for fLast  (100.)
    111121//
    112122// Initializes:
     
    121131      fPowerSpectrum(NULL),
    122132      fGraphEvents(NULL), fGraphPowerSpectrum(NULL),
     133      fNbins(100), fFirst(0.), fLast(100.), fPixId(-1),
    123134      fHGausHist(), fEvents(0),
    124135      fFGausFit(NULL), fFExpFit(NULL)
     
    133144  SetProbLimit();
    134145  SetNDFLimit();
     146  SetPickupLimit();
    135147  SetBinsAfterStripping();
    136148
     
    183195// --------------------------------------------------------------------------
    184196//
     197// Default Clear(), can be overloaded.
     198//
    185199// Sets:
    186200// - all other pointers to NULL
    187 // - all variables to 0., except fEventFrequency
     201// - all variables to 0., except fPixId to -1 and keep fEventFrequency
    188202// - all flags to kFALSE
    189203//
     
    197211  SetExpFitOK         ( kFALSE );
    198212  SetFourierSpectrumOK( kFALSE );
     213  SetExcluded         ( kFALSE );
    199214
    200215  fMean              = 0.;
     
    205220 
    206221  fCurrentSize       = 0;
     222  fPixId             = -1;
    207223
    208224  if (fHPowerProbability)
     
    247263
    248264// --------------------------------------------------------------------------
     265//
     266// Default Reset(), can be overloaded.
    249267//
    250268// Executes:
     
    264282// --------------------------------------------------------------------------
    265283//
     284// Default InitBins, can be overloaded.
     285//
     286// Executes:
     287// - fHGausHist.SetBins(fNbins,fFirst,fLast)
     288//
     289void MHGausEvents::InitBins()
     290{
     291  fHGausHist.SetBins(fNbins,fFirst,fLast);
     292}
     293
     294// --------------------------------------------------------------------------
     295//
    266296// Executes:
    267297// - FillArray()
     
    303333}
    304334
     335// --------------------------------------------------------------------------
     336//
     337// Set Excluded bit from outside
     338//
     339void MHGausEvents::SetExcluded(const Bool_t b)
     340{
     341    b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded);
     342}
     343
    305344
    306345const Double_t MHGausEvents::GetChiSquare()  const
     
    366405{
    367406  return TESTBIT(fFlags,kExpFitOK);
     407}
     408
     409const Bool_t MHGausEvents::IsExcluded() const
     410{
     411  return TESTBIT(fFlags,kExcluded);
    368412}
    369413
     
    409453  if (fFExpFit)
    410454    return;
     455
     456  if (fEvents.GetSize() < 8)
     457    {
     458      *fLog << warn << "Cannot create Fourier spectrum in pixel: " << fPixId
     459            << ". Number of events smaller than 8 " << endl;
     460      return;
     461    }
     462 
    411463
    412464  //
     
    424476  // This cuts only the non-used zero's, but MFFT will later cut the rest
    425477  MArray::StripZeros(fEvents);
     478
     479  if (fEvents.GetSize() < 8)
     480    {
     481      *fLog << warn << "Cannot create Fourier spectrum. " << endl;
     482      *fLog << warn << "Number of events (after stripping of zeros) is smaller than 8 "
     483            << "in pixel: " << fPixId << endl;
     484      return;
     485    }
    426486
    427487  MFFT fourier;
     
    506566  if (!fFGausFit)
    507567    {
    508     *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl;
     568    *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit "
     569          << "in pixel: " << fPixId << endl;
    509570    return kFALSE;
    510571    }
     
    528589  // The fit result is accepted under condition:
    529590  // 1) The results are not nan's
    530   // 2) The NDF is not smaller than fNDFLimit (5)
    531   // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%)
     591  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
     592  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
    532593  //
    533594  if (   TMath::IsNaN(fMean)
     
    545606
    546607// -----------------------------------------------------------------------------
     608//
     609// If flag IsGausFitOK() is set (histogram already successfully fitted),
     610// returns kTRUE
     611//
     612// If both fMean and fSigma are still zero, call FitGaus()
     613//
     614// Repeats the Gauss fit in a smaller range, defined by:
     615//
     616// min = GetMean() - fPickupLimit * GetSigma();
     617// max = GetMean() + fPickupLimit * GetSigma();
     618//
     619// The fit results are retrieved and stored in class-own variables. 
     620//
     621// A flag IsGausFitOK() is set according to whether the fit probability
     622// is smaller or bigger than fProbLimit, whether the NDF is bigger than
     623// fNDFLimit and whether results are NaNs.
     624//
     625Bool_t MHGausEvents::RepeatFit(const Option_t *option)
     626{
     627
     628  if (IsGausFitOK())
     629    return kTRUE;
     630
     631  if ((fMean == 0.) && (fSigma == 0.))
     632    return FitGaus();
     633
     634  //
     635  // Get new fitting ranges
     636  //
     637  Axis_t rmin = fMean - fPickupLimit * fSigma;
     638  Axis_t rmax = fMean + fPickupLimit * fSigma;
     639
     640  Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst());
     641  Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ;
     642
     643  fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax);
     644
     645  fHGausHist.Fit(fFGausFit,option);
     646
     647  fMean     = fFGausFit->GetParameter(1);
     648  fMeanErr  = fFGausFit->GetParameter(2);
     649  fSigma    = fFGausFit->GetParError(1) ;
     650  fSigmaErr = fFGausFit->GetParError(2) ;
     651  fProb     = fFGausFit->GetProb()      ;     
     652
     653  //
     654  // The fit result is accepted under condition:
     655  // 1) The results are not nan's
     656  // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit)
     657  // 3) The Probability is greater than fProbLimit (default: fgProbLimit)
     658  //
     659  if (   TMath::IsNaN ( fMean     )
     660      || TMath::IsNaN ( fMeanErr  )
     661      || TMath::IsNaN ( fProb     )   
     662      || TMath::IsNaN ( fSigma    )
     663      || TMath::IsNaN ( fSigmaErr )
     664      || fFGausFit->GetNDF() < fNDFLimit
     665      || fProb < fProbLimit )
     666    return kFALSE;
     667 
     668  SetGausFitOK(kTRUE);
     669  return kTRUE;
     670
     671}
     672
     673// -----------------------------------------------------------------------------
    547674//
    548675// Bypasses the Gauss fit by taking mean and RMS from the histogram
     
    559686  if (entries <= 0.)
    560687    {
    561       *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0" << endl;
     688      *fLog << warn << GetDescriptor()
     689            << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl;
    562690      return;
    563691    }
    564692 
    565   SetMean     ( fHGausHist.GetMean() );
    566   SetMeanErr  ( fHGausHist.GetRMS() / TMath::Sqrt(entries) );
    567   SetSigma    ( fHGausHist.GetRMS() );
    568   SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2. );
     693  fMean     = fHGausHist.GetMean();
     694  fMeanErr  = fHGausHist.GetRMS() / TMath::Sqrt(entries);
     695  fSigma    = fHGausHist.GetRMS() ;
     696  fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.;
    569697}
    570698
     
    577705 
    578706  *fLog << all                                                        << endl;
    579   *fLog << all << "Results of the Gauss Fit: "                        << endl;
     707  *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId     << endl;
    580708  *fLog << all << "Mean: "             << GetMean()                   << endl;
    581709  *fLog << all << "Sigma: "            << GetSigma()                  << endl;
     
    585713  *fLog << all                                                        << endl;
    586714 
     715}
     716
     717// --------------------------------------------------------------------------
     718//
     719// - Set fPixId to id
     720//
     721// Add id to names and titles of:
     722// - fHGausHist
     723//
     724void MHGausEvents::ChangeHistId(Int_t id)
     725{
     726
     727  fPixId = id;
     728
     729  fHGausHist.SetName(  Form("%s%d", fHGausHist.GetName(),  id));
     730  fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id));
     731
     732  fName  = Form("%s%d", fName.Data(),  id);
     733  fTitle = Form("%s%d", fTitle.Data(), id);
     734
    587735}
    588736
     
    766914}
    767915
    768 
     916const Double_t MHGausEvents::GetPickup() const
     917{
     918 
     919  if ((fMean == 0.) && (fSigma == 0.))
     920    return -1.;
     921
     922  const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma);
     923  const Int_t last  = fHGausHist.GetXaxis()->GetLast();
     924
     925  if (first >= last)
     926    return 0.;
     927 
     928  return fHGausHist.Integral(first, last, "width");
     929}
     930
     931
  • trunk/MagicSoft/Mars/mcalib/MHGausEvents.h

    r3625 r3636  
    2020private:
    2121
    22   const static Float_t  fgProbLimit;            // Default for fProbLimit (now 0.001)
    23   const static Int_t    fgNDFLimit;             // Default for fNDFLimit  (now 2)
    24   const static Int_t    fgPowerProbabilityBins; // Default for fPowerProbabilityBins (now 20)
    25   const static Int_t    fgBinsAfterStripping;   // Default for fBinsAfterStripping (now 40)
    26 
     22  const static Float_t  fgProbLimit;            //! Default for fProbLimit (now set to: 0.001)
     23  const static Int_t    fgNDFLimit;             //! Default for fNDFLimit  (now set to: 2)
     24  const static Float_t  fgPickupLimit;          //! Default for fPickupLimit (now set to: 5. )
     25  const static Int_t    fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20)
     26  const static Int_t    fgBinsAfterStripping;   //! Default for fBinsAfterStripping (now set to: 40)
     27 
    2728  Int_t    fPowerProbabilityBins;      // Bins for the projected power spectrum
    2829  Int_t    fBinsAfterStripping;        // Bins for the Gauss Histogram after stripping off the zeros at both ends
     
    4142  Double_t fProb;                      //  Probability of the Gauss fit (derived from Chi-Square and NDF
    4243
    43   enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK }; // Bits to hold information about fit results
     44  enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK,
     45         kExcluded };                  //  Bits to hold information about fit results
    4446 
    45   Byte_t fFlags;                       //   Byte to hold the bits fit result bits
     47  Byte_t fFlags;                       //  Byte to hold the bits fit result bits
    4648 
    47   Int_t fCurrentSize;                  //   Current size of the array fEvents
    48  
     49  Int_t fCurrentSize;                  //  Current size of the array fEvents
     50
    4951protected:
    5052
     53  Int_t   fNbins;                      // Number histogram bins for fHGausHist (used by Init())
     54  Axis_t  fFirst;                      // Lower histogram edge  for fHGausHist (used by Init())
     55  Axis_t  fLast;                       // Upper histogram edge  for fHGausHist (used by Init())
     56  Int_t   fPixId;                      //  Pixel ID
     57 
    5158  TH1F    fHGausHist;                  // Histogram which should hold the Gaussian distribution
    5259  TArrayF fEvents;                     // Array which holds the entries of GausHist
     
    5764  Float_t  fProbLimit;                 // Probability limit for judgement if fit is OK
    5865  Int_t    fNDFLimit;                  // NDF limit for judgement if fit is OK
     66  Float_t  fPickupLimit;               // Upper number of sigmas from the mean until events are considered as pickup
    5967
    6068  Float_t *CreateEventXaxis(Int_t n);  //   Create an x-axis for the Event TGraphs
     
    7583  virtual void Clear(Option_t *o="");
    7684  virtual void Reset(); 
    77 
     85  virtual void InitBins();
     86 
    7887  // Setters
    7988  void  SetEventFrequency(const Float_t f)   { fEventFrequency = f; }
     
    8493  void  SetSigmaErr( const Double_t d )   { fSigmaErr = d;   }
    8594  void  SetProb    ( const Double_t d )   { fProb     = d;   }
     95  void  SetPixId    ( const Int_t i    )   { fPixId    = i;   }
     96
     97  void  SetNbins    ( const Int_t i    )   { fNbins    = i;   } 
     98  void  SetFirst   ( const Double_t d )   { fFirst    = d;   }
     99  void  SetLast    ( const Double_t d )   { fLast     = d;   }
    86100 
     101  void  SetNDFLimit(  const Int_t   lim=fgNDFLimit  ) {  fNDFLimit = lim;  } 
     102  void  SetPickupLimit( const Float_t  lim =fgPickupLimit)  { fPickupLimit  = lim;   }
    87103  void  SetProbLimit( const Float_t lim=fgProbLimit ) {  fProbLimit = lim; }
    88   void  SetNDFLimit(  const Int_t   lim=fgNDFLimit  ) {  fNDFLimit = lim;  } 
    89104
    90105  // Setters ONLY for MC:
    91   void  SetGausFitOK(         const Bool_t b );
    92   void  SetExpFitOK(          const Bool_t b );
    93   void  SetFourierSpectrumOK( const Bool_t b );
     106  void  SetGausFitOK(         const Bool_t b=kTRUE );
     107  void  SetExpFitOK(          const Bool_t b=kTRUE );
     108  void  SetFourierSpectrumOK( const Bool_t b=kTRUE );
     109  void  SetExcluded         ( const Bool_t b=kTRUE ); 
    94110
    95111  // Getters
     
    101117  const Double_t GetProb()          const { return fProb;      }
    102118  const Int_t    GetNdf()           const;
     119  const Double_t GetPickup()        const;
     120  const Int_t    GetPixId()         const { return fPixId;     }
    103121
    104122  const Double_t GetSlope()        const;
     
    136154  const Bool_t IsEmpty()              const;
    137155  const Bool_t IsFourierSpectrumOK()  const;
     156  const Bool_t IsExcluded         ()  const;
    138157
    139158  // Fill
     
    145164  Bool_t FitGaus(Option_t *option="RQ0",
    146165                 const Double_t xmin=0., const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian
    147   void BypassFit();   // Take mean and RMS from the histogram
    148  
     166  Bool_t RepeatFit(const Option_t *option="RQ0");      // Repeat fit within limits defined by fPickupLimit
     167  void BypassFit();                                    // Take mean and RMS from the histogram
    149168 
    150169  // Draws
     
    155174 
    156175  // Miscelleaneous
     176  virtual void ChangeHistId(Int_t id);   // Changes names and titles of the histogram
    157177  void CreateFourierSpectrum();             // Create the fourier spectrum out of fEvents
    158178  void CreateGraphEvents();                 // Create the TGraph fGraphEvents of fEvents
Note: See TracChangeset for help on using the changeset viewer.