Ignore:
Timestamp:
03/10/04 16:50:25 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3446 r3455  
    138138#include "MCalibrationChargePINDiode.h"
    139139
    140 
    141140ClassImp(MCalibrationChargeCam);
    142141
    143142using namespace std;
    144143
     144const Float_t MCalibrationChargeCam::gkAverageQE                = 0.18;     
     145const Float_t MCalibrationChargeCam::gkAverageQEErr             = 0.02; 
     146const Float_t MCalibrationChargeCam::fgConvFFactorRelErrLimit   = 0.25;
    145147// --------------------------------------------------------------------------
    146148//
    … …  
    171173
    172174    Clear();
     175
     176    SetAverageQE();
     177    SetConvFFactorRelErrLimit();
    173178}
    174179
    … …  
    261266  fAverageOuterBadPix->Clear();
    262267
    263   fNumExcludedPixels                 = 0;
     268  fNumExcludedPixels            = 0;
     269  fMeanFluxPhesInnerPixel       = 0.;
     270  fMeanFluxPhesInnerPixelErr    = 0.;
     271  fMeanFluxPhesOuterPixel       = 0.;
     272  fMeanFluxPhesOuterPixelErr    = 0.;
    264273
    265274  CLRBIT(fFlags,kBlindPixelMethodValid);
     275  CLRBIT(fFlags,kFFactorMethodValid);
    266276  CLRBIT(fFlags,kPINDiodeMethodValid);
    267277
    268278  return;
     279}
     280
     281void MCalibrationChargeCam::SetFFactorMethodValid(const Bool_t b)
     282{
     283    b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
    269284}
    270285
    … …  
    567582      break;
    568583    case 9:
    569       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     584      if ((*this)[idx].IsExcluded()
     585          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     586          || !(*this)[idx].IsFFactorMethodValid())
    570587        return kFALSE;
    571588      val = (*this)[idx].GetPheFFactorMethod();
    572589      break;
    573590    case 10:
    574       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     591      if ((*this)[idx].IsExcluded()
     592          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     593          || !(*this)[idx].IsFFactorMethodValid())
    575594        return kFALSE;
    576595      val = (*this)[idx].GetPheFFactorMethodErr();
    577596      break;
    578597    case 11:
    579       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     598      if ((*this)[idx].IsExcluded()
     599          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     600          || !(*this)[idx].IsFFactorMethodValid())
    580601        return kFALSE;
    581602      val = (*this)[idx].GetMeanConversionFFactorMethod();
    582603      break;
    583604    case 12:
    584       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     605      if ((*this)[idx].IsExcluded()
     606          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     607          || !(*this)[idx].IsFFactorMethodValid())
    585608        return kFALSE;
    586609      val = (*this)[idx].GetConversionFFactorMethodErr();
    587610      break;
    588611    case 13:
    589       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     612      if ((*this)[idx].IsExcluded()
     613          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     614          || !(*this)[idx].IsFFactorMethodValid())
    590615        return kFALSE;
    591616      val = (*this)[idx].GetTotalFFactorFFactorMethod();
    592617      break;
    593618    case 14:
    594       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsFFactorMethodValid())
     619      if ((*this)[idx].IsExcluded()
     620          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     621          || !(*this)[idx].IsFFactorMethodValid())
    595622        return kFALSE;
    596623      val = (*this)[idx].GetTotalFFactorErrFFactorMethod();
    597624      break;
    598625    case 15:
    599       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     626      if ((*this)[idx].IsExcluded()
     627          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     628          || !(*this)[idx].IsBlindPixelMethodValid())
    600629        return kFALSE;
    601630      val = fBlindPixel->GetMeanFluxInsidePlexiglass()*area;
    602631      break;
    603632    case 16:
    604       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     633      if ((*this)[idx].IsExcluded()
     634          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     635          || !(*this)[idx].IsBlindPixelMethodValid())
    605636        return kFALSE;
    606637      val = fBlindPixel->GetMeanFluxErrInsidePlexiglass()*area;
    607638      break;
    608639    case 17:
    609       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     640      if ((*this)[idx].IsExcluded()
     641          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     642          || !(*this)[idx].IsBlindPixelMethodValid())
    610643        return kFALSE;
    611644      val = (*this)[idx].GetMeanConversionBlindPixelMethod();
    612645      break;
    613646    case 18:
    614       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     647      if ((*this)[idx].IsExcluded()
     648          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     649          || !(*this)[idx].IsBlindPixelMethodValid())
    615650        return kFALSE;
    616651      val = (*this)[idx].GetConversionBlindPixelMethodErr();
    617652      break;
    618653    case 19:
    619       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     654      if ((*this)[idx].IsExcluded()
     655          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     656          || !(*this)[idx].IsBlindPixelMethodValid())
    620657        return kFALSE;
    621658      val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
    622659      break;
    623660    case 20:
    624       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsBlindPixelMethodValid())
     661      if ((*this)[idx].IsExcluded()
     662          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     663          || !(*this)[idx].IsBlindPixelMethodValid())
    625664        return kFALSE;
    626665      val = (*this)[idx].GetTotalFFactorErrBlindPixelMethod();
    627666      break;
    628667    case 21:
    629       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     668      if ((*this)[idx].IsExcluded()
     669          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     670          || !(*this)[idx].IsPINDiodeMethodValid())
    630671        return kFALSE;
    631672      val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area;
    632673      break;
    633674    case 22:
    634       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     675      if ((*this)[idx].IsExcluded()
     676          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     677          || !(*this)[idx].IsPINDiodeMethodValid())
    635678        return kFALSE;
    636679      val = fPINDiode->GetMeanFluxErrOutsidePlexiglass()*area;
    637680      break;
    638681    case 23:
    639       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     682      if ((*this)[idx].IsExcluded()
     683          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     684          || !(*this)[idx].IsPINDiodeMethodValid())
    640685        return kFALSE;
    641686      val = (*this)[idx].GetMeanConversionPINDiodeMethod();
    642687      break;
    643688    case 24:
    644       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     689      if ((*this)[idx].IsExcluded()
     690          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     691          || !(*this)[idx].IsPINDiodeMethodValid())
    645692        return kFALSE;
    646693      val = (*this)[idx].GetConversionPINDiodeMethodErr();
    647694      break;
    648695    case 25:
    649       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     696      if ((*this)[idx].IsExcluded()
     697          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     698          || !(*this)[idx].IsPINDiodeMethodValid())
    650699        return kFALSE;
    651700      val = (*this)[idx].GetTotalFFactorPINDiodeMethod();
    652701      break;
    653702    case 26:
    654       if ((*this)[idx].IsExcluded() || !(*fBadPixels)[idx].IsCalibrationSignalOK() || !(*this)[idx].IsPINDiodeMethodValid())
     703      if ((*this)[idx].IsExcluded()
     704          || !(*fBadPixels)[idx].IsCalibrationSignalOK()
     705          || !(*this)[idx].IsPINDiodeMethodValid())
    655706        return kFALSE;
    656707      val = (*this)[idx].GetTotalFFactorErrPINDiodeMethod();
    … …  
    782833    (*this)[idx].DrawClone();
    783834}
     835
     836//
     837// Calculate the weighted mean of the phe's of all inner and outer pixels, respectively.
     838// Bad pixels are excluded from the calculation.
     839//
     840Bool_t MCalibrationChargeCam::CalcMeanFluxPhotonsFFactorMethod()
     841{
     842
     843  const Float_t avQERelErrSquare             = fAverageQEErr * fAverageQEErr
     844                                            / (fAverageQE    * fAverageQE );
     845
     846  Float_t sumweightsinner = 0.;
     847  Float_t sumphesinner    = 0.;
     848  Float_t sumweightsouter = 0.;
     849  Float_t sumphesouter    = 0.;
     850     
     851  TIter Next(fPixels);
     852  MCalibrationChargePix *pix;
     853  while ((pix=(MCalibrationChargePix*)Next()))
     854    {
     855
     856      if (!pix->IsFFactorMethodValid())
     857        continue;
     858
     859      const Int_t idx = pix->GetPixId();
     860
     861      if(!(*fBadPixels)[idx].IsCalibrationResultOK())
     862        continue;
     863     
     864      const Float_t nphe    = pix->GetPheFFactorMethod();
     865      const Float_t npheerr = pix->GetPheFFactorMethodErr();
     866      const Float_t ratio   = fGeomCam->GetPixRatio(idx);
     867
     868      if (npheerr > 0.)
     869        {
     870          //
     871          // first the inner pixels:
     872          //
     873          if (ratio == 1.)
     874            {
     875              const Float_t weight = 1./npheerr/npheerr;
     876              sumweightsinner += weight;
     877              sumphesinner    += weight*nphe;
     878            }
     879          else
     880            {
     881              //
     882              // now the outers
     883              //
     884              const Float_t weight = 1./npheerr/npheerr;
     885              sumweightsouter += weight;
     886              sumphesouter    += weight*nphe;
     887            }
     888        } /* if npheerr != 0 */
     889    } /* while ((pix=(MCalibrationChargePix*)Next())) */
     890 
     891  if (sumweightsinner <= 0. || sumphesinner <= 0.)
     892    {
     893      *fLog << warn << " Mean number of phe's from inner pixels cannot be calculated: "
     894            << " Sum of weights: " << sumweightsinner
     895            << " Sum of weighted phes: " << sumphesinner << endl;
     896      return kFALSE;
     897    }
     898  else
     899    {
     900      fMeanFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
     901      fMeanFluxPhesInnerPixelErr = TMath::Sqrt(1./sumweightsinner);
     902
     903    }
     904
     905  if (sumweightsouter <= 0. || sumphesouter <= 0.)
     906    {
     907      *fLog << warn << " Mean number of phe's from outer pixels cannot be calculated: "
     908            << " Sum of weights or sum of weighted phes is 0. " << endl;
     909    }
     910  else
     911    {
     912      fMeanFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
     913      fMeanFluxPhesOuterPixelErr = TMath::Sqrt(1./sumweightsouter);
     914    }
     915
     916
     917  const Float_t meanFluxPhotonsRelErrSquare = fMeanFluxPhesInnerPixelErr * fMeanFluxPhesInnerPixelErr
     918                                           / (fMeanFluxPhesInnerPixel    * fMeanFluxPhesInnerPixel);
     919
     920  fMeanFluxPhotonsInnerPixel    =  fMeanFluxPhesInnerPixel/fAverageQE;
     921  fMeanFluxPhotonsInnerPixelErr =  TMath::Sqrt(meanFluxPhotonsRelErrSquare + avQERelErrSquare)
     922                                 * fMeanFluxPhotonsInnerPixel;
     923
     924  fMeanFluxPhotonsOuterPixel    = 4.*fMeanFluxPhotonsInnerPixel;
     925  fMeanFluxPhotonsOuterPixelErr = 4.*fMeanFluxPhotonsInnerPixelErr; 
     926
     927  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
     928        << fMeanFluxPhesInnerPixel << " +- " << fMeanFluxPhesInnerPixelErr << endl;
     929
     930  *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): "
     931        << fMeanFluxPhotonsInnerPixel << " +- " << fMeanFluxPhotonsInnerPixelErr << endl;
     932
     933 
     934
     935  return kTRUE;
     936}
     937
     938void MCalibrationChargeCam::ApplyFFactorCalibration()
     939{
     940
     941  const Float_t meanphotRelErrSquare        = fMeanFluxPhotonsInnerPixelErr * fMeanFluxPhotonsInnerPixelErr
     942                                           /( fMeanFluxPhotonsInnerPixel    * fMeanFluxPhotonsInnerPixel );
     943
     944  TIter Next(fPixels);
     945  MCalibrationChargePix *pix;
     946  while ((pix=(MCalibrationChargePix*)Next()))
     947    {
     948
     949      if (!pix->IsFFactorMethodValid())
     950        continue;
     951
     952      const Int_t idx = pix->GetPixId();
     953
     954      if(!(*fBadPixels)[idx].IsCalibrationResultOK())
     955        continue;
     956
     957      const Float_t ratio   = fGeomCam->GetPixRatio(idx);
     958      //
     959      // Calculate the conversion factor between PHOTONS and FADC counts
     960      //
     961      // Nphot = Nphe / avQE
     962      // conv  = Nphot / FADC counts
     963      //
     964      Float_t conv;
     965     
     966      if (ratio == 1.)
     967        conv               = fMeanFluxPhotonsInnerPixel / pix->GetMeanCharge();
     968      else
     969        conv               = fMeanFluxPhotonsOuterPixel / pix->GetMeanCharge();
     970
     971      if (conv <= 0.)
     972        {
     973          pix->SetFFactorMethodValid(kFALSE);
     974          continue;
     975        }
     976     
     977      const Float_t chargeRelErrSquare       =   pix->GetMeanChargeErr() * pix->GetMeanChargeErr()
     978                                             / ( pix->GetMeanCharge()    * pix->GetMeanCharge());
     979      const Float_t rsigmaChargeRelErrSquare =   pix->GetRSigmaChargeErr() *  pix->GetRSigmaChargeErr()
     980                                              / (pix->GetRSigmaCharge()    * pix->GetRSigmaCharge()) ;
     981
     982      const Float_t convrelerr = TMath::Sqrt(meanphotRelErrSquare + chargeRelErrSquare);
     983
     984      if (convrelerr > fConvFFactorRelErrLimit)
     985        {
     986          *fLog << warn << GetDescriptor() << ": Conversion Factor F-Factor Method Rel. Error: "
     987                << convrelerr << " above limit of: " << fConvFFactorRelErrLimit
     988                << " in pixel: " << idx << endl;
     989          pix->SetFFactorMethodValid(kFALSE);
     990          continue;
     991        }
     992     
     993      //
     994      // Calculate the Total F-Factor of the camera (in photons)
     995      //
     996      const Float_t totalFFactor  =  (pix->GetRSigmaCharge()/pix->GetMeanCharge())
     997                                    *TMath::Sqrt(fMeanFluxPhotonsInnerPixel);
     998     
     999      //
     1000      // Calculate the error of the Total F-Factor of the camera ( in photons )
     1001      //
     1002      const Float_t totalFFactorErr = TMath::Sqrt(  rsigmaChargeRelErrSquare
     1003                                    + chargeRelErrSquare
     1004                                    + meanphotRelErrSquare );
     1005
     1006      pix->SetConversionFFactorMethod(conv,
     1007                                      convrelerr*conv,
     1008                                      totalFFactor*TMath::Sqrt(conv));
     1009
     1010      pix->SetTotalFFactorFFactorMethod(   totalFFactor   );
     1011      pix->SetTotalFFactorErrFFactorMethod(totalFFactorErr);     
     1012      pix->SetFFactorMethodValid();
     1013    }
     1014}
     1015
     1016
    7841017
    7851018void MCalibrationChargeCam::ApplyBlindPixelCalibration()
Note: See TracChangeset for help on using the changeset viewer.