Ignore:
Timestamp:
04/04/04 17:56:14 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars/manalysis
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/manalysis/MHPedestalCam.cc

    r3446 r3643  
    2828// MHPedestalCam                                                            //
    2929//                                                                         //
    30 // Hold the Pedestal information for all pixels in the camera              //
    31 //                                                                         //
     30// Fills the extracted pedestals of MExtractedSignalCam into
     31// the MHGausEvents-classes MHPedestalPix for every:
     32//
     33// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray 
     34//   or MHCalibrationCam::fHiGainArray, respectively, depending if
     35//   MArrivalTimePix::IsLoGainUsed() is set.
     36//
     37// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
     38//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
     39//   MHCalibrationCam::fAverageHiGainAreas
     40//
     41// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
     42//   stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
     43//   and MHCalibrationCam::fAverageHiGainSectors
     44//
     45// Every pedestal is directly taken from MExtractedSignalCam, filled by the
     46// appropriate extractor.
     47//
     48// The pedestals are filled into a histogram and an array, in order to perform
     49// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
     50// event-by-event basis and written into the corresponding average pixels.
     51//
     52// The histograms are fitted to a Gaussian, mean and sigma with its errors
     53// and the fit probability are extracted. If none of these values are NaN's and
     54// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
     55// the fit is declared valid.
     56// Otherwise, the fit is repeated within ranges of the previous mean
     57// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit())
     58// In case this does not make the fit valid, the histogram means and RMS's are
     59// taken directly (see MHGausEvents::BypassFit()).
     60//
     61// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas
     62// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup)
     63//
     64// The number of summed FADC slices is taken to re-normalize
     65// the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()):
     66// - Mean Pedestal        / slice = Mean Pedestal        / Number slices
     67// - Mean Pedestal Error  / slice = Mean Pedestal Error  / Number slices
     68// - Sigma Pedestal       / slice = Sigma Pedestal       / Sqrt (Number slices)
     69// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices)
     70//
     71// The class also fills arrays with the signal vs. event number, creates a fourier
     72// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
     73// projected fourier components follow an exponential distribution.
     74//
     75// This same procedure is performed for the average pixels.
     76//
     77// The following results are written into MPedestalCam:
     78//
     79// - MCalibrationPix::SetMean()
     80// - MCalibrationPix::SetMeanErr()
     81// - MCalibrationPix::SetSigma()
     82// - MCalibrationPix::SetSigmaErr()
     83// - MCalibrationPix::SetProb()
     84// - MCalibrationPix::SetNumPickup()
     85//
     86// For all averaged areas, the fitted sigma is multiplied with the square root of
     87// the number involved pixels in order to be able to compare it to the average of
     88// sigmas in the camera.
     89//
    3290/////////////////////////////////////////////////////////////////////////////
    3391#include "MHPedestalCam.h"
     
    46104#include "MPedestalPix.h"
    47105
     106#include "MGeomCam.h"
     107#include "MGeomPix.h"
     108
    48109ClassImp(MHPedestalCam);
    49110
     
    51112// --------------------------------------------------------------------------
    52113//
    53 // Default constructor. Creates a MHPedestalPix object for each pixel
     114// Default constructor.
     115//
     116// Initializes:
     117// - fExtractHiGainSlices to 1.
     118// - fExtractLoGainSlices to 1.
    54119//
    55120MHPedestalCam::MHPedestalCam(const char *name, const char *title)
    56     : fExtractSlices(0.)
     121    : fExtractHiGainSlices(1.), fExtractLoGainSlices(1.)
    57122{
    58123    fName  = name  ? name  : "MHPedestalCam";
    59124    fTitle = title ? title : "";
    60125
    61     //
    62     // loop over all Pixels and create two histograms
    63     // one for the Low and one for the High gain
    64     // connect all the histogram with the container fHist
    65     //
    66     fArray = new TObjArray;
    67     fArray->SetOwner();
    68 
    69 }
     126}
     127
     128
    70129
    71130// --------------------------------------------------------------------------
    72131//
    73 // Delete the array conatining the pixel pedest information
    74 //
    75 MHPedestalCam::~MHPedestalCam()
    76 {
    77   delete fArray;
    78 }
    79 
    80 // --------------------------------------------------------------------------
    81 //
    82 // Get i-th pixel (pixel number)
    83 //
    84 MHPedestalPix &MHPedestalCam::operator[](UInt_t i)
    85 {
    86   return *(MHPedestalPix*)(fArray->At(i));
    87 }
    88 
    89 // --------------------------------------------------------------------------
    90 //
    91 // Get i-th pixel (pixel number)
    92 //
    93 const MHPedestalPix &MHPedestalCam::operator[](UInt_t i) const
    94 {
    95   return *(MHPedestalPix*)(fArray->At(i));
    96 }
    97 
    98 
    99 // --------------------------------------------------------------------------
    100 //
    101 // Our own clone function is necessary since root 3.01/06 or Mars 0.4
    102 // I don't know the reason
    103 //
    104 TObject *MHPedestalCam::Clone(const char *) const
    105 {
    106 
    107   const Int_t n = fArray->GetSize();
    108  
    109   //
    110   // FIXME, this might be done faster and more elegant, by direct copy.
    111   //
    112   MHPedestalCam *cam = new MHPedestalCam;
    113  
    114   cam->fArray->Expand(n);
    115  
    116   for (int i=0; i<n; i++)
    117     {
    118       delete (*cam->fArray)[i];
    119       (*cam->fArray)[i] = (*fArray)[i]->Clone();
    120     }
    121 
    122   return cam;
    123 }
    124 
    125 
    126 
    127 // --------------------------------------------------------------------------
    128 //
    129 // To setup the object we get the number of pixels from a MGeomCam object
    130 // in the Parameter list.
    131 // MHPedestalPix sets its parameters to 0. (other than default which is -1.)
    132 //
    133 Bool_t MHPedestalCam::SetupFill(const MParList *pList)
    134 {
     132// Searches pointer to:
     133// - MPedestalCam
     134// - MExtractedSignalCam
     135//
     136// Retrieves from MExtractedSignalCam:
     137// - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
     138// - number of used Low  Gain FADC slices (MExtractedSignalCam::GetNumUsedLoGainFADCSlices())
     139//
     140// Initializes, if empty to MGeomCam::GetNumPixels():
     141// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
     142//
     143// Initializes, if empty to MGeomCam::GetNumAreas() for:
     144// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     145//
     146// Initializes, if empty to MGeomCam::GetNumSectors() for:
     147// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     148//
     149// Calls MHCalibrationCam::InitPedHists() for every entry in:
     150// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
     151// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     152// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     153//
     154// Sets Titles and Names for the Histograms
     155// - MHCalibrationCam::fAverageHiGainAreas
     156// - MHCalibrationCam::fAverageHiGainSectors
     157//
     158// Sets number of bins to MHCalibrationCam::fAverageNbins for:
     159// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
     160// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
     161//
     162Bool_t MHPedestalCam::ReInitHists(MParList *pList)
     163{
     164
     165  MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
     166  if (!signal)
     167    {
     168      gLog << err << "Cannot find MExtractedSignalCam... abort." << endl;
     169      return kFALSE;
     170    }
     171 
    135172
    136173  fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
     
    142179    }
    143180
    144   fArray->Delete();
     181  Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
     182  Float_t sliceslo = signal->GetNumUsedLoGainFADCSlices();
     183
     184  if (sliceshi == 0.)
     185    {
     186      gLog << err << "Number of used signal slices in MExtractedSignalCam is zero  ... abort."
     187           << endl;
     188      return kFALSE;
     189    }
     190
     191  if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
     192    {
     193      gLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam  ... abort."
     194           << endl;
     195      return kFALSE;
     196    }
     197
     198  if (fExtractLoGainSlices != 0. && sliceslo != fExtractLoGainSlices )
     199    {
     200      gLog << err << "Number of used Low Gain signal slices changed in MExtractedSignalCam  ... abort."
     201           << endl;
     202      return kFALSE;
     203    }
     204
     205  fExtractHiGainSlices = sliceshi;
     206  fExtractLoGainSlices = sliceslo;
     207
     208  const Int_t npixels  = fGeom->GetNumPixels();
     209  const Int_t nsectors = fGeom->GetNumSectors();
     210  const Int_t nareas   = fGeom->GetNumAreas();
     211
     212  if (fHiGainArray->GetEntries()==0)
     213  {
     214      fHiGainArray->Expand(npixels);
     215      for (Int_t i=0; i<npixels; i++)
     216      {
     217          (*fHiGainArray)[i] = new MHPedestalPix;
     218          InitPedHists((MHPedestalPix&)(*this)[i],i,fExtractHiGainSlices);
     219      }
     220  }
     221
     222  if (fLoGainArray->GetEntries()==0)
     223  {
     224      fLoGainArray->Expand(npixels);
     225      for (Int_t i=0; i<npixels; i++)
     226      {
     227          (*fLoGainArray)[i] = new MHPedestalPix;
     228          InitPedHists((MHPedestalPix&)(*this)(i),i,fExtractLoGainSlices);
     229      }
     230  }
     231
     232  if (fAverageHiGainAreas->GetEntries()==0)
     233  {
     234    fAverageHiGainAreas->Expand(nareas);
     235   
     236    for (Int_t j=0; j<nareas; j++)
     237      {
     238        (*fAverageHiGainAreas)[j] =
     239          new MHPedestalPix("AverageHiGainArea",
     240                                      "Average Pedestals area idx ");
     241
     242        InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices);
     243
     244        GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
     245        GetAverageHiGainArea(j).SetNbins(fAverageNbins);
     246      }
     247  }
     248
     249  if (fAverageLoGainAreas->GetEntries()==0)
     250  {
     251    fAverageLoGainAreas->Expand(nareas);
     252   
     253    for (Int_t j=0; j<nareas; j++)
     254      {
     255        (*fAverageLoGainAreas)[j] =
     256          new MHPedestalPix("AverageLoGainArea",
     257                                      "Pedestals average Area idx ");
     258
     259        InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices);
     260
     261        GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
     262        GetAverageLoGainArea(j).SetNbins(fAverageNbins);
     263      }
     264  }
     265
     266  if (fAverageHiGainSectors->GetEntries()==0)
     267  {
     268      fAverageHiGainSectors->Expand(nsectors);
     269
     270      for (Int_t j=0; j<nsectors; j++)
     271      {
     272          (*fAverageHiGainSectors)[j] =
     273            new MHPedestalPix("AverageHiGainSector",
     274                                        "Pedestals average sector ");
     275
     276          InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices);
     277
     278          GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
     279          GetAverageHiGainSector(j).SetNbins(fAverageNbins);
     280      }
     281  }
     282
     283  if (fAverageLoGainSectors->GetEntries()==0)
     284  {
     285      fAverageLoGainSectors->Expand(nsectors);
     286
     287      for (Int_t j=0; j<nsectors; j++)
     288      {
     289          (*fAverageLoGainSectors)[j] =
     290            new MHPedestalPix("AverageLoGainSector",
     291                                        "Pedestals average sector ");
     292
     293          InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices);
     294         
     295          GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
     296          GetAverageLoGainSector(j).SetNbins(fAverageNbins);
     297      }
     298  }
     299
    145300  return kTRUE;
    146301
    147302}
    148303
    149 // --------------------------------------------------------------------------
    150 //
    151 Bool_t MHPedestalCam::Fill(const MParContainer *par, const Stat_t w)
     304
     305// -------------------------------------------------------------
     306//
     307// If MBadPixelsPix::IsBad():
     308// - calls MHGausEvents::SetExcluded()
     309//
     310// Calls:
     311// - MHGausEvents::InitBins()
     312// - MHGausEvents::ChangeHistId(i)
     313// - MHGausEvents::SetEventFrequency(fPulserFrequency)
     314// - MHPedestalPix::SetNSlices(nslices)
     315//
     316void MHPedestalCam::InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices)
     317{
     318 
     319  hist.InitBins();
     320  hist.ChangeHistId(i);
     321  hist.SetEventFrequency(fPulserFrequency);
     322  hist.SetNSlices(nslices);
     323}
     324
     325
     326
     327// -------------------------------------------------------------------------------
     328//
     329// Retrieves pointer to MExtractedSignalCam:
     330//
     331// Retrieves from MGeomCam:
     332// - number of pixels
     333// - number of pixel areas
     334// - number of sectors
     335//
     336// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively
     337// with the signals MExtractedSignalCam::GetExtractedSignalHiGain() and
     338// MExtractedSignalCam::GetExtractedSignalLoGain(), respectively.
     339//
     340Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
    152341{
    153342
     
    159348    }
    160349 
    161 
    162   Float_t slices = signal->GetNumUsedHiGainFADCSlices();
    163 
    164   if (slices == 0.)
    165     {
    166       gLog << err << "Number of used signal slices in MExtractedSignalCam is zero  ... abort."
    167            << endl;
    168       return kFALSE;
    169     }
    170 
    171   if (fExtractSlices != 0. && slices != fExtractSlices )
    172     {
    173       gLog << err << "Number of used signal slices changed in MExtractedSignalCam  ... abort."
    174            << endl;
    175       return kFALSE;
    176     }
    177 
    178   fExtractSlices = slices;
    179 
    180   const Int_t n = signal->GetSize();
    181  
    182   if (fArray->GetEntries()==0)
    183     {
    184       fArray->Expand(n);
    185      
    186       for (Int_t i=0; i<n; i++)
    187         {
    188           (*fArray)[i] = new MHPedestalPix;
    189           MPedestalPix  &pix  = (*fPedestals)[i];
    190           pix.InitUseHists();
    191           MHPedestalPix &hist = (*this)[i];
    192           hist.ChangeHistId(i);
    193           hist.InitBins();
    194         }
    195     }
    196  
    197   if (fArray->GetEntries() != n)
    198     {
    199       gLog << err << "ERROR - Size mismatch... abort." << endl;
    200       return kFALSE;
    201     }
    202  
    203   for (Int_t i=0; i<n; i++)
    204     {
     350  const Int_t npixels  = fGeom->GetNumPixels();
     351  const Int_t nareas   = fGeom->GetNumAreas();
     352  const Int_t nsectors = fGeom->GetNumSectors();
     353
     354  Float_t sumareahi  [nareas],   sumarealo  [nareas];
     355  Float_t sumsectorhi[nsectors], sumsectorlo[nsectors];
     356  Int_t   numareahi  [nareas],   numarealo  [nareas];
     357  Int_t   numsectorhi[nsectors], numsectorlo[nsectors];
     358
     359  for (UInt_t j=0; j<nareas; j++)
     360    {
     361      sumareahi[j] = sumarealo[j] = 0.;
     362      numareahi[j] = numarealo[j] = 0;
     363    }
     364 
     365  for (UInt_t j=0; j<nsectors; j++)
     366    {
     367      sumsectorhi[j] = sumsectorlo[j] = 0.;
     368      numsectorhi[j] = numsectorlo[j] = 0;
     369    }
     370 
     371
     372  for (Int_t i=0; i<npixels; i++)
     373    {
     374
     375      MHGausEvents &histhi = (*this)[i];
     376      MHGausEvents &histlo = (*this)(i);
     377
     378      if (histhi.IsExcluded())
     379        continue;
    205380
    206381      const MExtractedSignalPix &pix = (*signal)[i];
    207382     
    208       const Float_t signal = pix.GetExtractedSignalHiGain();
    209 
    210       MHPedestalPix  &hist = (*this)[i];
    211       //
    212       // Don't fill signal per slice, we get completely screwed up
    213       // with the sigma. Better fill like it is and renorm later
    214       //
    215       hist.FillHistAndArray(signal);
     383      const Float_t pedhi = pix.GetExtractedSignalHiGain();
     384      const Float_t pedlo = pix.GetExtractedSignalLoGain();
     385
     386      const Int_t aidx   = (*fGeom)[i].GetAidx();
     387      const Int_t sector = (*fGeom)[i].GetSector();
     388
     389      histhi.FillHistAndArray(pedhi) ;
     390      sumareahi  [aidx]   += pedhi;
     391      numareahi  [aidx]   ++;
     392      sumsectorhi[sector] += pedhi;
     393      numsectorhi[sector] ++;
     394
     395      histlo.FillHistAndArray(pedlo);
     396      sumarealo  [aidx]   += pedlo;
     397      numarealo  [aidx]   ++;
     398      sumsectorlo[sector] += pedlo;
     399      numsectorlo[sector] ++;
     400
     401    }
     402 
     403  for (UInt_t j=0; j<nareas; j++)
     404    {
     405      MHGausEvents &histhi = GetAverageHiGainArea(j);
     406      histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
     407
     408      MHGausEvents &histlo = GetAverageLoGainArea(j);
     409      histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]);
     410    }
     411 
     412  for (UInt_t j=0; j<nsectors; j++)
     413    {
     414      MHGausEvents &histhi = GetAverageHiGainSector(j);
     415      histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
     416
     417      MHGausEvents &histlo = GetAverageLoGainSector(j);
     418      histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]);
    216419    }
    217420 
     
    219422}
    220423
    221 Bool_t MHPedestalCam::Finalize()
    222 {
    223     for (Int_t i=0; i<fArray->GetSize(); i++)
    224     {
    225 
    226         MHPedestalPix &hist = (*this)[i];
    227 
    228         //
    229         // 1) Return if the charge distribution is already succesfully fitted
    230         //    or if the histogram is empty
    231         //
    232         if (hist.IsGausFitOK() || hist.IsEmpty())
    233           continue;
    234 
    235         //
    236         // 2) Fit the Hi Gain histograms with a Gaussian
    237         //
    238         hist.FitGaus();
    239         hist.Renorm(fExtractSlices);
    240         hist.CreateFourierSpectrum();
    241 
    242     }
    243     return kTRUE;
    244 }
    245 
    246424// --------------------------------------------------------------------------
     425//
     426// Calls:
     427// - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
     428// - MHCalibrationCam::FitLoGainArrays() with Bad Pixels flags 0
     429//
     430Bool_t MHPedestalCam::FinalizeHists()
     431{
     432
     433  FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
     434                  MBadPixelsPix::kHiGainNotFitted,
     435                  MBadPixelsPix::kHiGainOscillating);
     436  FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels,
     437                  MBadPixelsPix::kLoGainNotFitted,
     438                  MBadPixelsPix::kLoGainOscillating);
     439
     440  return kTRUE;
     441 
     442}
     443
     444// ------------------------------------------------------------------
    247445//
    248446// The types are as follows:
     
    277475{
    278476
    279   if (fArray->GetSize() <= idx)
     477  if (fHiGainArray->GetSize() <= idx)
    280478    return kFALSE;
    281479
     
    344542}
    345543
     544// --------------------------------------------------------------------------
     545//
     546// Calls MHGausEvents::DrawClone() for pixel idx
     547//
    346548void MHPedestalCam::DrawPixelContent(Int_t idx) const
    347549{
  • trunk/MagicSoft/Mars/manalysis/MHPedestalCam.h

    r3642 r3643  
    3232  void DrawPixelContent(Int_t idx) const;
    3333 
    34   ClassDef(MHPedestalCam, 1)    // Storage Container for all pedestal information of the camera
     34  ClassDef(MHPedestalCam, 1)    // Histogram class for Charge Camera Pedestals
    3535};
    3636
  • trunk/MagicSoft/Mars/manalysis/MHPedestalPix.h

    r3642 r3643  
    1111private:
    1212
    13   static const Int_t   fgChargeNbins;        // Default for fNBins          (now set to: 450  )
     13  static const Int_t   fgChargeNbins;        // Default for fNbins          (now set to: 450  )
    1414  static const Axis_t  fgChargeFirst;        // Default for fFirst          (now set to: -0.5  )
    1515  static const Axis_t  fgChargeLast;         // Default for fLast           (now set to: 449.5)
     
    3131  void Renorm(); 
    3232
    33   ClassDef(MHPedestalPix, 1)     // Histograms for each calibrated pixel
     33  ClassDef(MHPedestalPix, 1)     // Histogram class for Charge Pedestal Pixel
    3434};
    3535
Note: See TracChangeset for help on using the changeset viewer.