Back to index

CSidcCalibrator.C

 
//----------------------------------------------------------------------------- 
//  $Header: /asis/offline/ceres/cool/project/RCS/CSidcCalibrator.C,v 3.0 1996/10/02 09:40:10 voigt Exp $ 
// 
//  COOL Program Library   
//  Copyright (C) CERES collaboration, 1996 
// 
//  Implementation of class CSidcCalibrator 
// 
//----------------------------------------------------------------------------- 
#include <math.h>   
#include "CSidcCalibrator.h"  
#include "CSidc.h"  
#include "CMinuitWrap.h" 
#include "CMCDigiHit.h" 
 
void edgefit(int*, double*, double*, double*, int*, void*); 
 
CSidcCalibrator::CSidcCalibrator(CSidc* newSidc)  
{ 
   // 
   // hold pointer to the SiDC 
   // 
   sidc = newSidc; 
    
   // 
   // initialize 'histogram' 
   // 
   minTimeBin = sidc->getSetup()->getMaxTimeBinReference() - 3.; 
   binWidth   = 0.1; 
   resetVector(); 
    
   // 
   // initialize MINUIT stuff 
   // 
   strcpy(names[0], "time     "); 
   strcpy(names[1], "slope    "); 
   strcpy(names[2], "total    "); 
   strcpy(names[3], "noise    "); 
   nParams = MaxParams; 
   fmin    = 0; 
   flag    = 0; 
   ierr    = 0; 
   istat   = 0; 
   for (int i=0; i<MaxParams; i++) { 
      grad  [i] = 0; 
      result[i] = 0; 
      val   [i] = 0; 
      step  [i] = 0; 
   } 
   dummy  = 0; 
   fedm   = 0; 
   errdef = 0; 
   idummy = 0; 
} 
 
CSidcCalibrator::~CSidcCalibrator()  
{ 
} 
 
void CSidcCalibrator::update() 
{ 
   CSidcCellCoord hitCenter; 
   double index; 
   for (int i=0; i<sidc->getHits().length(); i++) { 
      sidc->transform(hitCenter, sidc->getHits()(i)->getCenter()); 
      index = (hitCenter.getTimeBin() - minTimeBin) / binWidth; 
      if (index >= 0 && index < timeBinsInVector) 
	 timeBinVector(int(index)) += 1.0; 
   } 
   numberOfEvents++; 
} 
 
CBoolean CSidcCalibrator::doCalibration() 
{ 
   static CBoolean init = False; 
    
   if (!init) { 
      mninit(5, 6, 42);                             // initialize MINUIT 
      init = True; 
   } 
    
   val[0] = 35.0;      // position of half amplitude in timeBinVector 
   val[1] = 200.0;     // slope 
   val[2] = 300.0;     // total 
   val[3] = 20.0;      // noise 
    
   step[0] = 20.0; 
   step[1] = 50.0; 
   step[2] = 50.0;              
   step[3] = 10.0;         
    
   // 
   // initlialize fit parameters 
   // 
 
   double argList; 
   argList = -1; 
   mnexcm(&edgefit, "SET PRInt"   , &argList, 1, &ierr, (void *) 0);         // no printout 
   argList = 0; 
   mnexcm(&edgefit, "SET NOWarn"  , &argList, 0, &ierr, (void *) 0); 
   mnexcm(&edgefit, "CLEar"       , &argList, 0, &ierr, (void *) 0);    
   for (int i=0; i<MaxParams; i++) { 
      mnparm(i+1, names[i], val[i], step[i], 0, 0, &ierr); 
      if (ierr != 0) { 
         cerr << "CCSidcCalibrator::doCalibration():\n"; 
         cerr << "\tERROR\n"; 
         cerr << "\terror from mnparm for i= " << i << " " << ierr << endl; 
	 return False; 
      } 
   }    
        
   // 
   // Perform edgefit 
   // 
   mnexcm (&edgefit, "MINImize", 0, 0, &ierr, (void *) 0); 
   if (ierr != 0) { 
      cerr << "CCSidcCalibrator::doCalibration():\n"; 
      cerr << "\tERROR\n"; 
      cerr << "\terror from mnexcm: " << ierr << endl; 
   } 
    
   // 
   // Check fit 
   // 
   mnstat (&fmin, &fedm, &errdef, &idummy, &idummy, &istat); 
   if (istat < 2) { 
      cerr << "CCSidcCalibrator::doCalibration():\n"; 
      cerr << "\tERROR\n"; 
      cerr << "\tEdgefit failed: ierr = " << ierr << endl; 
      return False; 
   } 
    
   // 
   // Store fit results 
   //idummyidummy 
   mnpout(1, " ", &val[0], errdef, dummy, dummy, idummy); 
   mnpout(2, " ", &val[1], errdef, dummy, dummy, idummy); 
   mnpout(3, " ", &val[2], errdef, dummy, dummy, idummy);     
   mnpout(4, " ", &val[3], errdef, dummy, dummy, idummy);     
    
   return True; 
} 
 
void edgefit(int*, double [4], double *f, double qs[4], int*, void*) 
{ 
   *f = 0; 
    
   double functionValue;    
   for (int i=0; i<timeBinsInVector; i++) { 
      functionValue = qs[3] + qs[2] / (exp((double(i)-qs[0])/qs[1]) + 1); 
      *f += (timeBinVector(i)-functionValue) * (timeBinVector(i)-functionValue); 
   }       
} 
 
void CSidcCalibrator::resetVector() 
{ 
   numberOfEvents = 0; 
   timeBinVector = 0.0; 
} 
 
void CSidcCalibrator::printCalibration() 
{ 
   cout << "position of half amplitude : " << val[0] << endl; 
   cout << "slope : " << val[1] << endl; 
   cout << "total : " << val[2] << endl; 
   cout << "noise : " << val[3] << endl; 
} 

Back to index