//____________________________________________________________________
// 
//  TPC Helper class : Try to fit a TPC Cluster 
// 

//____________________________________________________________________
//
// $Id: BrTpcFitCluster.cxx,v 1.2 2001/06/22 17:45:24 cholm Exp $
// $Author: cholm $
// $Date: 2001/06/22 17:45:24 $
// $Copyright: (C) 2001 BRAHMS Collaboration <brahmlib@rhic.bnl.gov>
//
#ifndef BRAT_BrTpcFitCluster
#include "BrTpcFitCluster.h"
#endif
#ifndef ROOT_TMath
#include "TMath.h"
#endif
#ifndef BRAT_BrTpcSequence
#include "BrTpcSequence.h"
#endif
#ifndef WIN32
#include <iostream>
#else
#include <iostream.h>
#endif
#include <cassert>

//____________________________________________________________________
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, 
		Double_t *par, Int_t iflag)
{
  // static function so it can only operate on static objects, meaning
  // that all objects/variables must be known. Trick is to make a
  // static objectpointer
  // calculate chisquare
  
  f = BrTpcFitCluster::Instance()->ChiSquare(par);
}

// Set values of static members
//____________________________________________________________________
BrTpcFitCluster* BrTpcFitCluster::fgInstance = 0;

//____________________________________________________________________
Int_t BrTpcFitCluster::fgMaxFitPoints = 500;

//____________________________________________________________________

ClassImp(BrTpcFitCluster);

//____________________________________________________________________
 BrTpcFitCluster::BrTpcFitCluster()
{
  // Default constructor
  // Use this one 
  fNfitPoints = 0;
  fMinuit     = TVirtualFitter::Fitter(this, 5); // new TMinuit(5);  
  fMinuit->SetFCN(fcn);
  fZ          = new Double_t[fgMaxFitPoints];
  fEZ         = new Double_t[fgMaxFitPoints];
  fX          = new Double_t[fgMaxFitPoints];
  fY          = new Double_t[fgMaxFitPoints];
  
  for(Int_t i = 0; i < fgMaxFitPoints; i++) {
    fZ[i]  = 0;
    fEZ[i] = 0;
    fX[i]  = 0;
    fY[i]  = 0;
  }
  fNfitPoints = 0;

}

//____________________________________________________________________
 BrTpcFitCluster::~BrTpcFitCluster()
{
  // Dtor
  delete fMinuit;  
  
  delete [] fZ;
  delete [] fEZ;
  delete [] fX;
  delete [] fY;
}

//_____________________________________________________________________
void 
 BrTpcFitCluster::AddDataPoint(Double_t zVal, Double_t ezVal, 
			      Double_t xVal, Double_t yVal)
{
  if(fNfitPoints>=500) {
    cout << "The data array is full" << endl;
    return;
  }
  
  fZ[fNfitPoints]  = zVal;
  fEZ[fNfitPoints] = ezVal;
  fX[fNfitPoints]  = xVal;
  fY[fNfitPoints]  = yVal;

  fNfitPoints++;
}

//_____________________________________________________________________
 void BrTpcFitCluster::FitCluster(BrTpcCluster *cluster)
{
  // sets the data points in the sequences as data points in this
  // class and fits the cluster
  fCluster = cluster;
  
  Clear();

  // Set the datapoints by looping over cluster sequences
  const Int_t nSeqs = fCluster->GetNseqs();
  
  for(Int_t i = 0; i < nSeqs; i++) {
    
    BrTpcSequence *seq = fCluster->GetSeq(i);

    const Int_t nBins = seq->GetNseq();
    
    for(Int_t j = 0; j < nBins; j++) {
      
      AddDataPoint(seq->GetAdc(j), 1, seq->GetPad(), seq->GetTime());
    }
  }
  
  cout << "Number of fit points : " << fNfitPoints << endl;

  Fit();
}
  
//_____________________________________________________________________
 void BrTpcFitCluster::Clear()
{
  // Set the used array memebres to 0 and number of fit points to 0
  for(Int_t i = 0; i < fNfitPoints; i++) {
    fZ[i]  = 0;
    fEZ[i] = 0;
    fX[i]  = 0;
    fY[i]  = 0;
  }
  
  fNfitPoints = 0;
  
}

//_____________________________________________________________________
 void BrTpcFitCluster::Fit()
{
  fgInstance = this;
  
  Double_t arglist[10];
  Int_t    ierflg = 0;

  arglist[0] = 1;

  ierflg = fMinuit->ExecuteCommand("SET ERR", arglist ,1);
  
  // Set starting values and step sizes for parameters
  
  Double_t vstart[5] = { 
    fCluster->GetMaxADC(), fCluster->GetMeanPad(),
    fCluster->GetVarPad(), fCluster->GetMeanTime(),
    fCluster->GetVarTime()};

  Double_t step[5] = { 0.1 , 0.1, 0.01, 0.1, 0.01 };
  ierflg = fMinuit->SetParameter(0, "max adc" , vstart[0], step[0], 0, 0);
  ierflg = fMinuit->SetParameter(1, "pad pos" , vstart[1], step[1], 0, 0);
  ierflg = fMinuit->SetParameter(2, "pad var" , vstart[2], step[2], 0, 0);
  ierflg = fMinuit->SetParameter(3, "time pos", vstart[3], step[3], 0, 0);
  ierflg = fMinuit->SetParameter(4, "time var", vstart[4], step[4], 0, 0);
  
  // Now ready for minimization step
  arglist[0] = 500;
  arglist[1] = 1.;
  //  fMinuit->mnexcm("SIMPLEX", arglist ,2,ierflg);
  ierflg = fMinuit->ExecuteCommand("MIGRAD", arglist ,2);
  
  // Print results
  Double_t amin,  edm,   errdef; 
  Int_t    nvpar, nparx, icstat;
  icstat = fMinuit->GetStats(amin,edm,errdef,nvpar,nparx);
  fMinuit->PrintResults(3,amin);
  
}

//_________________________________________________________________________
 Double_t BrTpcFitCluster::ChiSquare(Double_t *par)
{
  // This function calculates the chisquare for all fit points
  Double_t chisq = 0;
   
  for (Int_t i = 0; i < fNfitPoints; i++) {

    Double_t delta  = (Func(fX[i], fY[i], par) - fZ[i])/fEZ[i];
    chisq += delta*delta;
  }
  
  return chisq;
}

//_________________________________________________________________________
 Double_t BrTpcFitCluster::Func(Double_t x, Double_t y, Double_t *par)
{
  // This method calculates the function value in each points
  // Function is assumed to be a doubble Gaus
  Double_t value= par[0]*TMath::Gaus(x, par[1], par[2])
    *TMath::Gaus(y, par[3], par[4]);
  
  return value;
}

Copyright ; 2002 BRAHMS Collaboration <brahmlib@rcf.rhic.bnl.gov>
Last Update on by

Validate HTML
Validate CSS