/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Emc/EmcCalib/EmcBhaCalib/EmcBhaCalib-00-00-34/src/EmcBhaCalib.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//
// Description:
//      Class EmcBhaCalib - performs calibration of EMC Xtals with Bhabha
//      events and a Chi^2 fit, the resulting system of linear equations 
//      of the fit is solved with the SLAP Algebra package
//
// Environment:
//      Software developed for the BESIII Detector at the BEPCII.
//
// Author List:
//      Chunxiu Liu
//
// Copyright Information:
//      Copyright (C) 2005               IHEP
//
//------------------------------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "EmcBhaCalib/EmcBhaCalib.h"

//-------------
// C Headers --
//-------------
extern "C" {
#include "slap/dlap.h"
}
#include <iostream>
#include <fstream>
#include <cmath>
#include <cassert>
#include <cstdlib>
//-------------------------------
// Collaborating Class Headers --
//-------------------------------

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/DataObject.h"

#include "CLHEP/Vector/ThreeVector.h"
using namespace std;

using CLHEP::Hep3Vector;

//--------------------
#include "GaudiKernel/MsgStream.h"


#include "TROOT.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1F.h"
#include "TObjArray.h"



//----------------
// Constructors --
//----------------
EmcBhaCalib::EmcBhaCalib(const std::string& name, ISvcLocator* pSvcLocator)
  :Algorithm(name, pSvcLocator),  
   m_dirHitsCut(200),
   m_convCrit(0.000001),
   m_askForMatrixInversion(true),
   m_fitPolynom(true),  //no used now
   m_writeToFile(false),
   m_readDataFromDB(false),
   m_equationSolver("GMR"),
   m_fileExt(""),
   m_fileDir("/ihepbatch/besdata/public/liucx/Calib/"),
   m_nrNonZeros(0), 
   m_nrXtalsEnoughHits(0),
   m_runNumberFile("runnumbers.dat"),
   m_MsgFlag(0)

{

  // Declare the properties  
  declareProperty("equationSolver", m_equationSolver);
  declareProperty("dirHitsCut", m_dirHitsCut);
  declareProperty("writeToFile", m_writeToFile);
  declareProperty("fileExt", m_fileExt);
  declareProperty("fileDir", m_fileDir);
  declareProperty("readDataFromDB", m_readDataFromDB);
  declareProperty("askForMatrixInversion", m_askForMatrixInversion);
  declareProperty("fitPolynom", m_fitPolynom);
  declareProperty("convCrit",  m_convCrit);
  declareProperty("runNumberFile", m_runNumberFile);
  declareProperty("MsgFlag", m_MsgFlag);

  m_calibConst = new double[6240];
  m_calibConstUnred = new double[6240];
  m_absoluteConstants = new double[6240];
  m_oldConstants = new double[6240];


  //ntuple
  m_tuple1=0;
 


}

//--------------
// Destructor --
//--------------
EmcBhaCalib::~EmcBhaCalib() {
  if ( 0 != m_calibConst) {
    delete [] m_calibConst;
    m_calibConst = 0;
  }
  if ( 0 != m_calibConstUnred) {
    delete [] m_calibConstUnred;
    m_calibConstUnred = 0;
  }
  if ( 0 != m_absoluteConstants) {
    delete [] m_absoluteConstants;
    m_absoluteConstants = 0;
  }
  if ( 0 != m_oldConstants) {
    delete [] m_oldConstants;
    m_oldConstants = 0;
  }
  if ( 0 != m_myCalibData) {
    delete m_myCalibData;
    m_myCalibData = 0;
  }

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode EmcBhaCalib::initialize(){
  
  MsgStream log(msgSvc(), name());
  log << MSG::INFO << "in initialize()" << endreq;

  m_myCalibData = new EmcBhaCalibData(6240, m_MsgFlag);

  m_calibrationSuccessful = false;  
  
  StatusCode status1;
  
  NTuplePtr nt1(ntupleSvc(),"FILE302/n1");
  if ( nt1 ) m_tuple1 = nt1;
  else {
    m_tuple1=ntupleSvc()->book("FILE302/n1",CLID_ColumnWiseTuple,"Calib");
    if( m_tuple1 ) {
      
         status1 = m_tuple1->addItem ("ixtal",ixtal); 
         status1 = m_tuple1->addItem ("gi0",gi0);
         status1 = m_tuple1->addItem ("calibConst",calibConst);
         status1 = m_tuple1->addItem ("err",err);
         status1 = m_tuple1->addItem ("nhitxtal",nhitxtal);
      
    }
    else    {   // did not manage to book the N tuple....
      log << MSG::ERROR <<"Cannot book N-tuple:" << long(m_tuple1) << endmsg;
      return StatusCode::FAILURE;
    }
  }

  // use EmcCalibConstSvc
  StatusCode scCalib; 
  scCalib = Gaudi::svcLocator() -> service("EmcCalibConstSvc", m_emcCalibConstSvc);
  if( scCalib != StatusCode::SUCCESS){
    log << MSG::ERROR << "can not use EmcCalibConstSvc" << endreq;
  } 
  else {
    std::cout << "Test EmcCalibConstSvc   DigiCalibConst(0)=  " 
              << m_emcCalibConstSvc -> getDigiCalibConst(0) << std::endl;
  }
  
   //init starting values for calibration constants from file or set to 1
  initCalibConst();  

  //digiConstCor();

  return StatusCode::SUCCESS;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode EmcBhaCalib::execute() {
  
  MsgStream log(msgSvc(), name());
  log << MSG::DEBUG << "in execute()" << endreq;   
 //read in accumulated matrix/matrices
  if ( m_readDataFromDB ) {
    m_calibrationSuccessful = readInFromDB();

    log << MSG::INFO << "read in accumulated matrix from DB"<<endreq;

  } else {
    m_calibrationSuccessful = readInFromFile();

    log << MSG::INFO <<"read in accumulated matrix from file"<<endreq; 

  }
  
  if ( m_calibrationSuccessful == true ) {
    
    //ask first if to do a matrix inversion
    if (m_askForMatrixInversion==true) {

      m_calibrationSuccessful = false;
      log << MSG::INFO << " Well, we have " 
          << m_nrXtalsEnoughHits << " crystals whith more "
          << "than " <<  m_dirHitsCut 
          << " direct hits. Do you want do "
          << "a matrix inversion ? [N]" << endreq; 


      m_calibrationSuccessful = true;

    }
   
    if ( m_calibrationSuccessful == true ) {
      
      //write added matrix and vector to files
      if ( m_writeToFile == true) {     
        writeOut();
      }
      //cout <<"writeout"<<endl;
        m_myCalibData->printVec(2);

      //reduce to continious array of only non zeros elements for SLAP
      if ( m_calibrationSuccessful = m_myCalibData->reduce() ) {

        cout<<"m_calibrationSuccessful="<<m_calibrationSuccessful<<endl;

        if ( m_myCalibData->nXtalsHit() != m_myCalibData->nXtals() ){
          log << MSG::INFO << " Reduced number of Xtals for calibration: " 
              << m_myCalibData->nXtalsHit() 
              << endreq;
        }
        cout<<"m_myCalibData->nXtalsHit()="<<m_myCalibData->nXtalsHit()
            <<"m_myCalibData->nXtals()="<<m_myCalibData->nXtals()<<endl;
        m_myCalibData->printVec(10);

        //solve matrix equation
        if ( m_calibrationSuccessful = solveEqua() ) {
          
          for (int i=0; i<m_myCalibData->nXtalsHit(); i++){
            //fill an array of constants for storage in database
            m_calibConstUnred[m_myCalibData->xtalInd(i)] 
              = m_calibConst[i];
            
            // if (m_myCalibData->xtalHitsDir(i)>10)
              // cout<<"Index,calibconst="<<"  "<<i <<"  "<<m_myCalibData->xtalInd(i)
              //  <<"  "<<m_calibConstUnred[m_myCalibData->xtalInd(i)]
              //  <<"  "<<m_calibConstUnred[m_myCalibData->xtalInd(i)]*
              //  m_oldConstants[m_myCalibData->xtalInd(i)]<<endl;
            
          }
          //do validation, fit polynom if necessary, fill CalList         
          prepareConstants();

          //if wanted write constants to plain ASCII file 
          if ( m_writeToFile==true){ 
            writeOutConst();
          }

          ntupleOut();

        }
        
      } else {
        log << MSG::WARNING << " Reduction of the Emc Bhabha calibration matrix FAILED !"
            << endl
            << " Will NOT perform Emc Bhabha calibration !"
            << endreq;
        return( StatusCode::FAILURE);
      }
    }
  } else {
    log << MSG::WARNING << " ERROR in reading in Emc Bhabha calibration data !"
        << endl
        << " Will NOT perform Emc Bhabha calibration !"
        << endreq;
    return( StatusCode::FAILURE);
  }

  return StatusCode::SUCCESS;
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode EmcBhaCalib::finalize() {

  MsgStream log(msgSvc(), name());


  log << MSG::INFO << "in endRun()" << endreq;


  return StatusCode::SUCCESS;
}


void
EmcBhaCalib::help() {
  
  MsgStream log(msgSvc(), name());
  
  log << MSG::INFO<< "Performs the Chi square fit of the accumulated " 
      << endreq;
}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void
EmcBhaCalib::initCalibConst( ) {

  
  MsgStream log(msgSvc(), name());
  
  for ( int i = 0; i< m_myCalibData->nXtals(); i++ ) {
    m_calibConst[i] = 1.;
  }  

  ifstream inConst("EmcCalib.Const");  
  if ( !inConst ) {
    log << MSG::VERBOSE  
        << " No starting values for calibration constants found ! "
        << "(EmcCalib.Const)" 
        << endl
        << "Set them to 1.0 ! " << endreq;
  } else {
    log << MSG::VERBOSE << " Read in starting values of calibration constants !"
        << endreq;
    
    int nConst;
    inConst >> nConst;
    log << MSG::VERBOSE << " Have starting values for " 
        << nConst << " Xtals !" << endl
        << "Set the others to 1.0 ! " << endmsg;
    
    int numberXtal;
    for ( int i = 0; i< nConst; i++ ) {
      inConst >> numberXtal >> m_calibConst[numberXtal];
    }
  }
  
  int nConstEmc;
      
  nConstEmc= m_emcCalibConstSvc -> getDigiCalibConstNo() ;

  if ( nConstEmc!=6240) cout<<"number of calibconst="<< nConstEmc<<endl;
  
  // log << MSG::VERBOSE << " Have starting values for " 
  //   << nConstEmc << " Xtals !" << endl
  //   << "Set the others to 1.0 ! " << endmsg;
  
  for ( int i = 0; i< nConstEmc; i++ ) {
    //cout<<i<<"  "
    //<<m_emcCalibConstSvc->getThetaIndex(i)<<"  "
    //<<m_emcCalibConstSvc->getPhiIndex(i)
    //<<"  "<<m_emcCalibConstSvc->getEmcID(i)<<endl;  
    m_calibConst[i] = m_emcCalibConstSvc -> getDigiCalibConst(i);
    //m_calibConst[i] =5.0;
    m_oldConstants[i]=m_emcCalibConstSvc -> getDigiCalibConst(i);

    // initialize m_absoluteConstants as  m_oldConstants.
    //m_absoluteConstants[i] =m_emcCalibConstSvc -> getDigiCalibConst(i);
    m_absoluteConstants[i] =1.0;
    // initialize m_calibConstUnred as 1.
    m_calibConstUnred[i] =1.0;
  }


}


//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
bool
EmcBhaCalib::solveEqua() {
  
  MsgStream log(msgSvc(), name());
  //-----------------------------------------------------
  //  solve system of equations with SLAP package 
  //-----------------------------------------------------
  log << MSG::VERBOSE << " Solve Matrix equation with method " 
      << m_equationSolver
      << endl
      << "Xtals hit: " << m_myCalibData->nXtalsHit() 
      << endl
      << "Nr of non Zeros: " << m_nrNonZeros << endreq;
  
  //input parameters for SLAP
  long int isym = 1;     //matrix M is symmetric
  long int nsave = 20;
  long int itol = 0;     //type of convergence criterion
  double convCriterion = m_convCrit;    //convergence crtiterion
  long int maxIter = 1000;    //maximum number of iterations
  long int errUnit = 6;     //unit on which to write error estimation
  // at each iteration  

  //working arrays needed by SLAP
  long int lengthDoubleWork;
  double* doubleWork;
  long int lengthIntWork;
  long int* intWork;

  //output parameters of slap
  long int iters=0;           //number of iterations required
  double errorv=1000;        //error estimate
  long int errorFlag=9999;  

  //sparse Ax=b solver.
  //uses the generalized minimum residual method 
  //The preconditioner  is   diagonal scaling.
  if ( m_equationSolver == "GMR" ) {

    cout<<m_equationSolver<<endl;
 
    cout<<"errorFlag="<<errorFlag<<endl;

    //workspaces
    lengthDoubleWork = (1 + m_myCalibData->nXtals()*(nsave+7) 
                         + nsave*(nsave+3))
                         + 50000;
    doubleWork = new double[lengthDoubleWork];
    lengthIntWork = 50;
    intWork = new long int [lengthIntWork];
  
    dsdgmr_ ( &(m_myCalibData->nXtalsHit()),
              m_myCalibData->getVectorR(),
              m_calibConst,
              &m_nrNonZeros,
              m_myCalibData->getMatrixM()->row(),
              m_myCalibData->getMatrixM()->column(),
              m_myCalibData->getMatrixM()->matrix(),
              &isym,
              &nsave,
              &itol,
              &convCriterion,
              &maxIter,
              &iters,
              &errorv,
              &errorFlag,
              &errUnit,
              doubleWork,
              &lengthDoubleWork,
              intWork,
              &lengthIntWork
             );   
    cout<<"errorFlag="<<errorFlag<<endl;

          
    log << MSG::VERBOSE << " Error flag: " <<  errorFlag << endreq;
    if ( errorFlag < 0 ) errorFlag = labs(errorFlag) + 2;
    switch ( errorFlag ) {
    case 0: log << MSG::VERBOSE << " all went well" 
                << endreq; break;
    case 1: log << MSG::ERROR << " insufficient storage allocated for _doubleWork or _intWork" 
                << endreq; break;
    case 2: log << MSG::ERROR << " method failed to reduce norm of current residual" 
                << endreq; break;
    case 3: log << MSG::ERROR << " insufficient length for _doubleWork" 
                << endreq; break;
    case 4: log << MSG::ERROR << " inconsistent _itol and  values" 
                << endreq; break;
    default: log << MSG::ERROR << " unknown flag" << endreq; 
    }
    log << MSG::VERBOSE << " Integer workspace used = " << intWork[8] << endl
        << " Double  workspace used = " << intWork[9] << endreq;
  }

  //Routine to solve a  symmetric positive definite  linear
  //system Ax  =  b  using   the  Preconditioned  Conjugate
  //Gradient   method.    The preconditioner  is   diagonal scaling.
  if ( m_equationSolver == "PCG" ) {
    
    cout<<m_equationSolver<<endl;

    itol = 1; 

    //workspaces
    lengthDoubleWork = 5 *m_myCalibData->nXtals() + 50000;
    doubleWork = new double[lengthDoubleWork];
    lengthIntWork = 50;
    intWork = new long int [lengthIntWork];
       
    dsdcg_( &(m_myCalibData->nXtalsHit()),
            m_myCalibData->getVectorR(),
            m_calibConst,
            &m_nrNonZeros,
            m_myCalibData->getMatrixM()->row(),
            m_myCalibData->getMatrixM()->column(),
            m_myCalibData->getMatrixM()->matrix(),
            &isym,
            &itol,
            &convCriterion,
            &maxIter,
            &iters,
            &errorv,
            &errorFlag,
            &errUnit,
            doubleWork,
            &lengthDoubleWork,
            intWork,
            &lengthIntWork
            );
      
    switch ( errorFlag ) {
    case 0: log << MSG::VERBOSE << "all went well" << endreq; break;
    case 1: log << MSG::ERROR << " insufficient storage allocated for WORK or IWORK" 
                << endreq; break;
    case 2: log << MSG::ERROR << " method failed to to converge in maxIter steps." 
                << endreq; break;
    case 3:log << MSG::ERROR << " Error in user input.  Check input value of _nXtal,_itol." 
               << endreq; break;
    case 4:log << MSG::ERROR << " User error tolerance set too tight. "
               << "Reset to 500.0*D1MACH(3).  Iteration proceeded." 
               << endreq; break;
    case 5:log << MSG::ERROR << " Preconditioning matrix, M,  is not Positive Definite. "
               << endreq; break;
    case 6: log << MSG::ERROR << " Matrix A is not Positive Definite."
                << endreq; break;
    default: log << MSG::ERROR << " unknown flag" << endreq; 
    }
     log << MSG::VERBOSE << " Integer workspace used = " << intWork[9] << endl
                    << "Double  workspace used = " << intWork[10] << endreq;
  }

   log << MSG::VERBOSE << "------ Calibration fit statistics ------- " << endl
                  << "maximum number of iterations =" << maxIter << endl
                  << "        number of iterations =" << iters << endl
                  << "error estimate of error in final solution =" 
                  << errorv << endreq;

  if ( 0 != doubleWork) delete [] doubleWork;     
  if ( 0 != intWork) delete [] intWork;

  if ( errorFlag != 0 ) return false;
  else return true; 

  return true;
}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void
EmcBhaCalib::writeOut() {

  ofstream vectorOut;
  std::string vectorFile = m_fileDir;
  vectorFile += m_fileExt;
  vectorFile += "allCalibVector.dat";
  vectorOut.open(vectorFile.c_str());

  ofstream matrixOut;
  std::string matrixFile = m_fileDir;
  matrixFile += m_fileExt;
  matrixFile += "allCalibMatrix.dat";
  matrixOut.open(matrixFile.c_str());

  MsgStream log(msgSvc(), name());

  log << MSG::VERBOSE << " Write out files " 
      << vectorFile << "   " 
      << matrixFile
      << endreq;

  m_myCalibData->writeOut(matrixOut,vectorOut);

  vectorOut.close();
  matrixOut.close();

}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
void
EmcBhaCalib::writeOutConst() {
  
  ofstream constOut;  

  std::string constFile = m_fileDir;
  constFile += m_fileExt;
  constFile += "CalibConst.dat";

  constOut.open(constFile.c_str());
  
  constOut << "#crystalIndex relative-constant absolute-constant" << endl;

  int chanIndex;
   
  //output constants to file
  for ( int i=0; i < m_myCalibData->nXtalsHit(); i++) {

     chanIndex=m_myCalibData->xtalInd(i);  //xtalInd(i) array of xtal indices
 
    //---- get the new absolute constant ----
    //set it to 0 if we have not enough hits -> we'll keep the old constant
    if ( m_myCalibData->xtalHitsDir(i) < m_dirHitsCut )
      m_absoluteConstants[chanIndex] = m_oldConstants[chanIndex];
    else 
      m_absoluteConstants[chanIndex] =
        m_oldConstants[chanIndex] * m_calibConstUnred[chanIndex];

  }

  //output constants to file
  for ( int i=0; i < m_myCalibData->nXtals(); i++) {

    long Xtal_Index = i;
      if(m_calibConstUnred[Xtal_Index]>0){
        constOut << Xtal_Index << " " 
                 << m_calibConstUnred[Xtal_Index] << " " 
                 << m_oldConstants[Xtal_Index] << " " 
                 << m_absoluteConstants[Xtal_Index]
                 << endl; 
        
        
      }
  }
  constOut.close();

}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
bool
EmcBhaCalib::readInFromFile() {

  MsgStream log(msgSvc(), name());

  log << MSG::INFO << " Read Bhabha calibration data from files !"
                  << endreq;

  std::string runNumberFile = m_fileDir;
  runNumberFile += m_fileExt;
  runNumberFile += m_runNumberFile;

  bool success = false;
  log << MSG::INFO << " Get file names from input file : " 
      << runNumberFile
      << endreq;
  
  std::string vectorFile;
  std::string matrixFile;

  log << MSG::INFO  << "Runnumber    non-zeros    xtals"
      << endreq;   

  ifstream datafile(runNumberFile.c_str());

  //cout<<"datafile="<<runNumberFile.c_str()<<""<<datafile.good()<<endl;

  if (datafile.good() > 0 ) {
    
    while( !datafile.eof() ) {
      
      ifstream vectorIn;
      ifstream matrixIn;
      
      std::string num;
      datafile >> num;
      
      if ( num.length() > 0 ) {
        
        vectorFile = m_fileDir;
        vectorFile += m_fileExt;
        vectorFile += num;
        vectorFile += "CalibVector.dat"; 

        matrixFile = m_fileDir;
        matrixFile += m_fileExt;
        matrixFile += num;
        matrixFile += "CalibMatrix.dat";

        vectorIn.open(vectorFile.c_str());
        matrixIn.open(matrixFile.c_str());

        if ( vectorIn.good() > 0 && matrixIn.good() > 0 ) {
          
          m_myCalibData->readIn(matrixIn,vectorIn);
          
          log << MSG::INFO << num 
              << "          "     
              << m_myCalibData->getMatrixM()->num_nonZeros()
              << "       "
              << m_myCalibData->nXtalsHit() 
              << endreq;

          success = true;

        }
        else {
          if ( !vectorIn )
            log << MSG::ERROR << " ERROR: Vector input file "
                << vectorFile
                <<  " doesn't exist !" << endreq;
          if ( !matrixIn )
            log << MSG::ERROR << " ERROR: Matrix input file "
                << matrixFile
                << " doesn't exist !" << endreq;
        }
        vectorIn.close();
        matrixIn.close();
      }
    }
  }

  if ( success == true) {

    for (int i=0; i < m_myCalibData->nXtals(); i++) {

      if ( m_myCalibData->xtalHitsDir(i) > m_dirHitsCut ) 
        {
          m_nrXtalsEnoughHits++;
        }
    }
    m_nrNonZeros = m_myCalibData->getMatrixM()->num_nonZeros();
    log << MSG::INFO<< " Final matrix :  "
        << "Number of non zero elements: " << m_nrNonZeros 
        << "       "
        << m_myCalibData->nXtalsHit()
        << endreq;

  }


  return success;
}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
bool
EmcBhaCalib::readInFromDB() {

  bool success = true;

  return success;
}

//* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
bool
EmcBhaCalib::prepareConstants() {

  bool successfull=false;

  //the old constant
  //float theCalConst = 1.0;
  int chanIndex;

  for ( int i = 0; i< m_myCalibData->nXtalsHit(); i++ ) {

    chanIndex=m_myCalibData->xtalInd(i);  //xtalInd(i) array of xtal indices
   
    // m_calibConstUnred[chanIndex]=1;
    //if (chanIndex==689) m_calibConstUnred[chanIndex]=1.47;
    //if (chanIndex==689)
    //  m_oldConstants[chanIndex]=m_oldConstants[chanIndex]*1.47;
    //---- get the new absolute constant ----
    //set it to 0 if we have not enough hits -> we'll keep the old constant
    if ( m_myCalibData->xtalHitsDir(i) < m_dirHitsCut )
      m_absoluteConstants[chanIndex] = m_oldConstants[chanIndex];
    else 
      m_absoluteConstants[chanIndex] =
        m_oldConstants[chanIndex] * m_calibConstUnred[chanIndex];
  }


  double DigiConst[6240];
  int IxtalNumber[6240];

  for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) {

    DigiConst[ind]=m_absoluteConstants[ind];
    // cout<<"ind="<<ind<<"\t"<< DigiConst[ind]<<endl;
  }

  int ixtal1,ixtal2,ixtal3;
  ixtal1=m_emcCalibConstSvc->getIndex(1,20,26);
  ixtal2=m_emcCalibConstSvc->getIndex(1,23,26);
  ixtal3=m_emcCalibConstSvc->getIndex(1,24,26);
  //cout<<ixtal1<<"   "<<ixtal2<<"   "<<ixtal3<<"   "<<endl;
  for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) {

    IxtalNumber[ind]=-1;
    /*    
    if (ind==ixtal1) IxtalNumber[ind]=ixtal3;
    if (ind==ixtal2) IxtalNumber[ind]=ixtal1;
    if (ind==ixtal3) IxtalNumber[ind]=ixtal2;
    */
  }

 TFile fconst("EmcCalibConst.root", "recreate");

 //define tree fill the absolute digicalibconsts into the root-file
  TTree* constr= new TTree("DigiCalibConst","DigiCalibConst");
  constr->Branch("DigiCalibConst",DigiConst,"DigiConst[6240]/D"); 
  constr->Branch("IxtalNumber",IxtalNumber,"IxtalNumber[6240]/I"); 

  constr->Fill();
  /* 

  double aa[6240];
  int Iaa[6240];
  constr->SetBranchAddress("DigiCalibConst", aa);
  constr->SetBranchAddress("IxtalNumber", Iaa);
  constr->GetEntry(0);
  cout<<Iaa[10]<<endl;  
  cout<<aa[10]<<endl;

  cout<<constr->GetEntry(0)<<endl;
  
  */

  constr->Write();
 
  delete constr;
  
  fconst.Close();

  // end tree
   
  successfull=true;
  return successfull;

}

void 
EmcBhaCalib::ntupleOut(){


  for (int i=0; i < m_myCalibData->nXtalsHit(); i++) {
        int xtalIndex=m_myCalibData->xtalInd(i);

        gi0 = m_oldConstants[xtalIndex];

        if (gi0<0) cout<<"gi0="<<gi0<<endl;
        err = (m_calibConstUnred[xtalIndex]-gi0)/gi0;
        calibConst=m_calibConstUnred[xtalIndex];
        ixtal=xtalIndex;
        nhitxtal=m_myCalibData->xtalHitsDir(i); 
        m_tuple1->write();   
  }

}


void
EmcBhaCalib::printInfo(std::string fileName) {

  ofstream xtalHitsDirOut;
  std::string xtalHitsDirFile = m_fileDir;
  xtalHitsDirFile += m_fileExt;
  xtalHitsDirFile += fileName;
  xtalHitsDirOut.open(xtalHitsDirFile.c_str());

  //nXtalsHit() is number of xtals hit
  for (int i=0; i < m_myCalibData->nXtalsHit(); i++) {
    if ( m_myCalibData->xtalHitsDir(i) > m_dirHitsCut ) 
      {
        int chanindex=m_myCalibData->xtalInd(i);
        xtalHitsDirOut<<chanindex<<endl;
      }
  }
  xtalHitsDirOut.close();

}

void 
EmcBhaCalib::digiConstCor(){

  // ofstream Fuout;
  //std::string Fuoutfile; 

  //Fuoutfile="emccalibconst.txt";
  //Fuout.open(Fuoutfile.c_str());

  double digiConst[6240];

  for(int ind=0; ind<m_myCalibData->nXtals(); ind++ ) {

    digiConst[ind]=m_oldConstants[ind];
    // Fuout<<m_oldConstants[ind]<<endl;
  }

  // Fuout.close();

  ifstream ExpEneIn;
  std::string ExpEneFile; 
  ExpEneFile = "cor.conf";
  ExpEneIn.open(ExpEneFile.c_str());

  double ExpEne[56];
  double Exp_barrel[44],Exp_east[6],Exp_west[6];

  for(int i=0;i<56;i++) {

    ExpEneIn>>ExpEne[i];
    ExpEne[i]=ExpEne[i]*1.843;
    if (i<6) Exp_east[i]=ExpEne[i];
    if (i>=6&&i<50) Exp_barrel[i-6]=ExpEne[i];
    if (i>=50) Exp_west[55-i]=ExpEne[i];

    //cout<<ExpEne[i]<<endl;
  }
  ExpEneIn.close();

  ifstream EDepEneIn;
  std::string EDepEneFile = "allxtal.dat";
  EDepEneIn.open(EDepEneFile.c_str());

  double CorDigiConst[6240];
  int ipart,ithe,iphi;
  double peak,sigma;
  double coeff;

  for(int i=0;i<6240;i++) {

    EDepEneIn>>ipart>>ithe>>iphi>>peak>>sigma;
    int ix=m_emcCalibConstSvc->getIndex(ipart,ithe,iphi);

    if (ipart==0) coeff=peak/Exp_east[ithe];
    if (ipart==1) coeff=peak/Exp_barrel[ithe];
    if (ipart==2) coeff=peak/Exp_west[ithe];
    cout<<coeff<<endl;
    //CorDigiConst[ix]=digiConst[ix]/coeff;
    CorDigiConst[ix]=coeff;
  }
  EDepEneIn.close();


  //define tree fill the absolute digicalibconsts into the root-file

  TTree* constr= new TTree("DigiCalibConst","DigiCalibConst");
    constr->Branch("DigiCalibConst",CorDigiConst,"CorDigiConst[6240]/D");
  constr->Fill();

  double aa[6240];
   constr->SetBranchAddress("DigiCalibConst", aa);
   constr->GetEntry(0);

  TFile fconst("EmcCalibConstCor.root", "recreate");
  constr->Write();

   fconst.Close();
   delete constr;


}

---