Changeset 3511


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

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3510 r3511  
    4242     of sigmas to calculation of variances which saves all the useless
    4343     square rooting.
     44
     45   - took out pointers to MCalibraitonChargeBlindPix and
     46     MCalibrationChargePINDiode in MCalibrationChargeCam.
    4447
    4548

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

    r3496 r3511  
    279279        return kFALSE;
    280280
    281     fCam->SetPINDiode(fPINDiode);
    282     fCam->SetBlindPixel(fBlindPixel);
    283 
    284281    fEvtTime = (MTime*)pList->FindObject("MTime");
    285282
    … …  
    558555      {
    559556          fCam->SetBlindPixelMethodValid(kTRUE);
    560           fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels);
     557          fCam->ApplyBlindPixelCalibration(*fGeom,*fBadPixels, *fBlindPixel);
    561558      }
    562559  }
    … …  
    582579      {
    583580          fCam->SetPINDiodeMethodValid(kTRUE);
    584           fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels);
     581          fCam->ApplyPINDiodeCalibration(*fGeom,*fBadPixels, *fPINDiode);
    585582      }
    586583  }

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

    r3496 r3511  
    141141const Float_t MCalibrationChargeCam::gkAverageQE                = 0.18;     
    142142const Float_t MCalibrationChargeCam::gkAverageQEErr             = 0.02; 
    143 const Float_t MCalibrationChargeCam::fgConvFFactorRelErrLimit   = 0.25;
    144 const Float_t MCalibrationChargeCam::fgPheFFactorRelLimit       = 7.;
     143const Float_t MCalibrationChargeCam::fgConvFFactorRelErrLimit   = 0.35;
     144const Float_t MCalibrationChargeCam::fgPheFFactorRelErrLimit    = 5.;
    145145// --------------------------------------------------------------------------
    146146//
    … …  
    170170    SetAverageQE();
    171171    SetConvFFactorRelErrLimit();
    172     SetPheFFactorRelLimit();
     172    SetPheFFactorRelErrLimit();
    173173}
    174174
    … …  
    263263  fNumExcludedPixels            = 0;
    264264  fMeanFluxPhesInnerPixel       = 0.;
    265   fMeanFluxPhesInnerPixelErr    = 0.;
     265  fMeanFluxPhesInnerPixelVar    = 0.;
    266266  fMeanFluxPhesOuterPixel       = 0.;
    267   fMeanFluxPhesOuterPixelErr    = 0.;
     267  fMeanFluxPhesOuterPixelVar    = 0.;
    268268
    269269  CLRBIT(fFlags,kBlindPixelMethodValid);
    … …  
    287287{
    288288    b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid);
     289}
     290
     291Float_t MCalibrationChargeCam::GetMeanFluxPhesInnerPixelErr()   const
     292{
     293  if (fMeanFluxPhesInnerPixelVar <= 0.)
     294    return -1.;
     295  return TMath::Sqrt(fMeanFluxPhesInnerPixelVar);
     296}
     297
     298Float_t MCalibrationChargeCam::GetMeanFluxPhesOuterPixelErr()   const
     299{
     300  if (fMeanFluxPhesOuterPixelVar <= 0.)
     301    return -1.;
     302  return TMath::Sqrt(fMeanFluxPhesOuterPixelVar);
     303}
     304
     305Float_t MCalibrationChargeCam::GetMeanFluxPhotonsInnerPixelErr()   const
     306{
     307  if (fMeanFluxPhotonsInnerPixelVar <= 0.)
     308    return -1.;
     309  return TMath::Sqrt(fMeanFluxPhotonsInnerPixelVar);
     310}
     311
     312Float_t MCalibrationChargeCam::GetMeanFluxPhotonsOuterPixelErr()   const
     313{
     314  if (fMeanFluxPhotonsOuterPixelVar <= 0.)
     315    return -1.;
     316  return TMath::Sqrt(fMeanFluxPhotonsOuterPixelVar);
    289317}
    290318
    … …  
    506534        return kFALSE;
    507535      if ((*this)[idx].GetRSigmaCharge() == -1.)
     536        return kFALSE;
     537      if ((*this)[idx].GetMeanCharge() == 0.)
     538        return kFALSE;
     539      val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetMeanCharge();
     540      break;
     541    case 8:
     542      if ((*this)[idx].IsExcluded())
     543        return kFALSE;
     544      if ((*this)[idx].GetRSigmaCharge() <= 0.)
    508545          return kFALSE;
    509       val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetMeanCharge();
    510       break;
    511     case 8:
    512       if ((*this)[idx].IsExcluded())
    513         return kFALSE;
    514       if ((*this)[idx].GetRSigmaCharge() == -1.)
     546      if ((*this)[idx].GetMeanCharge() <= 0.)
     547        return kFALSE;
     548      if ((*this)[idx].GetRSigmaChargeErr() <= 0.)
    515549          return kFALSE;
     550      if ((*this)[idx].GetMeanChargeErr() <= 0.)
     551        return kFALSE;
    516552      // relative error RsigmaCharge square
    517553      val =    (*this)[idx].GetRSigmaChargeErr()* (*this)[idx].GetRSigmaChargeErr()
    … …  
    553589      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsFFactorMethodValid())
    554590        return kFALSE;
    555       val = (*this)[idx].GetTotalFFactorErrFFactorMethod();
     591      val = (*this)[idx].GetTotalFFactorFFactorMethodErr();
    556592      break;
    557593    case 15:
    558594      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
    559595        return kFALSE;
    560       val = fBlindPixel->GetMeanFluxInsidePlexiglass()*area;
     596      //      val = fBlindPixel->GetMeanFluxInsidePlexiglass()*area;
     597      val = 1.;
    561598      break;
    562599    case 16:
    563600      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
    564601        return kFALSE;
    565       val = fBlindPixel->GetMeanFluxErrInsidePlexiglass()*area;
     602      //      val = fBlindPixel->GetMeanFluxErrInsidePlexiglass()*area;
     603      val = 1.;
    566604      break;
    567605    case 17:
    … …  
    583621      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsBlindPixelMethodValid())
    584622        return kFALSE;
    585       val = (*this)[idx].GetTotalFFactorErrBlindPixelMethod();
     623      val = (*this)[idx].GetTotalFFactorBlindPixelMethodErr();
    586624      break;
    587625    case 21:
    588626      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
    589627        return kFALSE;
    590       val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area;
     628      //      val = fPINDiode->GetMeanFluxOutsidePlexiglass()*area;
     629      val = 1.;
    591630      break;
    592631    case 22:
    593632      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
    594633        return kFALSE;
    595       val = fPINDiode->GetMeanFluxErrOutsidePlexiglass()*area;
     634      //      val = fPINDiode->GetMeanFluxErrOutsidePlexiglass()*area;
     635      val = 1.;
    596636      break;
    597637    case 23:
    … …  
    613653      if ((*this)[idx].IsExcluded() || !(*this)[idx].IsPINDiodeMethodValid())
    614654        return kFALSE;
    615       val = (*this)[idx].GetTotalFFactorErrPINDiodeMethod();
     655      val = (*this)[idx].GetTotalFFactorPINDiodeMethodErr();
    616656      break;
    617657    case 27:
    … …  
    725765// which are fPheFFactorRelLimit sigmas from the mean.
    726766//
    727 Bool_t MCalibrationChargeCam::CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, const MBadPixelsCam &bad)
    728 {
    729 
    730   const Float_t avQERelErrSquare             = fAverageQEErr * fAverageQEErr
    731                                             / (fAverageQE    * fAverageQE );
     767Bool_t MCalibrationChargeCam::CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, MBadPixelsCam &bad)
     768{
     769
     770  const Float_t avQERelVar = fAverageQEVar / (fAverageQE * fAverageQE );
    732771
    733772  Float_t sumweightsinner = 0.;
    … …  
    757796      if (npheerr > 0.)
    758797        {
     798          const Float_t weight = nphe*nphe/npheerr/npheerr;
    759799          //
    760800          // first the inner pixels:
    … …  
    762802          if (ratio == 1.)
    763803            {
    764               const Float_t weight = 1./npheerr/npheerr;
    765804              sumweightsinner += weight;
    766805              sumphesinner    += weight*nphe;
    … …  
    772811              // now the outers
    773812              //
    774               const Float_t weight = 1./npheerr/npheerr;
    775813              sumweightsouter += weight;
    776814              sumphesouter    += weight*nphe;
    … …  
    790828    {
    791829      fMeanFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
    792       fMeanFluxPhesInnerPixelErr = TMath::Sqrt(1./sumweightsinner);
    793 
     830      fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel;
    794831    }
    795832
    … …  
    802839    {
    803840      fMeanFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
    804       fMeanFluxPhesOuterPixelErr = TMath::Sqrt(1./sumweightsouter);
    805     }
    806 
    807   Float_t meanFluxPhotonsRelErrSquare = fMeanFluxPhesInnerPixelErr * fMeanFluxPhesInnerPixelErr
    808                                      / (fMeanFluxPhesInnerPixel    * fMeanFluxPhesInnerPixel);
     841      fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel;
     842    }
     843
     844  Float_t meanFluxPhotonsRelVar  = fMeanFluxPhesInnerPixelVar
     845                                / (fMeanFluxPhesInnerPixel * fMeanFluxPhesInnerPixel);
    809846
    810847  fMeanFluxPhotonsInnerPixel    =  fMeanFluxPhesInnerPixel/fAverageQE;
    811   fMeanFluxPhotonsInnerPixelErr =  TMath::Sqrt(meanFluxPhotonsRelErrSquare + avQERelErrSquare)
    812                                              * fMeanFluxPhotonsInnerPixel;
    813 
    814   fMeanFluxPhotonsOuterPixel    = 4.*fMeanFluxPhotonsInnerPixel;
    815   fMeanFluxPhotonsOuterPixelErr = 4.*fMeanFluxPhotonsInnerPixelErr; 
     848  fMeanFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
     849                                 * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel;
     850
     851  fMeanFluxPhotonsOuterPixel    = 4. *fMeanFluxPhotonsInnerPixel;
     852  fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar; 
    816853 
     854  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
     855        << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl;
     856
     857  *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): "
     858        << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl;
     859
    817860  //
    818861  // Here starts the second loop discarting pixels out of the range:
    819862  //
    820   const Float_t innererr = TMath::Sqrt((Float_t)validinner)*fPheFFactorRelLimit*fMeanFluxPhesInnerPixelErr;
    821   const Float_t outererr = TMath::Sqrt((Float_t)validouter)*fPheFFactorRelLimit*fMeanFluxPhesOuterPixelErr;
     863  const Float_t innervar = (Float_t)validinner*fPheFFactorRelVarLimit*fMeanFluxPhesInnerPixelVar;
     864  const Float_t outervar = (Float_t)validouter*fPheFFactorRelVarLimit*fMeanFluxPhesOuterPixelVar;
     865 
     866  Float_t innererr;
     867  Float_t outererr;
     868
     869  if (innervar > 0.)
     870    innererr = TMath::Sqrt(innervar);
     871  else
     872    {
     873      *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for inner pixels "
     874            << " is smaller than 0. " << endl;
     875      SetFFactorMethodValid(kFALSE);
     876      return kFALSE;
     877    }
     878 
     879  if (outervar > 0.)
     880    outererr = TMath::Sqrt(outervar);
     881  else
     882    {
     883      *fLog << warn << GetDescriptor() << ": Variance of mean number of photo-electrons for outer pixels "
     884            << " is smaller than 0. " << endl;
     885      SetFFactorMethodValid(kFALSE);
     886      return kFALSE;
     887    }
    822888
    823889  const Float_t lowerpheinnerlimit = fMeanFluxPhesInnerPixel - innererr;
    … …  
    851917      if (npheerr > 0.)
    852918        {
     919          const Float_t weight = nphe*nphe/npheerr/npheerr;
    853920          //
    854921          // first the inner pixels:
    … …  
    859926              if (nphe < lowerpheinnerlimit || nphe > upperpheinnerlimit)
    860927                {
    861                   pix2->SetFFactorMethodValid(kFALSE);
     928                  bad[idx].SetUnsuitable(MBadPixelsPix::kUnreliableRun);
    862929                  continue;
    863930                }
    864931
    865               const Float_t weight = 1./npheerr/npheerr;
    866932              sumweightsinner += weight;
    867933              sumphesinner    += weight*nphe;
    … …  
    874940              if (nphe < lowerpheouterlimit || nphe > upperpheouterlimit)
    875941                {
    876                   pix2->SetFFactorMethodValid(kFALSE);
     942                  bad[idx].SetUnsuitable(MBadPixelsPix::kUnreliableRun);
    877943                  continue;
    878944                }
    879945
    880               const Float_t weight = 1./npheerr/npheerr;
    881946              sumweightsouter += weight;
    882947              sumphesouter    += weight*nphe;
    … …  
    895960    {
    896961      fMeanFluxPhesInnerPixel    = sumphesinner/sumweightsinner;
    897       fMeanFluxPhesInnerPixelErr = TMath::Sqrt(1./sumweightsinner);
     962      fMeanFluxPhesInnerPixelVar = (1./sumweightsinner)*fMeanFluxPhesInnerPixel*fMeanFluxPhesInnerPixel;
    898963
    899964    }
    … …  
    907972    {
    908973      fMeanFluxPhesOuterPixel    = sumphesouter/sumweightsouter;
    909       fMeanFluxPhesOuterPixelErr = TMath::Sqrt(1./sumweightsouter);
    910     }
    911 
    912   meanFluxPhotonsRelErrSquare = fMeanFluxPhesInnerPixelErr * fMeanFluxPhesInnerPixelErr
    913                              / (fMeanFluxPhesInnerPixel    * fMeanFluxPhesInnerPixel);
     974      fMeanFluxPhesOuterPixelVar = (1./sumweightsouter)*fMeanFluxPhesOuterPixel*fMeanFluxPhesOuterPixel;
     975    }
     976
     977  meanFluxPhotonsRelVar = fMeanFluxPhesInnerPixelVar
     978                       / (fMeanFluxPhesInnerPixel    * fMeanFluxPhesInnerPixel);
    914979
    915980  fMeanFluxPhotonsInnerPixel    =  fMeanFluxPhesInnerPixel/fAverageQE;
    916   fMeanFluxPhotonsInnerPixelErr =  TMath::Sqrt(meanFluxPhotonsRelErrSquare + avQERelErrSquare)
    917                                  * fMeanFluxPhotonsInnerPixel;
    918 
    919   fMeanFluxPhotonsOuterPixel    = 4.*fMeanFluxPhotonsInnerPixel;
    920   fMeanFluxPhotonsOuterPixelErr = 4.*fMeanFluxPhotonsInnerPixelErr; 
     981  fMeanFluxPhotonsInnerPixelVar =  (meanFluxPhotonsRelVar + avQERelVar)
     982                                 * fMeanFluxPhotonsInnerPixel * fMeanFluxPhotonsInnerPixel;
     983
     984  fMeanFluxPhotonsOuterPixel    = 4. *fMeanFluxPhotonsInnerPixel;
     985  fMeanFluxPhotonsOuterPixelVar = 16.*fMeanFluxPhotonsInnerPixelVar; 
    921986
    922987  *fLog << inf << " Mean number of photo-electrons from inner pixels (F-Factor Method): "
    923         << fMeanFluxPhesInnerPixel << " +- " << fMeanFluxPhesInnerPixelErr << endl;
     988        << fMeanFluxPhesInnerPixel << " +- " << GetMeanFluxPhesInnerPixelErr() << endl;
    924989
    925990  *fLog << inf << " Mean number of photons from inner pixels (F-Factor Method): "
    926         << fMeanFluxPhotonsInnerPixel << " +- " << fMeanFluxPhotonsInnerPixelErr << endl;
     991        << fMeanFluxPhotonsInnerPixel << " +- " << GetMeanFluxPhotonsInnerPixelErr() << endl;
    927992
    928993  return kTRUE;
    … …  
    932997{
    933998
    934   const Float_t meanphotRelErrSquare        = fMeanFluxPhotonsInnerPixelErr * fMeanFluxPhotonsInnerPixelErr
    935                                            /( fMeanFluxPhotonsInnerPixel    * fMeanFluxPhotonsInnerPixel );
     999  const Float_t meanphotRelVar  = fMeanFluxPhotonsInnerPixelVar
     1000                               /( fMeanFluxPhotonsInnerPixel    * fMeanFluxPhotonsInnerPixel );
    9361001
    9371002  TIter Next(fPixels);
    … …  
    9711036        }
    9721037     
    973       const Float_t chargeRelErrSquare       =   pix->GetMeanChargeErr() * pix->GetMeanChargeErr()
    974                                              / ( pix->GetMeanCharge()    * pix->GetMeanCharge());
    975       const Float_t rsigmaChargeRelErrSquare =   pix->GetRSigmaChargeErr() *  pix->GetRSigmaChargeErr()
    976                                               / (pix->GetRSigmaCharge()    * pix->GetRSigmaCharge()) ;
    977 
    978       const Float_t convrelerr = TMath::Sqrt(meanphotRelErrSquare + chargeRelErrSquare);
    979 
    980       if (convrelerr > fConvFFactorRelErrLimit)
     1038      const Float_t chargeRelVar       =   pix->GetMeanChargeErr() * pix->GetMeanChargeErr()
     1039                                       / ( pix->GetMeanCharge()    * pix->GetMeanCharge());
     1040      const Float_t rsigmaChargeRelVar =   pix->GetRSigmaChargeErr() *  pix->GetRSigmaChargeErr()
     1041                                        / (pix->GetRSigmaCharge()    * pix->GetRSigmaCharge()) ;
     1042
     1043      const Float_t convrelvar = meanphotRelVar + chargeRelVar;
     1044
     1045      if (convrelvar > fConvFFactorRelVarLimit)
    9811046        {
    982           *fLog << warn << GetDescriptor() << ": Conversion Factor F-Factor Method Rel. Error: "
    983                 << convrelerr << " above limit of: " << fConvFFactorRelErrLimit
     1047          *fLog << warn << GetDescriptor() << ": Conversion Factor F-Factor Method Rel. Variance: "
     1048                << convrelvar << " above limit of: " << fConvFFactorRelVarLimit
    9841049                << " in pixel: " << idx << endl;
    9851050          pix->SetFFactorMethodValid(kFALSE);
    … …  
    9961061      // Calculate the error of the Total F-Factor of the camera ( in photons )
    9971062      //
    998       const Float_t totalFFactorErr = TMath::Sqrt(  rsigmaChargeRelErrSquare
    999                                     + chargeRelErrSquare
    1000                                     + meanphotRelErrSquare );
    1001 
    1002       pix->SetConversionFFactorMethod(conv,
    1003                                       convrelerr*conv,
    1004                                       totalFFactor*TMath::Sqrt(conv));
     1063      const Float_t totalFFactorVar = rsigmaChargeRelVar + chargeRelVar + meanphotRelVar;
     1064
     1065      if (convrelvar > 0. && conv > 0.)
     1066        pix->SetConversionFFactorMethod(conv,
     1067                                        TMath::Sqrt(convrelvar)*conv,
     1068                                        totalFFactor*TMath::Sqrt(conv));
    10051069
    10061070      pix->SetTotalFFactorFFactorMethod(   totalFFactor   );
    1007       pix->SetTotalFFactorErrFFactorMethod(totalFFactorErr);     
     1071
     1072      if (totalFFactorVar > 0.)
     1073        pix->SetTotalFFactorFFactorMethodErr(TMath::Sqrt(totalFFactorVar));     
     1074
    10081075      pix->SetFFactorMethodValid();
    10091076    }
    … …  
    10121079
    10131080
    1014 void MCalibrationChargeCam::ApplyBlindPixelCalibration(const MGeomCam &geom, const MBadPixelsCam &bad)
    1015 {
    1016 
    1017   Float_t flux    = fBlindPixel->GetMeanFluxInsidePlexiglass();
    1018   Float_t fluxerr = fBlindPixel->GetMeanFluxErrInsidePlexiglass();
     1081void MCalibrationChargeCam::ApplyBlindPixelCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargeBlindPix &blindpix)
     1082{
     1083
     1084  Float_t flux    = blindpix.GetMeanFluxInsidePlexiglass();
     1085  Float_t fluxerr = blindpix.GetMeanFluxErrInsidePlexiglass();
    10191086
    10201087  TIter Next(fPixels);
    … …  
    10561123}
    10571124
    1058 void MCalibrationChargeCam::ApplyPINDiodeCalibration(const MGeomCam &geom, const MBadPixelsCam &bad)
    1059 {
    1060 
    1061   Float_t flux    = fPINDiode->GetMeanFluxOutsidePlexiglass();
    1062   Float_t fluxerr = fPINDiode->GetMeanFluxErrOutsidePlexiglass();
     1125void MCalibrationChargeCam::ApplyPINDiodeCalibration(const MGeomCam &geom, const MBadPixelsCam &bad, const MCalibrationChargePINDiode &pindiode)
     1126{
     1127
     1128  Float_t flux    = pindiode.GetMeanFluxOutsidePlexiglass();
     1129  Float_t fluxerr = pindiode.GetMeanFluxErrOutsidePlexiglass();
    10631130
    10641131  TIter Next(fPixels);
    … …  
    11371204// and the number of photons reaching the plexiglass for one Inner Pixel
    11381205//
    1139 // FIXME: The PINDiode is still not working and so is the code
    1140 //
    11411206Bool_t MCalibrationChargeCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
    11421207{
    … …  
    11571222// and the number of photons reaching one Inner Pixel.
    11581223// The procedure is not yet defined.
    1159 //
    1160 // FIXME: The PINDiode is still not working and so is the code
    11611224//
    11621225Bool_t MCalibrationChargeCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)

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

    r3479 r3511  
    2626 
    2727  static const Float_t fgConvFFactorRelErrLimit; // The limit for acceptance of the rel. error of the conversion factor with the FFactor method
    28   static const Float_t fgPheFFactorRelLimit;     // The rel. limit for acceptance of a calculated number of phe's w.r.t the mean number (in sigma of the error)
     28  static const Float_t fgPheFFactorRelErrLimit;  // The rel. limit for acceptance of a calculated number of phe's w.r.t the mean number (in sigma of the error)
    2929 
    3030  Float_t fAverageQE;                // The average quantum efficieny (see Class description)
    31   Float_t fAverageQEErr;             // The error of the average quantum efficieny (see Class description)
     31  Float_t fAverageQEVar;             // The error of the average quantum efficieny (see Class description)
    3232 
    33   Float_t fConvFFactorRelErrLimit;   // The limit for acceptance of the rel. error of the conversion factor with the FFactor method
    34   Float_t fPheFFactorRelLimit;       // The rel. limit for acceptance of a calculated number of phe's w.r.t the mean number (in sigma of the error).
     33  Float_t fConvFFactorRelVarLimit;   // The limit for acceptance of the rel. error of the conversion factor with the FFactor method
     34  Float_t fPheFFactorRelVarLimit;    // The rel. limit for acceptance of a calculated number of phe's w.r.t the mean number (in variances of the error).
    3535 
    3636  Int_t fNumPixels;
    … …  
    4343  MBadPixelsPix         *fAverageOuterBadPix;    //-> Average Pixel of all events
    4444 
    45   const MCalibrationChargeBlindPix *fBlindPixel; //! Pointer to the Blind Pixel with fit results
    46   const MCalibrationChargePINDiode *fPINDiode;   //! Pointer to the PIN Diode with fit results
    47 
    4845  TH1D* fOffsets;                                //!
    4946  TH1D* fSlopes;                                 //!
    … …  
    5855
    5956  Float_t fMeanFluxPhesInnerPixel;        //  The mean number of photo-electrons in an INNER PIXEL
    60   Float_t fMeanFluxPhesInnerPixelErr;     //  The uncertainty about the number of photo-electrons INNER PIXEL 
     57  Float_t fMeanFluxPhesInnerPixelVar;     //  The variance of the number of photo-electrons INNER PIXEL 
    6158  Float_t fMeanFluxPhesOuterPixel;        //  The mean number of photo-electrons in an INNER PIXEL
    62   Float_t fMeanFluxPhesOuterPixelErr;     //  The uncertainty about the number of photo-electrons INNER PIXEL 
     59  Float_t fMeanFluxPhesOuterPixelVar;     //  The variance of the number of photo-electrons INNER PIXEL 
    6360
    6461  Float_t fMeanFluxPhotonsInnerPixel;        //  The mean number of photo-electrons in an INNER PIXEL
    65   Float_t fMeanFluxPhotonsInnerPixelErr;     //  The uncertainty about the number of photo-electrons INNER PIXEL 
     62  Float_t fMeanFluxPhotonsInnerPixelVar;     //  The variance of the number of photo-electrons INNER PIXEL 
    6663  Float_t fMeanFluxPhotonsOuterPixel;        //  The mean number of photo-electrons in an INNER PIXEL
    67   Float_t fMeanFluxPhotonsOuterPixelErr;     //  The uncertainty about the number of photo-electrons INNER PIXEL 
     64  Float_t fMeanFluxPhotonsOuterPixelVar;     //  The variance of the number of photo-electrons INNER PIXEL 
    6865 
    6966public:
    … …  
    7875  void SetAverageQE(          const Float_t qe= gkAverageQE,
    7976                              const Float_t err=gkAverageQEErr)         { fAverageQE    = qe;           
    80                                                                           fAverageQEErr = err;        }
     77                                                                          fAverageQEVar = err*err; }
    8178  void SetNumPixelsExcluded(  const UInt_t n )            {  fNumExcludedPixels = n; }
    82   void SetConvFFactorRelErrLimit( const Float_t f=fgConvFFactorRelErrLimit ) { fConvFFactorRelErrLimit = f; }
    83   void SetPheFFactorRelLimit (  const Float_t f=fgPheFFactorRelLimit )     { fPheFFactorRelLimit = f;    } 
    84 
    85   void SetPINDiode  ( const MCalibrationChargePINDiode *d ) {  fPINDiode   = d;      }
    86   void SetBlindPixel( const MCalibrationChargeBlindPix *b ) {  fBlindPixel = b;      }
     79  void SetConvFFactorRelErrLimit( const Float_t f=fgConvFFactorRelErrLimit ) { fConvFFactorRelVarLimit = f*f; }
     80  void SetPheFFactorRelErrLimit (  const Float_t f=fgPheFFactorRelErrLimit )     { fPheFFactorRelVarLimit = f*f;    } 
    8781
    8882  void SetFFactorMethodValid(    const Bool_t b = kTRUE );
    … …  
    10094
    10195  Float_t GetMeanFluxPhesInnerPixel()     const { return fMeanFluxPhesInnerPixel;     }
    102   Float_t GetMeanFluxPhesInnerPixelErr()  const { return fMeanFluxPhesInnerPixelErr;  }
     96  Float_t GetMeanFluxPhesInnerPixelErr()  const;
    10397  Float_t GetMeanFluxPhesOuterPixel()     const { return fMeanFluxPhesOuterPixel;     }
    104   Float_t GetMeanFluxPhesOuterPixelErr()  const { return fMeanFluxPhesOuterPixelErr;  }
     98  Float_t GetMeanFluxPhesOuterPixelErr()  const;
    10599
    106100  Float_t GetMeanFluxPhotonsInnerPixel()     const { return fMeanFluxPhotonsInnerPixel;     }
    107   Float_t GetMeanFluxPhotonsInnerPixelErr()  const { return fMeanFluxPhotonsInnerPixelErr;  }
     101  Float_t GetMeanFluxPhotonsInnerPixelErr()  const;
    108102  Float_t GetMeanFluxPhotonsOuterPixel()     const { return fMeanFluxPhotonsOuterPixel;     }
    109   Float_t GetMeanFluxPhotonsOuterPixelErr()  const { return fMeanFluxPhotonsOuterPixelErr;  }
     103  Float_t GetMeanFluxPhotonsOuterPixelErr()  const;
    110104
    111105  Bool_t IsBlindPixelMethodValid()   const;
    … …  
    138132  Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const;
    139133
    140   Bool_t CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, const MBadPixelsCam &bad); 
     134  Bool_t CalcMeanFluxPhotonsFFactorMethod(const MGeomCam &geom, MBadPixelsCam &bad); 
    141135
    142   void   ApplyPINDiodeCalibration(const MGeomCam &geom, const MBadPixelsCam &bad);
    143   void   ApplyBlindPixelCalibration(const MGeomCam &geom, const MBadPixelsCam &bad);
     136  void   ApplyPINDiodeCalibration(const MGeomCam &geom,
     137                                  const MBadPixelsCam &bad,
     138                                  const MCalibrationChargePINDiode &pindiode);
     139  void   ApplyBlindPixelCalibration(const MGeomCam &geom,
     140                                    const MBadPixelsCam &bad,
     141                                    const MCalibrationChargeBlindPix &blindpix);
    144142  void   ApplyFFactorCalibration(const MGeomCam &geom, const MBadPixelsCam &bad);
    145143

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

    r3496 r3511  
    8383//  since we might be applying it to PMTs in the totally wrong direction.
    8484//
     85//
     86//  Error of all variables are calculated by error-propagation. Note that internally,
     87//  all error variables contain Variances in order to save the CPU-intensive square rooting
    8588//
    8689/////////////////////////////////////////////////////////////////////////////
    … …  
    109112const Float_t MCalibrationChargePix::fgTimeLimit                = 1.5;
    110113const Float_t MCalibrationChargePix::fgTimeErrLimit             = 3.;
     114const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit    = 5.;
    111115// --------------------------------------------------------------------------
    112116//
    … …  
    136140  SetTimeLimit();
    137141  SetTimeErrLimit();
     142  SetPheFFactorMethodLimit();
     143 
    138144}
    139145
    … …  
    156162
    157163  fHiGainMeanCharge                 =  -1.;
    158   fHiGainMeanChargeErr              =  -1.;
     164  fHiGainMeanChargeVar              =  -1.;
    159165  fHiGainSigmaCharge                =  -1.;
    160   fHiGainSigmaChargeErr             =  -1.;
     166  fHiGainSigmaChargeVar             =  -1.;
    161167  fHiGainChargeProb                 =  -1.;
    162168
    163169  fLoGainMeanCharge                 =  -1.;
    164   fLoGainMeanChargeErr              =  -1.;
     170  fLoGainMeanChargeVar              =  -1.;
    165171  fLoGainSigmaCharge                =  -1.;
    166   fLoGainSigmaChargeErr             =  -1.;
     172  fLoGainSigmaChargeVar             =  -1.;
    167173  fLoGainChargeProb                 =  -1.;
    168174
    169175  fRSigmaCharge                     =  -1.;
    170   fRSigmaChargeErr                  =  -1.;
     176  fRSigmaChargeVar                  =  -1.;
    171177 
    172178  fHiGainNumPickup                  =  -1;
    … …  
    177183  fPed                              =  -1.;
    178184  fPedRms                           =  -1.;
    179   fPedErr                           =  -1.;
     185  fPedVar                           =  -1.;
    180186
    181187  fLoGainPedRms                     =  -1.;
    182   fLoGainPedRmsErr                  =  -1.;
     188  fLoGainPedRmsVar                  =  -1.;
    183189
    184190  fTimeFirstHiGain                  =   0 ;
    … …  
    191197
    192198  fPheFFactorMethod                 =  -1.;
    193   fPheFFactorMethodErr              =  -1.;
     199  fPheFFactorMethodVar              =  -1.;
    194200
    195201  fMeanConversionFFactorMethod      =  -1.;
    … …  
    198204  fMeanConversionCombinedMethod     =  -1.;
    199205
    200   fConversionFFactorMethodErr       =  -1.;
    201   fConversionBlindPixelMethodErr    =  -1.;
    202   fConversionPINDiodeMethodErr      =  -1.;
    203   fConversionCombinedMethodErr      =  -1.;
     206  fConversionFFactorMethodVar       =  -1.;
     207  fConversionBlindPixelMethodVar    =  -1.;
     208  fConversionPINDiodeMethodVar      =  -1.;
     209  fConversionCombinedMethodVar      =  -1.;
    204210
    205211  fSigmaConversionFFactorMethod     =  -1.;
    … …  
    231237  fPed       = ped;   
    232238  fPedRms    = pedrms;
    233   fPedErr    = pederr;
     239  fPedVar    = pederr*pederr;
    234240}
    235241 
    … …  
    246252{
    247253    if (IsHiGainSaturation())
    248         fLoGainMeanChargeErr = f;
     254        fLoGainMeanChargeVar = f*f;
    249255    else
    250         fHiGainMeanChargeErr = f;
     256        fHiGainMeanChargeVar = f*f;
    251257
    252258}
    … …  
    264270{
    265271    if (IsHiGainSaturation())
    266         fLoGainSigmaChargeErr = f;
     272        fLoGainSigmaChargeVar = f*f;
    267273    else
    268         fHiGainSigmaChargeErr = f;
    269 
    270 }
    271 
    272 
     274        fHiGainSigmaChargeVar = f*f;
     275
     276}
    273277
    274278// --------------------------------------------------------------------------
    … …  
    279283{
    280284  fMeanConversionFFactorMethod  = c;
    281   fConversionFFactorMethodErr   = err;
     285  fConversionFFactorMethodVar   = err*err;
    282286  fSigmaConversionFFactorMethod = sig;
    283287}
    … …  
    290294{
    291295  fMeanConversionCombinedMethod  = c;
    292   fConversionCombinedMethodErr   = err;
     296  fConversionCombinedMethodVar   = err*err;
    293297  fSigmaConversionCombinedMethod = sig;
    294298}
    … …  
    302306{
    303307  fMeanConversionBlindPixelMethod  = c;
    304   fConversionBlindPixelMethodErr   = err;
     308  fConversionBlindPixelMethodVar   = err*err;
    305309  fSigmaConversionBlindPixelMethod = sig;
    306310}
    … …  
    313317{
    314318  fMeanConversionPINDiodeMethod  = c ;
    315   fConversionPINDiodeMethodErr   = err;
     319  fConversionPINDiodeMethodVar   = err*err;
    316320  fSigmaConversionPINDiodeMethod = sig;
    317321}
    … …  
    408412void MCalibrationChargePix::SetAbsTimeBordersLoGain(const Byte_t f, const Byte_t l)
    409413{
    410 
    411414  fTimeFirstLoGain = f;
    412415  fTimeLastLoGain  = l;
    413  
    414416}
    415417
    … …  
    421423Float_t  MCalibrationChargePix::GetPedRmsErr()  const
    422424{
    423   return IsHiGainSaturation() ? fLoGainPedRmsErr : fPedErr/2.;
     425  return IsHiGainSaturation() ? TMath::Sqrt(fLoGainPedRmsVar) : TMath::Sqrt(fPedVar)/2.;
     426}
     427
     428Float_t  MCalibrationChargePix::GetPedErr()  const
     429{
     430  return  TMath::Sqrt(fPedVar);
    424431}
    425432
    426433Float_t  MCalibrationChargePix::GetMeanCharge()   const
    427434{
    428     return IsHiGainSaturation() ? fLoGainMeanCharge : fHiGainMeanCharge ;
     435    return IsHiGainSaturation() ? GetLoGainMeanCharge() : GetHiGainMeanCharge() ;
    429436}
    430437
    431438Float_t  MCalibrationChargePix::GetMeanChargeErr()   const
    432439{
    433     return IsHiGainSaturation() ? fLoGainMeanChargeErr : fHiGainMeanChargeErr ;
     440    return IsHiGainSaturation() ? GetLoGainMeanChargeErr() : GetHiGainMeanChargeErr() ;
    434441}
    435442
    … …  
    441448Float_t  MCalibrationChargePix::GetSigmaCharge()   const
    442449{
    443     return IsHiGainSaturation() ? fLoGainSigmaCharge : fHiGainSigmaCharge ;
     450    return IsHiGainSaturation() ? GetLoGainSigmaCharge() : GetHiGainSigmaCharge() ;
    444451}
    445452
    446453Float_t  MCalibrationChargePix::GetSigmaChargeErr()   const
    447454{
    448     return IsHiGainSaturation() ? fLoGainSigmaChargeErr : fHiGainSigmaChargeErr ;
     455    return IsHiGainSaturation() ? GetLoGainSigmaChargeErr() : GetHiGainSigmaChargeErr() ;
     456}
     457
     458Float_t MCalibrationChargePix::GetHiGainMeanChargeErr()  const
     459{
     460  return TMath::Sqrt(fHiGainMeanChargeVar);
     461}
     462
     463Float_t MCalibrationChargePix::GetLoGainMeanCharge()  const
     464{
     465  return fLoGainMeanCharge * fConversionHiLo;
     466}
     467
     468Float_t MCalibrationChargePix::GetLoGainMeanChargeErr()  const
     469{
     470
     471  const Float_t chargeRelVar     =  fLoGainMeanChargeVar
     472                                 /( fLoGainMeanCharge * fLoGainMeanCharge );
     473
     474  const Float_t conversionRelVar =  fConversionHiLoVar
     475                                 /( fConversionHiLo   * fConversionHiLo   );
     476
     477  return TMath::Sqrt(chargeRelVar+conversionRelVar) * GetLoGainMeanCharge();
     478}
     479
     480Float_t MCalibrationChargePix::GetLoGainSigmaCharge()  const
     481{
     482  return fLoGainSigmaCharge * fConversionHiLo;
     483}
     484
     485Float_t MCalibrationChargePix::GetLoGainSigmaChargeErr()  const
     486{
     487
     488  const Float_t sigmaRelVar     =  fLoGainSigmaChargeVar
     489                                /( fLoGainSigmaCharge * fLoGainSigmaCharge );
     490
     491  const Float_t conversionRelVar =  fConversionHiLoVar
     492                                 /( fConversionHiLo   * fConversionHiLo   );
     493
     494  return TMath::Sqrt(sigmaRelVar+conversionRelVar) * GetLoGainSigmaCharge();
     495}
     496
     497Float_t MCalibrationChargePix::GetHiGainSigmaChargeErr()  const
     498{
     499  return TMath::Sqrt(fHiGainSigmaChargeVar);
     500}
     501
     502Float_t MCalibrationChargePix::GetRSigmaCharge()  const
     503{
     504  return IsHiGainSaturation() ? fRSigmaCharge*fConversionHiLo : fRSigmaCharge ;
     505}
     506
     507Float_t MCalibrationChargePix::GetRSigmaChargeErr()  const
     508{
     509 if (IsHiGainSaturation())
     510    {
     511      const Float_t rsigmaRelVar     =  fRSigmaChargeVar
     512                                    /( fRSigmaCharge * fRSigmaCharge );
     513      const Float_t conversionRelVar =  fConversionHiLoVar
     514                                     /( fConversionHiLo   * fConversionHiLo   );
     515      return TMath::Sqrt(rsigmaRelVar+conversionRelVar) * GetRSigmaCharge();
     516    }
     517 else
     518   return TMath::Sqrt(fRSigmaChargeVar);
     519
     520}
     521
     522Float_t MCalibrationChargePix::GetConversionHiLoErr()  const
     523{
     524  if (fConversionHiLoVar < 0.)
     525    return -1.;
     526  return TMath::Sqrt(fConversionHiLoVar);
     527}
     528
     529Float_t MCalibrationChargePix::GetPheFFactorMethodErr()  const
     530{
     531  if (fPheFFactorMethodVar < 0.)
     532    return -1.;
     533  return TMath::Sqrt(fPheFFactorMethodVar);
     534}
     535
     536Float_t MCalibrationChargePix::GetConversionCombinedMethodErr()  const
     537{
     538  if (fConversionCombinedMethodVar < 0.)
     539    return -1.;
     540  return TMath::Sqrt(fConversionCombinedMethodVar);
     541}
     542
     543Float_t MCalibrationChargePix::GetConversionPINDiodeMethodErr()  const
     544{
     545  if (fConversionPINDiodeMethodVar < 0.)
     546    return -1.;
     547  return TMath::Sqrt(fConversionPINDiodeMethodVar);
     548}
     549
     550Float_t MCalibrationChargePix::GetConversionBlindPixelMethodErr()  const
     551{
     552  if (fConversionBlindPixelMethodVar < 0.)
     553    return -1.;
     554  return TMath::Sqrt(fConversionBlindPixelMethodVar);
     555}
     556
     557Float_t MCalibrationChargePix::GetConversionFFactorMethodErr()  const
     558{
     559  if (fConversionFFactorMethodVar < 0.)
     560    return -1.;
     561  return TMath::Sqrt(fConversionFFactorMethodVar);
     562}
     563
     564Float_t MCalibrationChargePix::GetTotalFFactorCombinedMethodErr()  const
     565{
     566  if (fTotalFFactorCombinedMethodVar < 0.)
     567    return -1.;
     568  return TMath::Sqrt(fTotalFFactorCombinedMethodVar);
     569}
     570
     571Float_t MCalibrationChargePix::GetTotalFFactorPINDiodeMethodErr()  const
     572{
     573  if (fTotalFFactorPINDiodeMethodVar < 0.)
     574    return -1.;
     575  return TMath::Sqrt(fTotalFFactorPINDiodeMethodVar);
     576}
     577
     578Float_t MCalibrationChargePix::GetTotalFFactorBlindPixelMethodErr()  const
     579{
     580  if (fTotalFFactorBlindPixelMethodVar < 0.)
     581    return -1.;
     582  return TMath::Sqrt(fTotalFFactorBlindPixelMethodVar);
     583}
     584
     585Float_t MCalibrationChargePix::GetTotalFFactorFFactorMethodErr()  const
     586{
     587  if (fTotalFFactorFFactorMethodVar < 0.)
     588    return -1.;
     589  return TMath::Sqrt(fTotalFFactorFFactorMethodVar);
    449590}
    450591
    … …  
    454595}
    455596
    456 
    457597Bool_t MCalibrationChargePix::IsExcluded()     const
    458598{
    … …  
    505645//
    506646// 1) Pixel has a fitted charge greater than fChargeLimit*PedRMS
    507 // 2) Pixel has a fit error greater than fChargeErrLimit
    508 // 3) Pixel has a fitted charge greater its fChargeRelErrLimit times its charge error
     647// 2) Pixel has a fit error greater than fChargeVarLimit
     648// 3) Pixel has a fitted charge greater its fChargeRelVarLimit times its charge error
    509649// 4) Pixel has a charge sigma bigger than its Pedestal RMS
    510650//
    … …  
    520660    }
    521661 
    522   if (GetMeanChargeErr() < fChargeErrLimit)
    523     {
    524       *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
    525             << fChargeErrLimit << " in Pixel  " << fPixId << endl;
     662  const Float_t meanchargevar = IsHiGainSaturation() ? fLoGainMeanChargeVar : fHiGainMeanChargeVar;
     663 
     664  if (meanchargevar < fChargeVarLimit)
     665    {
     666      *fLog << warn << "WARNING: Variance of Fitted Charge is smaller than "
     667            << meanchargevar << " in Pixel  " << fPixId << endl;
    526668      bad->SetChargeErrNotValid();
    527669      bad->SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    528670    }
    529671     
    530    if (GetMeanCharge() < fChargeRelErrLimit*GetMeanChargeErr())
     672   if (GetMeanCharge()*GetMeanCharge() < fChargeRelVarLimit*meanchargevar)
    531673    {
    532674      *fLog << warn << "WARNING: Fitted Charge is smaller than "
    533             << fChargeRelErrLimit << "* its error in Pixel  " << fPixId << endl;
     675            << TMath::Sqrt(fChargeRelVarLimit) << "* its error in Pixel  " << fPixId << endl;
    534676      bad->SetChargeRelErrNotValid();
    535677      bad->SetUnsuitable(MBadPixelsPix::kUnsuitableRun);
    … …  
    577719{
    578720
    579     Float_t pedRmsSquare          = fPedRms * fPedRms;
    580     Float_t pedRmsSquareErrSquare = fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
    581     Float_t pedRmsSquareErr;
     721    Float_t pedRmsSquare      = fPedRms * fPedRms;
     722    Float_t pedRmsSquareVar   = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
    582723   
    583724    //
    … …  
    587728    // We extract the pure NSB contribution:
    588729    //
    589     const Float_t elecRmsSquare          =    fElectronicPedRms    * fElectronicPedRms;
    590     const Float_t elecRmsSquareErrSquare = 4.*fElectronicPedRmsErr * fElectronicPedRmsErr * elecRmsSquare;
     730    const Float_t elecRmsSquare    =    fElectronicPedRms    * fElectronicPedRms;
     731    const Float_t elecRmsSquareVar = 4.*fElectronicPedRmsVar * elecRmsSquare;
    591732   
    592     Float_t nsbSquare             =  pedRmsSquare          - elecRmsSquare;
    593     Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
     733    Float_t nsbSquare             =  pedRmsSquare    - elecRmsSquare;
     734    Float_t nsbSquareRelVar       = (pedRmsSquareVar + elecRmsSquareVar)
    594735                                        / (nsbSquare * nsbSquare) ;
    595736   
    596737    if (nsbSquare < 0.)
    597         nsbSquare = 0.;
     738      nsbSquare = 0.;
    598739   
    599740    //
    … …  
    601742    // add it quadratically to the electronic noise
    602743    //
    603     const Float_t conversionSquare             =    fConversionHiLo    * fConversionHiLo;
    604     const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoErr * fConversionHiLoErr / conversionSquare;
     744    const Float_t conversionSquare        =    fConversionHiLo    * fConversionHiLo;
     745    const Float_t convertedNsbSquare      =    nsbSquare       / conversionSquare;
     746    const Float_t convertedNsbSquareVar   =    nsbSquareRelVar
     747                                            * convertedNsbSquare * convertedNsbSquare;
    605748   
    606     const Float_t convertedNsbSquare          =  nsbSquare             / conversionSquare;
    607     const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
    608                                                 * convertedNsbSquare * convertedNsbSquare;
     749    pedRmsSquare     = convertedNsbSquare    + elecRmsSquare;
     750    pedRmsSquareVar  = convertedNsbSquareVar + elecRmsSquareVar;
    609751   
    610     pedRmsSquare     = convertedNsbSquare           + elecRmsSquare;
    611     pedRmsSquareErr  = TMath::Sqrt( convertedNsbSquareErrSquare  + elecRmsSquareErrSquare );
    612    
    613     //
    614     // Correct also for the conversion to Hi-Gain:
    615     //
    616     fLoGainPedRms    = TMath::Sqrt(pedRmsSquare)               * fConversionHiLo;
    617     fLoGainPedRmsErr = 0.5 * pedRmsSquareErr /  fLoGainPedRms  * fConversionHiLo;
    618 }
    619 
     752    fLoGainPedRms    = TMath::Sqrt(pedRmsSquare);
     753    fLoGainPedRmsVar = 0.25 * pedRmsSquareVar /  pedRmsSquare;
     754
     755}
     756
     757//
     758//
     759//
    620760Bool_t MCalibrationChargePix::CalcReducedSigma()
    621761{
    622762
    623   const Float_t sigmaSquare               =    GetSigmaCharge()   * GetSigmaCharge();
    624   const Float_t sigmaSquareErrSquare      = 4.*GetSigmaChargeErr()* GetSigmaChargeErr() * sigmaSquare;
    625 
    626   Float_t pedRmsSquare;         
    627   Float_t pedRmsSquareErrSquare;
     763  const Float_t sigmacharge    = IsHiGainSaturation() ? fLoGainSigmaCharge    : fHiGainSigmaCharge   ;
     764  const Float_t sigmachargevar = IsHiGainSaturation() ? fLoGainSigmaChargeVar : fHiGainSigmaChargeVar;
     765
     766  const Float_t sigmaSquare     =      sigmacharge     * sigmacharge;
     767  const Float_t sigmaSquareVar  = 4.*  sigmachargevar  * sigmaSquare;
     768
     769  Float_t pedRmsSquare ;         
     770  Float_t pedRmsSquareVar;
    628771 
    629772  if (IsHiGainSaturation())
    630   {
    631       pedRmsSquare             =      fLoGainPedRms    * fLoGainPedRms;
    632       pedRmsSquareErrSquare    =  4.* fLoGainPedRmsErr * fLoGainPedRmsErr * pedRmsSquare;
     773    {
     774      pedRmsSquare     =      fLoGainPedRms    * fLoGainPedRms;
     775      pedRmsSquareVar  =  4.* fLoGainPedRmsVar * pedRmsSquare;
    633776    }
    634777  else
    635778  {
    636 
    637       pedRmsSquare           =       fPedRms * fPedRms;                                         
    638       pedRmsSquareErrSquare  =       fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
     779      pedRmsSquare       =  fPedRms * fPedRms;                                         
     780      pedRmsSquareVar    =  fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
    639781  }
    640782  //
    … …  
    642784  //
    643785  const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
    644 
    645786  if (rsigmachargesquare <= 0.)
    646787    {
    … …  
    651792    }
    652793 
     794
    653795  fRSigmaCharge    = TMath::Sqrt(rsigmachargesquare);
    654   fRSigmaChargeErr = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare) / 2. / fRSigmaCharge;
     796  fRSigmaChargeVar = 0.25 * (sigmaSquareVar + pedRmsSquareVar) / rsigmachargesquare;
    655797
    656798  return kTRUE;
    … …  
    670812    }
    671813 
     814  const Float_t charge    = IsHiGainSaturation() ? fLoGainMeanCharge    : fHiGainMeanCharge   ;
     815  const Float_t chargevar = IsHiGainSaturation() ? fLoGainMeanChargeVar : fHiGainMeanChargeVar;
     816
    672817  //
    673818  // Square all variables in order to avoid applications of square root
    … …  
    675820  // First the relative error squares
    676821  //
    677   const Float_t chargeSquare              =     GetMeanCharge()    * GetMeanCharge();
    678   const Float_t chargeSquareRelErrSquare  = 4.* GetMeanChargeErr() * GetMeanChargeErr() / chargeSquare;
    679 
    680   const Float_t ffactorsquare             =    gkFFactor    * gkFFactor;
    681   const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
    682 
    683 
    684   const Float_t rsigmaSquare              =     fRSigmaCharge    * fRSigmaCharge;
    685   const Float_t rsigmaSquareRelErrSquare  = 4.* fRSigmaChargeErr * fRSigmaChargeErr / rsigmaSquare;
    686 
     822  const Float_t chargeSquare        =     charge * charge;
     823  const Float_t chargeSquareRelVar  = 4.* chargevar/ chargeSquare;
     824
     825  const Float_t ffactorsquare       =    gkFFactor    * gkFFactor;
     826  const Float_t ffactorsquareRelVar = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
     827
     828  const Float_t rsigmaSquare        =     fRSigmaCharge    * fRSigmaCharge;
     829  const Float_t rsigmaSquareRelVar  = 4.* fRSigmaChargeVar / rsigmaSquare;
    687830
    688831  //
    … …  
    692835  fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmaSquare;
    693836
    694   if (fPheFFactorMethod < 0.)
     837  if (fPheFFactorMethod < fPheFFactorMethodLimit)
    695838    {
    696839      SetFFactorMethodValid(kFALSE);
    … …  
    701844  // Calculate the Error of Nphe
    702845  //
    703   const Float_t pheFFactorRelErrSquare =  ffactorsquareRelErrSquare
    704                                          + chargeSquareRelErrSquare
    705                                          + rsigmaSquareRelErrSquare;
    706   fPheFFactorMethodErr                 =  TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
     846  fPheFFactorMethodVar =  (ffactorsquareRelVar + chargeSquareRelVar + rsigmaSquareRelVar)
     847                         * fPheFFactorMethod * fPheFFactorMethod;
    707848
    708849  SetFFactorMethodValid(kTRUE);
    709 
    710850  return kTRUE;
    711851}
    … …  
    715855{
    716856 
    717   const Float_t chargeRelErrSquare     =    fLoGainMeanChargeErr * fLoGainMeanChargeErr
    718                                          /( fLoGainMeanCharge    * fLoGainMeanCharge  );
    719   const Float_t sigmaRelErrSquare      =    fLoGainSigmaChargeErr * fLoGainSigmaChargeErr
    720                                          /( fLoGainSigmaCharge    * fLoGainSigmaCharge );
    721   const Float_t conversionRelErrSquare =    fConversionHiLoErr * fConversionHiLoErr
    722                                          /( fConversionHiLo    * fConversionHiLo  );
    723  
    724   fLoGainMeanCharge         *= fConversionHiLo;
    725   fLoGainMeanChargeErr       = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fLoGainMeanCharge;
    726  
    727   fLoGainSigmaCharge        *= fConversionHiLo;
    728   fLoGainSigmaChargeErr      =  TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fLoGainSigmaCharge;
    729  
    730   fElectronicPedRms     = gkElectronicPedRms    * TMath::Sqrt(fNumLoGainSamples);
    731   fElectronicPedRmsErr  = gkElectronicPedRmsErr * TMath::Sqrt(fNumLoGainSamples);
     857  fElectronicPedRms       = gkElectronicPedRms    * TMath::Sqrt(fNumLoGainSamples);
     858  fElectronicPedRmsVar    = gkElectronicPedRmsErr * gkElectronicPedRmsErr * fNumLoGainSamples;
    732859 
    733860  CalcLoGainPed();
    734 
    735 }
    736 
    737 
     861}
     862
     863
Note: See TracChangeset for help on using the changeset viewer.