XtCalib.cpp

Go to the documentation of this file.

#include "include/XtCalib.h"
#include "include/fun.h"
 
#include "TMinuit.h"
 
XtCalib::XtCalib() { cout << "Calibration type: XtCalib" << endl; }
 
XtCalib::~XtCalib() {}
 
void XtCalib::init( TObjArray* hlist, MdcCosGeom* pGeom ) {
  CalibBase::init( hlist, pGeom );
 
  m_fdXt = new TFolder( "mfdxt", "fdxt" );
  hlist->Add( m_fdXt );
 
  char hname[200];
  for ( int lay = 0; lay < 43; lay++ )
  {
    for ( int iEntr = 0; iEntr < gNEntr[lay]; iEntr++ )
    {
      for ( int lr = 0; lr < NLR; lr++ )
      {
        for ( int bin = 0; bin <= NXTBIN; bin++ )
        {
          sprintf( hname, "mHxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin );
          m_hxt[lay][iEntr][lr][bin] = new TH1D( hname, "", 300, -1.5, 1.5 );
          m_fdXt->Add( m_hxt[lay][iEntr][lr][bin] );
        }
      }
    }
  }
}
 
void XtCalib::mergeHist( TFile* fhist ) {
  CalibBase::mergeHist( fhist );
 
  char     hname[200];
  TFolder* fd = (TFolder*)fhist->Get( "fdXt" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    for ( int iEntr = 0; iEntr < gNEntr[lay]; iEntr++ )
    {
      for ( int lr = 0; lr < NLR; lr++ )
      {
        for ( int bin = 0; bin <= NXTBIN; bin++ )
        {
          sprintf( hname, "Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin );
          TH1F* hist = (TH1F*)fd->FindObjectAny( hname );
          m_hxt[lay][iEntr][lr][bin]->Add( hist );
          //             if((0==lay)&&(0==iEntr)&&(0==lr)&&(50==bin)){
          //                  cout << setw(15) << hist->GetEntries()
          //                   << setw(15) <<
          // m_hxt[lay][iEntr][lr][bin]->GetEntries() << endl;
          //             }
        }
      }
    }
  }
}
 
void XtCalib::calib( MdcCalibConst* calconst, TObjArray* newXtList, TObjArray* r2tList ) {
  CalibBase::calib( calconst, newXtList, r2tList );
 
  Int_t    ierflg;
  Int_t    istat;
  Int_t    nvpar;
  Int_t    nparx;
  Double_t fmin;
  Double_t edm;
  Double_t errdef;
  Double_t arglist[10];
 
  TMinuit* gmxt = new TMinuit( 6 );
  gmxt->SetPrintLevel( -1 );
  gmxt->SetFCN( fcnXT );
  gmxt->SetErrorDef( 1.0 );
  gmxt->mnparm( 0, "xtpar0", 0, 0.1, 0, 0, ierflg );
  gmxt->mnparm( 1, "xtpar1", 0, 0.1, 0, 0, ierflg );
  gmxt->mnparm( 2, "xtpar2", 0, 0.1, 0, 0, ierflg );
  gmxt->mnparm( 3, "xtpar3", 0, 0.1, 0, 0, ierflg );
  gmxt->mnparm( 4, "xtpar4", 0, 0.1, 0, 0, ierflg );
  gmxt->mnparm( 5, "xtpar5", 0, 0.1, 0, 0, ierflg );
  arglist[0] = 0;
  gmxt->mnexcm( "SET NOW", arglist, 0, ierflg );
 
  TMinuit* gmxtEd = new TMinuit( 1 );
  gmxtEd->SetPrintLevel( -1 );
  gmxtEd->SetFCN( fcnXtEdge );
  gmxtEd->SetErrorDef( 1.0 );
  gmxtEd->mnparm( 0, "xtpar0", 0, 0.1, 0, 0, ierflg );
  arglist[0] = 0;
  gmxtEd->mnexcm( "SET NOW", arglist, 0, ierflg );
 
  //      double xtpar;
  int      i;
  Stat_t   histEntry;
  double   xtpar;
  double   xterr;
  double   tbcen;
  double   deltx;
  double   xcor;
  double   xerr;
  double   xtini[8];
  double   xtfit[8];
  ofstream fxtlog( "xtlog" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    if ( 0 == gFgCalib[lay] ) continue;
    for ( int iEntr = 0; iEntr < gNEntr[lay]; iEntr++ )
    {
      for ( int iLR = 0; iLR < NLR; iLR++ )
      {
        fxtlog << "Layer " << setw( 3 ) << lay << setw( 3 ) << iEntr << setw( 3 ) << iLR
               << endl;
        for ( int ord = 0; ord < NXTPAR; ord++ )
        {
          if ( 0 == iEntr ) xtpar = calconst->getXtpar( lay, 8, iLR, ord );
          else if ( 1 == iEntr ) xtpar = calconst->getXtpar( lay, 9, iLR, ord );
          xtini[ord] = xtpar;
          xtfit[ord] = xtpar;
        }
        Tmax = xtini[6];
 
        for ( int bin = 0; bin <= NXTBIN; bin++ )
        {
          histEntry = (int)( m_hxt[lay][iEntr][iLR][bin]->GetEntries() );
          if ( histEntry > 100 )
          {
            deltx = m_hxt[lay][iEntr][iLR][bin]->GetMean();
            xerr  = m_hxt[lay][iEntr][iLR][bin]->GetRMS();
          }
          else { continue; }
 
          if ( bin < NXTBIN ) tbcen = ( (double)bin + 0.5 ) * gTbinw;
          else tbcen = xtini[6]; // m_tm[lay][iEntr][iLR];
          xcor = xtFun( tbcen, xtini ) - deltx;
 
          if ( ( tbcen <= Tmax ) || ( bin == NXTBIN ) )
          {
            TBINCEN.push_back( tbcen );
            XMEAS.push_back( xcor );
            ERR.push_back( xerr );
          }
          else
          {
            TBINCENED.push_back( tbcen );
            XMEASED.push_back( xcor );
            ERRED.push_back( xerr );
          }
          fxtlog << setw( 3 ) << bin << setw( 15 ) << deltx << setw( 15 ) << xcor << setw( 15 )
                 << tbcen << setw( 15 ) << xerr << endl;
        } // end of bin loop
 
        if ( XMEAS.size() < 12 )
        {
          TBINCEN.clear();
          XMEAS.clear();
          ERR.clear();
 
          TBINCENED.clear();
          XMEASED.clear();
          ERRED.clear();
 
          continue;
        }
 
        for ( int ord = 0; ord <= 5; ord++ )
        {
          arglist[0] = ord + 1;
          arglist[1] = xtini[ord];
          gmxt->mnexcm( "SET PARameter", arglist, 2, ierflg );
        }
 
        // fix the xtpar[0] at 0
        if ( 1 == gfixXtC0 )
        {
          arglist[0] = 1;
          arglist[1] = 0.0;
          gmxt->mnexcm( "SET PARameter", arglist, 2, ierflg );
          gmxt->mnexcm( "FIX", arglist, 1, ierflg );
        }
 
        arglist[0] = 1000;
        arglist[1] = 0.1;
        gmxt->mnexcm( "MIGRAD", arglist, 2, ierflg );
        gmxt->mnstat( fmin, edm, errdef, nvpar, nparx, istat );
 
        fxtlog << "Xtpar: " << endl;
        if ( ( 0 == ierflg ) && ( istat >= 2 ) )
        {
          for ( int ord = 0; ord <= 5; ord++ )
          {
            gmxt->GetParameter( ord, xtpar, xterr );
            //                calconst -> resetXtpar(lay, iEntr, iLR, ord, xtpar);
            xtfit[ord] = xtpar;
 
            if ( 1 == gNEntr[lay] )
            {
              for ( i = 0; i < 18; i++ ) calconst->resetXtpar( lay, i, iLR, ord, xtpar );
            }
            else if ( 2 == gNEntr[lay] )
            {
              if ( 0 == iEntr )
              {
                for ( i = 0; i < 9; i++ ) // entr<0
                  calconst->resetXtpar( lay, i, iLR, ord, xtpar );
              }
              else
              {
                for ( i = 9; i < 18; i++ ) // entr>0
                  calconst->resetXtpar( lay, i, iLR, ord, xtpar );
              }
            }
            fxtlog << setw( 15 ) << xtpar << setw( 15 ) << xterr << endl;
          }
        }
        else
        {
          for ( int ord = 0; ord <= 5; ord++ )
          { fxtlog << setw( 15 ) << xtini[ord] << setw( 15 ) << "0" << endl; }
        }
        fxtlog << setw( 15 ) << Tmax << setw( 15 ) << "0" << endl;
 
        //   release the first parameter
        if ( 1 == gfixXtC0 )
        {
          arglist[0] = 1;
          gmxt->mnexcm( "REL", arglist, 1, ierflg );
        }
 
        Dmax = xtFun( Tmax, xtfit );
 
        if ( XMEASED.size() >= 3 )
        {
          // fit xt in the edge area
          arglist[0] = 1;
          arglist[1] = xtini[7];
          gmxtEd->mnexcm( "SET PARameter", arglist, 2, ierflg );
 
          arglist[0] = 1000;
          arglist[1] = 0.1;
          gmxtEd->mnexcm( "MIGRAD", arglist, 2, ierflg );
          gmxtEd->mnstat( fmin, edm, errdef, nvpar, nparx, istat );
 
          if ( ( 0 == ierflg ) && ( istat >= 2 ) )
          {
            gmxtEd->GetParameter( 0, xtpar, xterr );
            if ( xtpar < 0.0 ) xtpar = 0.0;
            //                calconst -> resetXtpar(lay, iEntr, iLR, 7, xtpar);
 
            if ( 1 == gNEntr[lay] )
            {
              for ( i = 0; i < 18; i++ ) calconst->resetXtpar( lay, i, iLR, 7, xtpar );
            }
            else if ( 2 == gNEntr[lay] )
            {
              if ( 0 == iEntr )
              {
                for ( i = 0; i < 9; i++ ) calconst->resetXtpar( lay, i, iLR, 7, xtpar );
              }
              else
              {
                for ( i = 9; i < 18; i++ ) calconst->resetXtpar( lay, i, iLR, 7, xtpar );
              }
            }
            fxtlog << setw( 15 ) << xtpar << setw( 15 ) << xterr << endl;
          }
          else { fxtlog << setw( 15 ) << xtini[7] << setw( 15 ) << "0" << endl; }
        }
        else { fxtlog << setw( 15 ) << xtini[7] << setw( 15 ) << "0" << endl; }
        fxtlog << "Tm " << setw( 15 ) << Tmax << "  Dmax " << setw( 15 ) << Dmax << endl;
 
        TBINCEN.clear();
        XMEAS.clear();
        ERR.clear();
        TBINCENED.clear();
        XMEASED.clear();
        ERRED.clear();
      } // lr loop
    }   // entrance loop
  }     // layer loop
  fxtlog.close();
 
  renameHist();
  delete gmxt;
  delete gmxtEd;
}
 
void XtCalib::renameHist() {
  char hname[200];
  m_fdXt->SetName( "fdXt" );
  for ( int lay = 0; lay < 43; lay++ )
  {
    for ( int iEntr = 0; iEntr < gNEntr[lay]; iEntr++ )
    {
      for ( int lr = 0; lr < NLR; lr++ )
      {
        for ( int bin = 0; bin <= NXTBIN; bin++ )
        {
          sprintf( hname, "Hxt%02d_E%02d_LR%01d_B%02d", lay, iEntr, lr, bin );
          m_hxt[lay][iEntr][lr][bin]->SetName( hname );
        }
      }
    }
  }
}