MdcCalibAlg/share/distcalib/src/fun.cpp

Go to the documentation of this file.

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <iostream>
#include <sstream>
#include <vector>
 
#include "TFile.h"
#include "TTree.h"
 
#include "include/MdcCalibConst.h"
#include "include/fun.h"
 
using namespace std;
 
vector<double> XMEAS;
vector<double> TBINCEN;
vector<double> ERR;
double         Tmax;
double         Dmax;
vector<double> XMEASED;
vector<double> TBINCENED;
vector<double> ERRED;
 
int gNEntr[43];
 
/* from calib config file */
double gTimeShift     = 0.0;    /* if T<0 after subtracting Tes, use this */
double gTesMin        = 0.0;    /* minimun Tes for calibration */
double gTesMax        = 9999.0; /* maximun Tes for calibration */
int    gFgIniCalConst = 2;      /* effective for IniMdcCalib */
bool   gPreT0SetTm    = true;   /* flag for updating Tm in PreT0Calib */
double gInitT0        = 50.0;   /* initial value of T0 fit */
double gT0Shift       = 0.0;    /* t0 shift based on leading edge fitting */
double gTminFitChindf = 20.0;   /* chisquare cut for Tmin fit */
double gTmaxFitChindf = 20.0;   /* chisquare cut for Tmax fit */
int    gResiType      = 0;      /* 0: including measurement point; 1: excluding */
int    gCalSigma      = 1;      /* flag for update sigma */
int    gFixXtC0       = 0;      /* 1: fix c0 at 0 */
int    gFgCalib[NLAYER];
double gTminFitRange[NLAYER][2];
double gTmaxFitRange[NLAYER][2];
double gInitTm[NLAYER]; /* initial value of Tm fit */
double gQmin[NLAYER];   /* QT fit range */
double gQmax[NLAYER];   /* QT fit range */
 
double xtFun( double t, double xtpar[] ) {
  int    ord;
  double dist = 0.0;
  double tm   = xtpar[6];
 
  if ( t < tm )
  {
    for ( ord = 0; ord <= 5; ord++ ) { dist += xtpar[ord] * pow( t, ord ); }
  }
  else
  {
    for ( ord = 0; ord <= 5; ord++ ) { dist += xtpar[ord] * pow( tm, ord ); }
    dist += xtpar[7] * ( t - tm );
  }
 
  return dist;
}
 
void fcnXT( Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag ) {
  unsigned int i;
  int          ord;
  Double_t     deta;
  Double_t     chisq = 0.0;
  Double_t     dfit;
 
  for ( i = 0; i < TBINCEN.size(); i++ )
  {
    dfit = 0;
    for ( ord = 0; ord <= 5; ord++ ) { dfit += par[ord] * pow( TBINCEN[i], ord ); }
    deta = ( dfit - XMEAS[i] ) / ERR[i];
    chisq += deta * deta;
  }
 
  f = chisq;
}
 
void fcnXtEdge( Int_t& npar, Double_t* gin, Double_t& f, Double_t* par, Int_t iflag ) {
  unsigned int i;
  Double_t     deta;
  Double_t     chisq = 0.0;
  Double_t     dfit;
 
  for ( i = 0; i < TBINCENED.size(); i++ )
  {
    dfit = par[0] * ( TBINCENED[i] - Tmax ) + Dmax;
    deta = ( dfit - XMEASED[i] ) / ERRED[i];
    chisq += deta * deta;
  }
 
  f = chisq;
}
 
Double_t xtFitFun( Double_t* x, Double_t par[] ) {
  Double_t val = 0.0;
  for ( Int_t ord = 0; ord < 6; ord++ ) { val += par[ord] * pow( x[0], ord ); }
  return val;
}
 
Double_t xtFitEdge( Double_t* x, Double_t par[] ) {
  double val = Dmax + ( x[0] - Tmax ) * par[0];
  return val;
}
 
void writeConst( MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList ) {
  TFile fout( "MdcCalibConst_new.root", "recreate" );
 
  int    key;
  double xtpar;
  TTree* xttree = new TTree( "XtTree", "XtTree" );
  xttree->Branch( "xtkey", &key, "key/I" );
  xttree->Branch( "xtpar", &xtpar, "xtpar/D" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    for ( int entr = 0; entr < 18; entr++ )
    {
      for ( int lr = 0; lr < 3; lr++ )
      {
        for ( int ord = 0; ord < 8; ord++ )
        {
          key   = calconst->getXtKey( lay, entr, lr, ord );
          xtpar = calconst->getXtpar( lay, entr, lr, ord );
          xttree->Fill();
        }
      }
    }
  }
 
  double t0;
  double delt0;
  TTree* t0tree = new TTree( "T0Tree", "T0Tree" );
  t0tree->Branch( "t0", &t0, "t0/D" );
  t0tree->Branch( "delt0", &delt0, "delt0/D" );
  for ( int wid = 0; wid < 6796; wid++ )
  {
    t0    = calconst->getT0( wid );
    delt0 = calconst->getDelT0( wid );
    t0tree->Fill();
  }
 
  double qtval[2];
  TTree* qttree = new TTree( "QtTree", "QtTree" );
  qttree->Branch( "qtpar0", &( qtval[0] ), "qtpar0/D" );
  qttree->Branch( "qtpar1", &( qtval[1] ), "qtpar1/D" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    qtval[0] = calconst->getQtpar0( lay );
    qtval[1] = calconst->getQtpar1( lay );
    qttree->Fill();
  }
 
  double sdpar;
  TTree* sdtree = new TTree( "SdTree", "SdTree" );
  sdtree->Branch( "sdkey", &key, "key/I" );
  sdtree->Branch( "sdpar", &sdpar, "sdpar/D" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    for ( int entr = 0; entr < 6; entr++ )
    {
      for ( int lr = 0; lr < 2; lr++ )
      {
        for ( int bin = 0; bin < 24; bin++ )
        {
          key   = calconst->getSdKey( lay, entr, lr, bin );
          sdpar = calconst->getSdpar( lay, entr, lr, bin );
          sdtree->Fill();
        }
      }
    }
  }
 
  fout.cd();
  xttree->Write();
  t0tree->Write();
  qttree->Write();
  sdtree->Write();
  if ( ( newXtList->GetEntries() ) > 0 ) newXtList->Write();
  if ( ( r2tList->GetEntries() ) > 0 ) r2tList->Write();
  fout.Close();
}
 
vector<string> getHistList() {
  vector<string> fnames;
 
  string command( "JobOutputDir=`/bin/ls -dt1 joboutput-* 2>/dev/null | head -1`\n"
                  "if [ -d \"${JobOutputDir}\" ]; then\n"
                  "  find ${JobOutputDir} -name hist.root\n"
                  "fi\n" );
 
  stringstream fnstream;
 
  char* fnbuf   = new char[1024];
  FILE* fstream = popen( command.c_str(), "r" );
 
  while ( fgets( fnbuf, 1024, fstream ) != NULL ) { fnstream << fnbuf; }
 
  string fname;
  while ( !( fnstream >> fname ).eof() ) { fnames.push_back( fname ); }
 
  pclose( fstream );
  delete[] fnbuf;
 
  if ( fnames.empty() )
  {
    cout << "WARNING: Failed to retrieve hist files in the current directory!" << endl;
    //       exit(1);
  }
  return fnames;
}
 
vector<string> getHistList( string path ) {
  vector<string>    fnames;
  string            newpath = path;
  string::size_type strl    = newpath.length();
  if ( ( strl > 1 ) && ( '/' == newpath[strl - 1] ) ) newpath.erase( strl - 1 );
 
  char com1[500];
  sprintf( com1, "JobOutputDir=`/bin/ls -dt1 %s/joboutput-* 2>/dev/null | head -1`\n",
           newpath.c_str() );
  string       command1( com1 );
  string       command2( "if [ -d \"${JobOutputDir}\" ]; then\n"
                               "  find ${JobOutputDir} -name hist.root\n"
                               "fi\n" );
  string       command = command1 + command2;
  stringstream fnstream;
 
  char* fnbuf   = new char[1024];
  FILE* fstream = popen( command.c_str(), "r" );
 
  while ( fgets( fnbuf, 1024, fstream ) != NULL ) { fnstream << fnbuf; }
 
  string fname;
  while ( !( fnstream >> fname ).eof() ) { fnames.push_back( fname ); }
 
  pclose( fstream );
  delete[] fnbuf;
 
  if ( fnames.empty() )
  {
    cout << "ERROR: Failed to retrieve hist files!" << endl;
    exit( 1 );
  }
  return fnames;
}