Changeset 3065 for trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc

Timestamp:

02/09/04 10:36:04 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc (modified) (8 diffs)

trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc

@@ -262 +262 @@
     {
       
-      if (pix->IsChargeFitValid() && !pix->IsExcluded()) 
+      if (pix->IsChargeValid() && !pix->IsExcluded()) 
         {
 
@@ -284 +284 @@
       {
         
-        if (!pix->IsChargeFitValid() && !pix->IsExcluded())
+        if (!pix->IsChargeValid() && !pix->IsExcluded())
           {
 
@@ -383 +383 @@
   Float_t arearatio = cam.GetPixRatio(idx);
 
-  if (arearatio == 0)
+ if (arearatio == 0)
     return kFALSE;
 
@@ -587 +587 @@
       if (!(*this)[idx].IsFitted())
         return kFALSE;
-      if (!(*this)[idx].IsChargeFitValid())
+      if (!(*this)[idx].IsChargeValid())
         val = 1;
       else
@@ -702 +702 @@
   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 +787 @@
   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 +845 @@
   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 +912 @@
   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;