Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3064) +++ trunk/MagicSoft/Mars/Changelog (revision 3065) @@ -4,4 +4,26 @@ -*-*- END OF LINE -*-*- + 2004/02/09: Markus Gaug + + * MMcCalibrationUpdate.cc + * MCalibrate.cc + * MCalibrationPix.cc + * MCalibrationPix.h + - replace MCalibrationPix->IsChargeFitValid() by IsChargeValid() + because we can calibrate with the histogram mean and RMS even + if the fit is not OK. Thus, there is a question: IsFitted() and + IsChargeValid(), i.e. does the charge make any sense? + + * MCalibrationCam.cc + MCalibrationConfig.h + MCalibrationCalc.cc + - found serious bug in calculation of conversion factors for the + outer pixels: + GetMeanConversionBlindPixelMethod and + GetMeanConversionPINDiodeMethod + did not correct for the area of the outer pixel, this did only + MCalibrationCam::GetPixelContent which is not used by MCalibrate + + 2004/02/06: Markus Gaug Index: trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc (revision 3064) +++ trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc (revision 3065) @@ -222,5 +222,5 @@ calpix.SetBlindPixelMethodValid(); - calpix.SetChargeFitValid(); + calpix.SetChargeValid(); calpix.SetConversionHiLo(fConversionHiLo); Index: trunk/MagicSoft/Mars/mcalib/MCalibrate.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrate.cc (revision 3065) @@ -169,5 +169,5 @@ MCalibrationPix &pix = (*fCalibrations)[pixidx]; - if (!pix.IsChargeFitValid()) + if (!pix.IsChargeValid()) continue; Index: trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3065) @@ -279,4 +279,7 @@ return kFALSE; } + + + fCalibrations->SetGeomCam(cam); fNumHiGainSamples = fSignals->GetNumUsedHiGainFADCSlices(); Index: trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc (revision 3065) @@ -262,5 +262,5 @@ { - if (pix->IsChargeFitValid() && !pix->IsExcluded()) + if (pix->IsChargeValid() && !pix->IsExcluded()) { @@ -284,5 +284,5 @@ { - if (!pix->IsChargeFitValid() && !pix->IsExcluded()) + if (!pix->IsChargeValid() && !pix->IsExcluded()) { @@ -383,5 +383,5 @@ Float_t arearatio = cam.GetPixRatio(idx); - if (arearatio == 0) + if (arearatio == 0) return kFALSE; @@ -587,5 +587,5 @@ if (!(*this)[idx].IsFitted()) return kFALSE; - if (!(*this)[idx].IsChargeFitValid()) + if (!(*this)[idx].IsChargeValid()) val = 1; else @@ -702,37 +702,84 @@ const Float_t mean = fBlindPixel->GetLambda(); const Float_t merr = fBlindPixel->GetErrLambda(); + + // + // Start calculation of number of photons + // + // The blind pixel has exactly one cm^2 area (with negligible error), + // thus we omi division by 1. + // + fMeanPhotInsidePlexiglass = mean * gkCalibrationInnerPixelArea; + + // Start calculation of number of photons relative Variance (!!) + fMeanPhotErrInsidePlexiglass = merr*merr/mean/mean; + fMeanPhotErrInsidePlexiglass += gkCalibrationInnerPixelAreaError*gkCalibrationInnerPixelAreaError + / gkCalibrationInnerPixelArea/gkCalibrationInnerPixelArea; switch (fColor) { case kECGreen: - fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEGreen) // real photons - *TMath::Power(10,gkCalibrationBlindPixelAttGreen) // correct for absorption - * gkCalibrationInnerPixelArea; // correct for area - - + fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEGreen; + fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenError*gkCalibrationBlindPixelQEGreenError + / gkCalibrationBlindPixelQEGreen / gkCalibrationBlindPixelQEGreen; + + fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption + // attenuation has negligible error break; case kECBlue: - fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEBlue ) - *TMath::Power(10,gkCalibrationBlindPixelAttBlue) - * gkCalibrationInnerPixelArea; + fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEBlue; + fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueError*gkCalibrationBlindPixelQEBlueError + / gkCalibrationBlindPixelQEBlue / gkCalibrationBlindPixelQEBlue; + + fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption + // attenuation has negligible error break; case kECUV: - fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEUV ) - *TMath::Power(10,gkCalibrationBlindPixelAttUV) - * gkCalibrationInnerPixelArea; + fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEUV; + fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEUVError*gkCalibrationBlindPixelQEUVError + / gkCalibrationBlindPixelQEUV / gkCalibrationBlindPixelQEUV; + + fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption + // attenuation has negligible error break; case kECCT1: default: - fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQECT1 ) - *TMath::Power(10,gkCalibrationBlindPixelAttCT1) - * gkCalibrationInnerPixelArea; - break; - } - - SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable); + fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQECT1; + fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Error*gkCalibrationBlindPixelQECT1Error + / gkCalibrationBlindPixelQECT1 / gkCalibrationBlindPixelQECT1; + + fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption + // attenuation has negligible error + break; + } *fLog << inf << endl; - *fLog << inf << "Mean number of Photons for an Inner Pixel (inside Plexiglass): " + *fLog << inf << " Mean number of Photons for an Inner Pixel (inside Plexiglass): " << fMeanPhotInsidePlexiglass << endl; + + if (fMeanPhotInsidePlexiglass > 0.) + SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable); + else + { + CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable); + return kFALSE; + } + + if (fMeanPhotErrInsidePlexiglass < 0.) + { + *fLog << warn << " Relative Variance on number of Photons for an Inner Pixel (inside Plexiglass): " + << fMeanPhotErrInsidePlexiglass << endl; + CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable); + return kFALSE; + } + + // Finish calculation of errors -> convert from relative variance to absolute error + fMeanPhotErrInsidePlexiglass = TMath::Sqrt(fMeanPhotErrInsidePlexiglass); + fMeanPhotErrInsidePlexiglass *= fMeanPhotInsidePlexiglass; + + *fLog << inf << " Error on number of Photons for an Inner Pixel (inside Plexiglass): " + << fMeanPhotErrInsidePlexiglass << endl; + *fLog << inf << endl; + + TIter Next(fPixels); @@ -740,10 +787,44 @@ while ((pix=(MCalibrationPix*)Next())) { - if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.)) + + if(pix->IsChargeValid()) { - Float_t conversion = fMeanPhotInsidePlexiglass/pix->GetCharge(); - Float_t conversionerr = 0.; - Float_t conversionsigma = 0.; + const Float_t charge = pix->GetCharge(); + const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId()); + const Float_t chargeerr = pix->GetErrCharge(); + + Float_t nphot; + Float_t nphoterr; + + if (ratio == 1.) + { + nphot = fMeanPhotInsidePlexiglass; + nphoterr = fMeanPhotErrInsidePlexiglass; + } + else + { + nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea + *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError + *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError; + nphoterr = TMath::Sqrt(nphoterr); + } + + + + const Float_t conversion = nphot/charge; + Float_t conversionerr; + + conversionerr = nphoterr/charge + * nphoterr/charge ; + conversionerr += chargeerr/charge + * chargeerr/charge + * conversion*conversion; + conversionerr = TMath::Sqrt(conversionerr); + + const Float_t conversionsigma = 0.; + pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); @@ -764,31 +845,65 @@ const Float_t mean = fPINDiode->GetCharge(); const Float_t merr = fPINDiode->GetErrCharge(); + + // Start calculation of number of photons + fMeanPhotOutsidePlexiglass = mean * gkCalibrationInnerPixelvsPINDiodeArea; + + // Start calculation of number of photons relative Variance (!!) + fMeanPhotErrOutsidePlexiglass = merr*merr/mean/mean; + fMeanPhotErrOutsidePlexiglass += gkCalibrationInnerPixelvsPINDiodeAreaError*gkCalibrationInnerPixelvsPINDiodeAreaError + / gkCalibrationInnerPixelvsPINDiodeArea/gkCalibrationInnerPixelvsPINDiodeArea; switch (fColor) { case kECGreen: - fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEGreen) // real photons - * gkCalibrationInnerPixelvsPINDiodeArea; // correct for area + fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEGreen; + fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEGreenError*gkCalibrationPINDiodeQEGreenError + / gkCalibrationPINDiodeQEGreen/gkCalibrationPINDiodeQEGreen; break; case kECBlue: - fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEBlue ) - * gkCalibrationInnerPixelvsPINDiodeArea; - break; + fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEBlue; + fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEBlueError*gkCalibrationPINDiodeQEBlueError + / gkCalibrationPINDiodeQEBlue/gkCalibrationPINDiodeQEBlue; + break; case kECUV: - fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEUV ) - * gkCalibrationInnerPixelvsPINDiodeArea; + fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEUV; + fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEUVError*gkCalibrationPINDiodeQEUVError + / gkCalibrationPINDiodeQEUV/gkCalibrationPINDiodeQEUV; break; case kECCT1: default: - fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQECT1 ) - * gkCalibrationInnerPixelvsPINDiodeArea; - break; - } - - SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); + fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQECT1; + fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQECT1Error*gkCalibrationPINDiodeQECT1Error + / gkCalibrationPINDiodeQECT1/gkCalibrationPINDiodeQECT1; + break; + } + *fLog << inf << endl; - *fLog << inf << mean << " Mean number of Photons for an Inner Pixel (outside Plexiglass): " + *fLog << inf << " Mean number of Photons for an Inner Pixel (outside Plexiglass): " << fMeanPhotOutsidePlexiglass << endl; + + if (fMeanPhotOutsidePlexiglass > 0.) + SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); + else + { + CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); + return kFALSE; + } + + if (fMeanPhotErrOutsidePlexiglass < 0.) + { + *fLog << warn << "Relative Variance on number of Photons for an Inner Pixel (outside Plexiglass): " + << fMeanPhotErrOutsidePlexiglass << endl; + CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable); + return kFALSE; + } + + // Finish calculation of errors -> convert from relative variance to absolute error + fMeanPhotErrOutsidePlexiglass = TMath::Sqrt(fMeanPhotErrOutsidePlexiglass); + fMeanPhotErrOutsidePlexiglass *= fMeanPhotOutsidePlexiglass; + + *fLog << inf << " Error on number of Photons for an Inner Pixel (outside Plexiglass): " + << fMeanPhotErrOutsidePlexiglass << endl; *fLog << inf << endl; @@ -797,7 +912,49 @@ while ((pix=(MCalibrationPix*)Next())) { - - if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.)) - pix->SetConversionPINDiodeMethod(fMeanPhotOutsidePlexiglass/pix->GetCharge(), 0., 0.); + + if (pix->IsChargeValid()) + { + + const Float_t charge = pix->GetCharge(); + const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId()); + const Float_t chargeerr = pix->GetErrCharge(); + + Float_t nphot; + Float_t nphoterr; + + if (ratio == 1.) + { + nphot = fMeanPhotInsidePlexiglass; + nphoterr = fMeanPhotErrInsidePlexiglass; + } + else + { + nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea + *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea; + nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError + *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError; + nphoterr = TMath::Sqrt(nphoterr); + } + + const Float_t conversion = nphot/charge; + Float_t conversionerr; + + conversionerr = nphoterr/charge + * nphoterr/charge ; + conversionerr += chargeerr/charge + * chargeerr/charge + * conversion*conversion; + if (conversionerr > 0.) + conversionerr = TMath::Sqrt(conversionerr); + + const Float_t conversionsigma = 0.; + + pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma); + + if (conversionerr/conversion < 0.1) + pix->SetBlindPixelMethodValid(); + + } } return kTRUE; Index: trunk/MagicSoft/Mars/mcalib/MCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationCam.h (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationCam.h (revision 3065) @@ -30,4 +30,6 @@ MCalibrationPINDiode *fPINDiode; //-> Pointer to the PIN Diode with fit results + MGeomCam *fGeomCam; //! Need geom cam to know which pixel in inner or outer + TH1D* fOffsets; //! TH1D* fSlopes; //! @@ -69,4 +71,6 @@ void SetBlindPixelMethodValid(const Bool_t b = kTRUE); void SetPINDiodeMethodValid(const Bool_t b = kTRUE); + + void SetGeomCam(MGeomCam *geom) { fGeomCam = geom; } // Getters Index: trunk/MagicSoft/Mars/mcalib/MCalibrationConfig.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationConfig.h (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationConfig.h (revision 3065) @@ -23,4 +23,10 @@ const Float_t gkCalibrationBlindPixelQEUV = 0.247; const Float_t gkCalibrationBlindPixelQECT1 = 0.247; + +// Average QE Error of Blind Pixel (three colours) +const Float_t gkCalibrationBlindPixelQEGreenError = 0.015; +const Float_t gkCalibrationBlindPixelQEBlueError = 0.02; +const Float_t gkCalibrationBlindPixelQEUVError = 0.02; +const Float_t gkCalibrationBlindPixelQECT1Error = 0.02; // Attenuation factor Blind Pixel (three colours) @@ -82,4 +88,11 @@ const Float_t gkCalibrationPINDiodeQECT1 = -1.0; +// +// Average QE of the PIN Diode +// +const Float_t gkCalibrationPINDiodeQEGreenError = -1.0; +const Float_t gkCalibrationPINDiodeQEBlueError = -1.0; +const Float_t gkCalibrationPINDiodeQEUVError = -1.0; +const Float_t gkCalibrationPINDiodeQECT1Error = -1.0; #endif /* MARS_MCalibrationConfig */ Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3065) @@ -107,5 +107,5 @@ CLRBIT(fFlags, kExcluded); CLRBIT(fFlags, kExcludeQualityCheck); - CLRBIT(fFlags, kChargeFitValid); + CLRBIT(fFlags, kChargeValid); CLRBIT(fFlags, kFitted); CLRBIT(fFlags, kOscillating); @@ -258,7 +258,7 @@ // Set the Excluded Bit from outside // -void MCalibrationPix::SetChargeFitValid(Bool_t b ) -{ - b ? SETBIT(fFlags, kChargeFitValid) : CLRBIT(fFlags, kChargeFitValid); +void MCalibrationPix::SetChargeValid(Bool_t b ) +{ + b ? SETBIT(fFlags, kChargeValid) : CLRBIT(fFlags, kChargeValid); } @@ -443,7 +443,7 @@ } -Bool_t MCalibrationPix::IsChargeFitValid() const -{ - return TESTBIT(fFlags, kChargeFitValid); +Bool_t MCalibrationPix::IsChargeValid() const +{ + return TESTBIT(fFlags, kChargeValid); } @@ -584,9 +584,9 @@ ApplyLoGainConversion(); - if (CheckChargeFitValidity()) - SETBIT(fFlags,kChargeFitValid); + if (CheckChargeValidity()) + SETBIT(fFlags,kChargeValid); else { - CLRBIT(fFlags,kChargeFitValid); + CLRBIT(fFlags,kChargeValid); return kFALSE; } @@ -751,5 +751,5 @@ // 5) Pixel has a charge sigma bigger than its Pedestal RMS // -Bool_t MCalibrationPix::CheckChargeFitValidity() +Bool_t MCalibrationPix::CheckChargeValidity() { Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h (revision 3064) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationPix.h (revision 3065) @@ -78,5 +78,5 @@ enum { kHiGainSaturation, kExcluded, kExcludeQualityCheck, - kChargeFitValid, kTimeFitValid, + kChargeValid, kTimeFitValid, kFitted, kOscillating, kBlindPixelMethodValid, kFFactorMethodValid, kPINDiodeMethodValid }; @@ -84,5 +84,5 @@ MHCalibrationPixel *fHist; // Pointer to the histograms performing the fits, etc. - Bool_t CheckChargeFitValidity(); + Bool_t CheckChargeValidity(); Bool_t CheckTimeFitValidity(); Bool_t CalcFFactorMethod(); @@ -148,5 +148,5 @@ Bool_t IsExcludeQualityCheck() const; Bool_t IsHiGainSaturation() const; - Bool_t IsChargeFitValid() const; + Bool_t IsChargeValid() const; Bool_t IsFitted() const; Bool_t IsOscillating(); @@ -169,5 +169,5 @@ void SetExcluded(Bool_t b = kTRUE); void SetExcludeQualityCheck(Bool_t b = kTRUE); - void SetChargeFitValid(Bool_t b = kTRUE); + void SetChargeValid(Bool_t b = kTRUE); void SetFitted(Bool_t b = kTRUE); void SetOscillating(Bool_t b = kTRUE);