merge.cpp

Go to the documentation of this file.

#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
#include <vector>
 
#include "TCanvas.h"
#include "TFile.h"
#include "TFolder.h"
#include "TLatex.h"
#include "TList.h"
#include "TObjArray.h"
#include "TPostScript.h"
#include "TProfile.h"
#include "TROOT.h"
#include "TSpline.h"
#include "TStyle.h"
#include "TTree.h"
 
#include "include/GrXtCalib.h"
#include "include/IniCalib.h"
#include "include/MdcCalibConst.h"
#include "include/MdcCosGeom.h"
#include "include/PreT0Calib.h"
#include "include/PreXtCalib.h"
#include "include/QtCalib.h"
#include "include/T0Calib.h"
#include "include/XtCalib.h"
#include "include/XtInteCalib.h"
#include "include/fun.h"
 
using namespace std;
 
int main( int argc, char* argv[] ) {
  char* jobname;
  char* strcal;
  int   fgCal = 1;
  if ( argc > 2 )
  {
    jobname = argv[1];
    strcal  = argv[2];
    sscanf( strcal, "%d", &fgCal );
  }
  else if ( argc > 1 ) { jobname = argv[1]; }
  else
  {
    cout << "bad argument" << endl;
    return -1;
  }
  if ( fgCal <= 0 ) cout << "do not calibrate " << endl;
 
  string path   = "";
  string strJob = jobname;
  cout << "strJob: " << strJob << endl;
  string::size_type ilast = strJob.find_last_of( "/" );
  if ( string::npos != ilast ) { path = strJob.substr( 0, ilast ); }
 
  int      calType;
  string   constname;
  string   confname;
  string   str;
  string   strtmp;
  ifstream fjob( jobname );
  if ( !fjob.is_open() )
  {
    cout << "ERROR: can not read jobOption: " << jobname << endl;
    return 0;
  }
  else
  {
    cout << "Open jobOption: " << jobname << endl;
    while ( getline( fjob, str ) )
    {
      if ( str.find( "//", 0 ) != string::npos ) { continue; }
      else if ( str.find( "CalibRootCnvSvc.Mdcrootfile", 0 ) != string::npos )
      {
        string::size_type i1 = str.find_first_of( "\"" );
        string::size_type i2 = str.find_last_of( "\"" );
        constname            = str.substr( i1 + 1, i2 - i1 - 1 );
      }
      else if ( str.find( "MdcCalibAlg.ConfigFile", 0 ) != string::npos )
      {
        string::size_type i1 = str.find_first_of( "\"" );
        string::size_type i2 = str.find_last_of( "\"" );
        confname             = str.substr( i1 + 1, i2 - i1 - 1 );
      }
      else if ( str.find( "MdcCalibAlg.MdcCalFlg", 0 ) != string::npos )
      {
        string::size_type i1 = str.find_first_of( "=" );
        string::size_type i2 = str.find_last_of( ";" );
        strtmp               = str.substr( i1 + 1, i2 - i1 - 1 );
        sscanf( strtmp.c_str(), "%d", &calType );
      }
    }
  }
 
  MdcCosGeom* pGeom = 0;
  pGeom             = new MdcCosGeom( "/home/bes/wulh/document/wireconf.txt",
                                      "/home/bes/wulh/calibConst/MdcAlignPar_ini.txt" );
  pGeom->initialize( 0.0 );
 
  int lay;
  for ( lay = 0; lay < 15; lay++ ) gNEntr[lay] = 1;
  for ( lay = 15; lay < 43; lay++ ) gNEntr[lay] = 2;
 
  // read config file
  string   strconfig;
  string   strcomment;
  ifstream fconfig( confname.c_str() );
  if ( !fconfig.is_open() )
  {
    cout << "ERROR: can not read config file" << endl;
    return 0;
  }
  else
  {
    cout << "Open config file: " << confname << endl;
    while ( fconfig >> strconfig )
    {
      if ( '#' == strconfig[0] ) { getline( fconfig, strcomment ); }
      else if ( "TimeShift" == strconfig ) { fconfig >> gTimeShift; }
      else if ( "TesMin" == strconfig ) { fconfig >> gTesMin; }
      else if ( "TesMax" == strconfig ) { fconfig >> gTesMax; }
      else if ( "FlagIniCalConst" == strconfig ) { fconfig >> gFgIniCalConst; }
      else if ( "FlagUpdateTmInPreT0" == strconfig ) { fconfig >> gPreT0SetTm; }
      else if ( "InitT0" == strconfig ) { fconfig >> gInitT0; }
      else if ( "T0Shift" == strconfig ) { fconfig >> gT0Shift; }
      else if ( "TminFitChindf" == strconfig ) { fconfig >> gTminFitChindf; }
      else if ( "TmaxFitChindf" == strconfig ) { fconfig >> gTmaxFitChindf; }
      else if ( "ResidualType" == strconfig ) { fconfig >> gResiType; }
      else if ( "UpdateSigma" == strconfig ) { fconfig >> gCalSigma; }
      else if ( "FixXtC0" == strconfig ) { fconfig >> gFixXtC0; }
      else if ( "RawTimeFitRange" == strconfig )
      {
        for ( lay = 0; lay < NLAYER; lay++ )
        {
          fconfig >> strtmp >> gFgCalib[lay] >> gTminFitRange[lay][0] >>
              gTminFitRange[lay][1] >> gTmaxFitRange[lay][0] >> gTmaxFitRange[lay][1] >>
              gInitTm[lay] >> strtmp >> gQmin[lay] >> gQmax[lay];
        }
      }
    }
  }
  fconfig.close();
 
  TObjArray* hlist = new TObjArray( 0 );
  CalibBase* pcal;
  if ( 0 == calType ) pcal = new IniCalib();
  else if ( 1 == calType ) pcal = new PreXtCalib();
  else if ( 2 == calType ) pcal = new PreT0Calib();
  else if ( 3 == calType ) pcal = new XtCalib();
  else if ( 4 == calType ) pcal = new GrXtCalib();
  else if ( 9 == calType ) pcal = new XtInteCalib();
  else if ( 5 == calType ) pcal = new T0Calib();
  else if ( 8 == calType ) pcal = new QtCalib();
  else
  {
    cout << "Error CalibType" << endl;
    return 0;
  }
  pcal->init( hlist, pGeom );
 
  vector<string> fhistname = getHistList();
  if ( 0 == fhistname.size() )
  {
    cout << "hist file path: " << path << endl;
    fhistname = getHistList( path );
  }
  for ( unsigned nf = 0; nf < fhistname.size(); nf++ )
  {
    TFile* fin = new TFile( fhistname[nf].c_str() );
    if ( !fin->IsOpen() ) { continue; }
    else
    {
      cout << "merge hist file " << nf << ": " << fhistname[nf] << endl;
      pcal->mergeHist( fin );
      fin->Close();
    }
  }
 
  // read calib const.
  MdcCalibConst* calconst = new MdcCalibConst();
  //      if(fgCal <= 0){constname =
  //      "/home/bes/wulh/calibConst/MdcCalibConst_11397_psi3770_652b_v1.root";}
  if ( constname == "" )
  { constname = "/home/bes/wulh/calibConst/MdcCalibConst_11397_psi3770_652b_v1.root"; }
  TFile fconst( constname.c_str() );
  if ( !fconst.IsOpen() )
  {
    cout << "ERROR: " << constname << " does not exist" << endl;
    return 0;
  }
  cout << "Open calib const file " << constname << endl;
  int    key;
  double xtpar;
  int    entry;
  TTree* xttree = (TTree*)fconst.Get( "XtTree" );
  xttree->SetBranchAddress( "xtkey", &key );
  xttree->SetBranchAddress( "xtpar", &xtpar );
  entry = (int)xttree->GetEntries();
  for ( int i = 0; i < entry; i++ )
  {
    xttree->GetEntry( i );
    calconst->fillXtpar( key, xtpar );
  }
 
  double t0;
  double delt0;
  TTree* t0tree = (TTree*)fconst.Get( "T0Tree" );
  t0tree->SetBranchAddress( "t0", &t0 );
  t0tree->SetBranchAddress( "delt0", &delt0 );
  entry = (int)t0tree->GetEntries();
  for ( int i = 0; i < entry; i++ )
  {
    t0tree->GetEntry( i );
    calconst->fillT0( t0 );
    calconst->fillDelT0( delt0 );
  }
 
  double qtpar0;
  double qtpar1;
  TTree* qttree = (TTree*)fconst.Get( "QtTree" );
  qttree->SetBranchAddress( "qtpar0", &qtpar0 );
  qttree->SetBranchAddress( "qtpar1", &qtpar1 );
  entry = (int)qttree->GetEntries();
  for ( int i = 0; i < entry; i++ )
  {
    qttree->GetEntry( i );
    calconst->fillQtpar0( qtpar0 );
    calconst->fillQtpar1( qtpar1 );
  }
 
  double sdpar;
  TTree* sdtree = (TTree*)fconst.Get( "SdTree" );
  sdtree->SetBranchAddress( "sdkey", &key );
  sdtree->SetBranchAddress( "sdpar", &sdpar );
  entry = sdtree->GetEntries();
  for ( int i = 0; i < entry; i++ )
  {
    sdtree->GetEntry( i );
    calconst->fillSdpar( key, sdpar );
  }
 
  TObjArray* newXtList = new TObjArray( 0 );
  TList*     list      = fconst.GetListOfKeys();
  int        listSize  = (int)( list->GetSize() );
  cout << "Number of trees in old calib const file: " << listSize << endl;
  if ( listSize > 4 )
  {
    gROOT->cd();
    for ( int i = 0; i < listSize; i++ )
    {
      TTree* tree = (TTree*)fconst.Get( list->At( i )->GetName() );
      string strName( tree->GetName() );
      if ( string::npos != strName.find( "trNewXt" ) )
      {
        TTree* t2 = tree->CopyTree( "" );
        newXtList->Add( t2 );
      }
    }
  }
  fconst.Close();
 
  gStyle->SetCanvasBorderMode( 0 );
  gStyle->SetCanvasColor( 10 );
  gStyle->SetOptFit( 0011 );
  gStyle->SetStatColor( 10 );
  gStyle->SetStatBorderSize( 1 );
  gStyle->SetStatFont( 42 );
  gStyle->SetStatFontSize( 0.04 );
  gStyle->SetStatX( 0.9 );
  gStyle->SetStatY( 0.9 );
  gStyle->SetPadColor( 10 );
  gStyle->SetFuncColor( 2 );
 
  TObjArray* r2tList = new TObjArray( 0 );
 
  // x-t calibration
  pcal->calib( calconst, newXtList, r2tList );
 
  TFile fhist( "histall.root", "recreate" );
  fhist.cd();
  hlist->Write();
  fhist.Close();
 
  // output new calib const file
  if ( fgCal > 0 ) writeConst( calconst, newXtList, r2tList );
 
  return 0;
}