Back to index

CCoGPulse.C

 
//----------------------------------------------------------------------------- 
//  $Header: CCoGPulse.C,v 1.6 97/04/21 11:07:05 messer Exp $ 
// 
//  COOL Program Library   
//  Copyright (C) CERES collaboration, 1996 
// 
//  Implementation of class ... 
// 
//----------------------------------------------------------------------------- 
#include <math.h> 
#include "CCoGPulse.h"  
#include "CSidc.h"  
#include "CMCServer.h" 
 
CCoGPulse::CCoGPulse()  
{ 
} 
 
CCoGPulse::CCoGPulse(int newAnode, int newTStart) : CPulse(newAnode, newTStart) 
{ 
   sumAmp = 0; 
} 
 
CCoGPulse::~CCoGPulse()  
{ 
} 
 
CBoolean CCoGPulse::calculateCoG(CSidc& sidc) 
{ 
   // 
   // For pulses with maxAmp not in saturation (threshold value from sidc setup) 
   // indeed only cog is calculated. If, however, there are one or more 
   // tbins in saturation, we do a gaussian regression a la Jan. 
   // 
   CBoolean gaussRegressionWasDone = False; 
   if (getMaxAmp() >= sidc.getSetup()->getSaturationAmplitude()) 
      gaussRegressionWasDone = doGaussRegression(sidc); 
 
   int time; 
   double thisamp, amp0, amp1, amp2, pulseAmp, pulsePos, pulseRms; 
   CCell *cell; 
       
   amp0 = 0;            // integral amplitude 
   amp1 = 0;            // first momentum 
   amp2 = 0;            // second momentum	  
 
   CBoolean subtractPedestals = sidc.getLookupItems().length()>0 ? True : False; 
    
   // 
   // Loop over timebins in pulse and sum up momenta. 
   // 
   for (time=getTStart(); time<=getTStop(); time++) {     
      cell = sidc.getCells()(getAnode(),time); 
      if (cell != 0) { 
	 if (subtractPedestals) 
	    thisamp = cell->getAmp() - sidc.getLookupItems()(getAnode())->getPedestal(); 
	 else 
	    thisamp = cell->getAmp(); 
	 amp0 += thisamp; 
	 amp1 += time * thisamp; 
	 amp2 += time * time * thisamp; 
      } 
   } 
 
   // 
   // Calculate pulse properties and apply gain correction. 
   // 
   pulseAmp = amp0 / sidc.getLookupItems()(getAnode())->getGain(); 
   pulsePos = amp1/amp0; 
   pulseRms = sqrt(amp2/amp0 - (amp1/amp0)*(amp1/amp0)); 
 
   setSumAmp(pulseAmp); 
    
   // 
   // Set pulse properties only if they make sence. 
   // 
   if (!gaussRegressionWasDone) { 
      if (pulseAmp>1. && pulseAmp<1.e5 && pulseRms>0.1) { 
	 setCenter(pulsePos); 
	 setRmsTime(pulseRms); 
	 setAmp(pulseAmp); 
	 setIsUsedForHit(False); 
      } 
      else { 
	 setCenter(-1); 
	 setRmsTime(-1); 
	 setAmp(-1); 
	 setSumAmp(-1); 
	 setIsUsedForHit(True); 
	 return False; 
      } 
   } 
    
   return True; 
} 
 
CBoolean CCoGPulse::doGaussRegression(CSidc& sidc) 
{ 
   double s0   = 0.0; double s1  = 0.0; double s2   = 0.0; double s3  = 0.0;   double s4  = 0.0; 
   double sy   = 0.0; double sxy = 0.0; double sx2y = 0.0; double sy2 = 0.0;  
   double w    = 0.0; double z   = 0.0; 
    
   double gaussSigma = -1.0; 
   double gaussTime0 = -1.0; 
   double gaussAmpl  = -1.0; 
 
   CBoolean subtractPedestals = sidc.getLookupItems().length()>0 ? True : False; 
       
   for(int i=getTStart(); i<=getTStop(); i++){ 
      if (sidc.getCells()(getAnode(), i) && 
	  sidc.getCells()(getAnode(), i)->getAmp() <= sidc.getSetup()->getSaturationAmplitude()) { 
	 w = sidc.getCells()(getAnode(), i)->getAmp(); 
	 if (subtractPedestals) w -= sidc.getLookupItems()(getAnode())->getPedestal(); 
	 w /= sidc.getLookupItems()(getAnode())->getGain(); 
	 if (w > 0.0) { 
	    z     = log(w); 
	    s1   += double(i)       * w; 
	    s2   += double(i*i)     * w; 
	    s3   += double(i*i*i)   * w; 
	    s4   += double(i*i*i*i) * w; 
	    sy   += z   * w; 
	    sy2  += z*z * w; 
	    sxy  += double(i)   * z * w; 
	    sx2y += double(i*i) * z * w; 
	    s0   += w; 
	 } 
      } 
   } 
    
   double s1kv = s1*s1; 
   double s12  = s1*s2; 
   double s13  = s1*s3; 
   double s14  = s1*s4; 
   double s2kv = s2*s2; 
   double s23  = s2*s3; 
   double s24  = s2*s4; 
   double s3kv = s3*s3; 
   double del  = s0*(s24-s3kv)+s1*(s23-s14)+s2*(s13-s2kv); 
   double aa, bb, cc; 
   if(del>0.0) { 
      cc=((s24-s3kv)*sy+(s23-s14)*sxy+(s13-s2kv)*sx2y)/del; 
      bb=((s23-s14)*sy+(s0*s4-s2kv)*sxy+(s12-s0*s3)*sx2y)/del; 
      aa=((s13-s2kv)*sy+(s12-s0*s3)*sxy+(s0*s2-s1kv)*sx2y)/del; 
      if (aa<0.0){ 
	 gaussSigma = sqrt(-0.5/aa); 
	 gaussTime0 = -bb/2.0/aa; 
	 double determ = (cc-pow(bb,2)/4./aa); 
	 if(abs(determ)>1.0e-10 && abs(determ)<50.0) gaussAmpl = exp(determ); 
      } 
   } 
   // 
   // Make sure results are sensible 
   // 
   if(gaussAmpl<10. || gaussAmpl>1000. || gaussTime0<1. || gaussTime0>1000. || 
      gaussSigma<.1 || gaussSigma>100.) { 
      return False; 
   } 
   else { 
      double gaussArea = sqrt(2*Pi) * gaussSigma * gaussAmpl; 
      setCenter(gaussTime0); 
      setAmp(gaussArea); 
      setMaxAmp(gaussAmpl); 
      setRmsTime(gaussSigma); 
      return True; 
   } 
} 
 
CCoGPulse* CCoGPulse::split(const CCyclicCollection<CCell>& cells) 
{ 
   // 
   // Scan this pulse for a local minimum. As soon as a minimum is found, the 
   // part before this is chopped off and returned as a new pulse. 
   // Note that the local minimum itself is not used for either pulse. 
   // If no minimum was found the zero pointer is returned. 
   // 
 
   double oldAmp = cells(getAnode(), getTStart())->getAmp();    
   CCoGPulse *newPulse = 0; 
   for (int time=getTStart()+1; time<=getTStop(); time++) { 
      // 
      // check if this cell is a local minimum in the pulse 
      // 
      if (cells(getAnode(),time-1)->getAmp() > cells(getAnode(),time)->getAmp() && 
	  cells(getAnode(),time+1)->getAmp() > cells(getAnode(),time)->getAmp()) {	  
	 newPulse = new CCoGPulse(getAnode(), getTStart()); 
	 newPulse->setTStop(time-1); 
	 setTStart(time+1); 
	 return newPulse; 
      }  
   } 
   return newPulse;    
} 

Back to index