Changeset 5812 for trunk/MagicSoft/Mars/mcalib

Timestamp:

01/12/05 20:41:09 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

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

trunk/MagicSoft/Mars/mcalib/MCalibrationIntensityRelTimeCam.cc

@@ -47 +47 @@
 #include "MGeomPix.h"
 
+#include "MHCamera.h"
+
 #include "MLogManip.h"
 
 #include <TOrdCollection.h>
 #include <TGraphErrors.h>
+#include <TH1D.h>
+#include <TF1.h>
 
 ClassImp(MCalibrationIntensityRelTimeCam);
@@ -177 +181 @@
   const Float_t sqrt2   = 1.414;
   const Float_t fadc2ns = 3.333;
-  
+  const Float_t norm    = fadc2ns / sqrt2;
+
   TArrayF res(size);
   TArrayF reserr(size);
@@ -185 +190 @@
   Int_t cnt = 0;
 
+  TH1D *h = 0;
+
   for (Int_t i=0;i<GetSize();i++)
     {
@@ -193 +200 @@
       MCalibrationChargeCam *cam = (MCalibrationChargeCam*)chargecam.GetCam(i);
 
-      if (col != MCalibrationCam::kNONE)
-        if (relcam->GetPulserColor() != col)
-          continue;
+      if (relcam->GetPulserColor() != col)
+        continue;
 
       const MCalibrationChargePix &apix = (MCalibrationChargePix&)cam->GetAverageArea(aidx);
@@ -208 +214 @@
       Double_t var    = 0.;
       Int_t    num    = 0;
-     //
+
+      MHCamera camres(geom,"CamRes","Time Resolution;Time Resolution [ns];channels");
+      //
       // Get the area calibration pix from the calibration cam
       //
-      for (Int_t i=0; i<cam->GetSize(); i++)
+      for (Int_t j=0; j<cam->GetSize(); j++)
         {
-          const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[i];
-          const MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[i];
+          const MCalibrationChargePix &pix = (MCalibrationChargePix&)(*cam)[j];
+          const MCalibrationRelTimePix &relpix = (MCalibrationRelTimePix&)(*relcam)[j];
           //
           // Don't use bad pixels
@@ -222 +230 @@
           //
           //
-          if (aidx != geom[i].GetAidx())
+          if (aidx != geom[j].GetAidx())
             continue;
           
-          resol  += relpix.GetTimePrecision();
-          resol2 += relpix.GetTimePrecision()*relpix.GetTimePrecision();
+          const Float_t pres = relpix.GetTimePrecision();
+
+          resol  += pres;
+          resol2 += pres;
           num++;
+          
+          camres.Fill(j,pres);
+          camres.SetUsed(j);
         }
       
@@ -235 +248 @@
           var  = (resol2 - resol*resol*num) / (num-1);
 
-          res[cnt] = resol * fadc2ns / sqrt2;
+          res[cnt] = resol;
           if (var > 0.)
-            reserr[cnt] = TMath::Sqrt(var)/ sqrt2 * fadc2ns;
+            reserr[cnt] = TMath::Sqrt(var);
           else
             reserr[cnt] = 0.;
+
+          //
+          // Make also a Gauss-fit to the distributions. The RMS can be determined by 
+          // outlier, thus we look at the sigma and the RMS and take the smaller one, afterwards.
+          // 
+          h = camres.ProjectionS(TArrayI(),TArrayI(1,&aidx),"_py",750);
+          h->SetDirectory(NULL);
+          h->Fit("gaus","QL");
+          TF1 *fit = h->GetFunction("gaus");
+
+          Float_t ci2   = fit->GetChisquare();
+          Float_t sigma = fit->GetParameter(2);
+
+          if (ci2 > 500. || sigma > reserr[cnt])
+            {
+              h->Fit("gaus","QLM");
+              fit = h->GetFunction("gaus");
+
+              ci2   = fit->GetChisquare();
+              sigma = fit->GetParameter(2);
+            }
+          
+          const Float_t mean  = fit->GetParameter(1);
+          const Float_t ndf   = fit->GetNDF();
+
+          
+          *fLog << inf << "Mean number photo-electrons: " << sig[cnt] << endl;
+          *fLog << inf << "Time Resolution area idx: " << aidx << " Results: " << endl;
+          *fLog << inf << "Mean: " << Form("%4.3f",mean) 
+                << "+-" << Form("%4.3f",fit->GetParError(1))
+                << "  Sigma: " << Form("%4.3f",sigma) << "+-" << Form("%4.3f",fit->GetParError(2)) 
+                << "  Chisquare: " << Form("%4.3f",fit->GetChisquare()) << "  NDF  : " << ndf << endl;          
+
+          delete h;
+          gROOT->GetListOfFunctions()->Remove(fit);
+
+          if (sigma < reserr[cnt] && ndf > 2)
+            {
+              res   [cnt] = mean;
+              reserr[cnt] = sigma;
+            }
+
+          res[cnt]    *= norm;
+          reserr[cnt] *= norm;
         }
       else
@@ -247 +304 @@
         }
       cnt++;
-    }
+   }
   
   TGraphErrors *gr = new TGraphErrors(size,