/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/MdcCalibAlg/MdcCalibAlg-00-09-02/src/IniMdcCalib.cxx

Go to the documentation of this file.

#include "MdcCalibAlg/IniMdcCalib.h"
 
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/Bootstrap.h"
 
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/PropertyMgr.h"

#include "EventModel/EventHeader.h"
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/Event.h"
#include "MdcRawEvent/MdcDigi.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "TStyle.h"

#include <iostream>
#include <fstream>
#include <iomanip>
#include <string>
#include <cstring>
#include <cmath>

#include "TF1.h"

using namespace Event;

IniMdcCalib::IniMdcCalib(){
}

IniMdcCalib::~IniMdcCalib(){
}

void IniMdcCalib::clear(){
     int iEs;
     for(int lay=0; lay<MdcCalNLayer; lay++){
          delete m_hlaymapT[lay];
          delete m_htdc[lay];
          delete m_htraw[lay];
          for(iEs=0; iEs<m_param.nEsFlag; iEs++){
               delete m_htdcTes[lay][iEs];
               delete m_htrawTes[lay][iEs];
          }

          delete m_hlaymapQ[lay];
          delete m_hqraw[lay];
     }

     for(int wir=0; wir<MdcCalTotCell; wir++){
          delete m_htrawCel[wir];
          delete m_hqrawCel[wir];
     }

     delete m_hLayerHitmapT;
     delete m_hWireHitMapT;

     delete m_hLayerHitmapQ;
     delete m_hWireHitMapQ;
     for(iEs=0; iEs<m_param.nEsFlag; iEs++) delete m_hTes[iEs];
     delete m_hTesAllFlag;
     delete m_hTesAll;
     delete m_hTesCal;
     delete m_hTesFlag;

     delete m_fdcom;
     delete m_fdTmap;
     delete m_fdTraw;
     delete m_fdTrawCel;
     delete m_fdTrawTes;
     delete m_fdQmap;
     delete m_fdQraw;
     delete m_fdQrawCel;
}

void IniMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                             IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "IniMdcCalib");
     log << MSG::INFO << "IniMdcCalib::initialize()" << endreq;

     m_hlist = hlist;
     m_mdcGeomSvc = mdcGeomSvc;
     m_mdcFunSvc = mdcFunSvc;
     m_mdcUtilitySvc = mdcUtilitySvc;

     int lay;
     int cel;
     int wir;
     int ncel;
     char hname[200];

     m_nWire = m_mdcGeomSvc -> getWireSize();
     m_nLayer = m_mdcGeomSvc -> getLayerSize();

     m_fdcom = new TFolder("Common", "Common");
     m_hlist->Add(m_fdcom);

     m_fdTmap = new TFolder("Thitmap", "Thitmap");
     m_hlist->Add(m_fdTmap);

     m_fdTraw = new TFolder("Traw", "Traw");
     m_hlist->Add(m_fdTraw);

     m_fdTrawCel = new TFolder("TrawCell", "TrawCell");
     m_hlist->Add(m_fdTrawCel);

     m_fdTrawTes = new TFolder("TrawTes", "TrawTes");
     m_hlist->Add(m_fdTrawTes);

     m_fdQmap = new TFolder("Qhitmap", "Qhitmap");
     m_hlist->Add(m_fdQmap);

     m_fdQraw = new TFolder("Qraw", "Qraw");
     m_hlist->Add(m_fdQraw);

     m_fdQrawCel = new TFolder("QrawCell", "QrawCell");
     m_hlist->Add(m_fdQrawCel);

     m_hLayerHitmapT = new TH1F("T_Hitmap_Layer", "", 43, -0.5, 42.5);
     m_fdcom->Add(m_hLayerHitmapT);
                
     m_hWireHitMapT = new TH1F("T_Hitmap_Wire", "", 6796, -0.5, 6795.5);
     m_fdcom->Add(m_hWireHitMapT);

     m_hLayerHitmapQ = new TH1F("Q_Hitmap_Layer", "", 43, -0.5, 42.5);
     m_fdcom->Add(m_hLayerHitmapQ);
                
     m_hWireHitMapQ = new TH1F("Q_Hitmap_Wire", "", 6796, -0.5, 6795.5);
     m_fdcom->Add(m_hWireHitMapQ);

     int iEs;
     for(iEs=0; iEs<m_param.nEsFlag; iEs++){
          sprintf(hname, "Tes_%d", m_param.esFlag[iEs]);
          m_hTes[iEs] = new TH1F(hname, "", 750, 0, 1500);
          m_fdcom->Add(m_hTes[iEs]);
     }

     m_hTesAllFlag = new TH1F("TesAllFlag", "", 300, -0.5, 299.5);
     m_fdcom -> Add(m_hTesAllFlag);

     m_hTesAll = new TH1F("TesAll", "", 750, 0, 1500);
     m_fdcom->Add(m_hTesAll);

     m_hTesCal = new TH1F("TesCal", "", 750, 0, 1500);
     m_fdcom->Add(m_hTesCal);

     m_hTesFlag = new TH1F("Tes_Flag", "", 300, -0.5, 299.5);
     m_fdcom->Add(m_hTesFlag);

     for(lay=0; lay<m_nLayer; lay++){
          ncel = m_mdcGeomSvc->Layer(lay)->NCell();

          sprintf(hname, "T_hitmap_Lay%02d", lay);
          m_hlaymapT[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
          m_fdTmap -> Add(m_hlaymapT[lay]);

          sprintf(hname, "TDC_Lay%02d", lay);
          m_htdc[lay] = new TH1F(hname, "", 800, 0, 20000);
          m_fdTraw -> Add(m_htdc[lay]);

          sprintf(hname, "Traw_Lay%02d", lay);
          m_htraw[lay] = new TH1F(hname, "", 500, 0, 1000);
          m_fdTraw -> Add(m_htraw[lay]);

          for(iEs=0; iEs<m_param.nEsFlag; iEs++){
               sprintf(hname, "TDC_Lay%02d_Tes%d", lay, m_param.esFlag[iEs]);
               m_htdcTes[lay][iEs] = new TH1F(hname, "", 800, 0, 20000);
               m_fdTrawTes -> Add(m_htdcTes[lay][iEs]);

               sprintf(hname, "Traw_Lay%02d_Tes%d", lay, m_param.esFlag[iEs]);
               m_htrawTes[lay][iEs] = new TH1F(hname, "", 300, 0, 600);
               m_fdTrawTes -> Add(m_htrawTes[lay][iEs]);
          }

          sprintf(hname, "Q_hitmap_Lay%02d", lay);
          m_hlaymapQ[lay] = new TH1F(hname, "", ncel, -0.5, (float)ncel-0.5);
          m_fdQmap -> Add(m_hlaymapQ[lay]);

          sprintf(hname, "Qraw_Lay%02d", lay);
          m_hqraw[lay] = new TH1F(hname, "", 2000, 0, 4000);
          m_fdQraw -> Add(m_hqraw[lay]);
     }

     for(wir=0; wir<MdcCalTotCell; wir++){
          lay = m_mdcGeomSvc -> Wire(wir) -> Layer();
          cel = m_mdcGeomSvc -> Wire(wir) -> Cell();

          sprintf(hname, "Traw_%02d_%03d_%04d", lay, cel, wir);
          m_htrawCel[wir] = new TH1F(hname, "", 300, 0, 600);
          m_fdTrawCel -> Add(m_htrawCel[wir]);

          sprintf(hname, "Qraw_%02d_%03d_%04d", lay, cel, wir);
          m_hqrawCel[wir] = new TH1F(hname, "", 2000, 0, 4000);
          m_fdQrawCel -> Add(m_hqrawCel[wir]);
     }
}

int IniMdcCalib::fillHist(MdcCalEvent* event){
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "IniMdcCalib");
     log << MSG::DEBUG << "IniMdcCalib::fillHist()" << endreq;

     int lay;
     int cel;
     int wir;
     int digiId;
     unsigned fgOverFlow;
     double tdc;
     double traw;
     double adc;
     double qraw;
     Identifier id; 

     IDataProviderSvc* eventSvc = NULL;
     Gaudi::svcLocator()->service("EventDataSvc", eventSvc);

     SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
     if (!eventHeader) {
          log << MSG::FATAL << "Could not find Event Header" << endreq;
          return( StatusCode::FAILURE);
     }
     int iRun = eventHeader->runNumber();
     int iEvt = eventHeader->eventNumber();

     // get EsTimeCol
     double tes = -9999.0;
     int esTimeflag = -1;
     SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc,"/Event/Recon/RecEsTimeCol");
     if( (!aevtimeCol) || (aevtimeCol->size()==0) ){
          tes = -9999.0;
          esTimeflag = -1;
     }else{
          RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
          for(; iter_evt!=aevtimeCol->end(); iter_evt++){
               tes = (*iter_evt)->getTest();
               esTimeflag = (*iter_evt)->getStat();
//             cout << "tes " << tes << endl;
          }
     }
     m_hTesAllFlag->Fill(esTimeflag);
     m_hTesAll->Fill(tes);
     m_hTesFlag->Fill(esTimeflag);
     int nTesFlag = -1;
     for(int iEs=0; iEs<m_param.nEsFlag; iEs++){
          if(esTimeflag == m_param.esFlag[iEs]){
               m_hTes[iEs]->Fill(tes);
               nTesFlag = iEs;
               break;
          }
     }
     bool fgFillTes = false;
     if( (nTesFlag > -1) && (tes > m_param.tesMin) && (tes < m_param.tesMax) ){
          m_hTesCal->Fill(tes);
          fgFillTes = true;
     }

//      cout << setw(10) << iRun << setw(10) << iEvt << setw(10) << tes << endl;

     // retrieve Mdc digi
     SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc,"/Event/Digi/MdcDigiCol");
     if (!mdcDigiCol) {
          log << MSG::FATAL << "Could not find event" << endreq;
     }

     MdcDigiCol::iterator iter = mdcDigiCol->begin();
     digiId = 0;
     for(; iter != mdcDigiCol->end(); iter++, digiId++) {
          MdcDigi *aDigi = (*iter);
          id = (aDigi)->identify();

          lay = MdcID::layer(id);
          cel = MdcID::wire(id);
          wir = m_mdcGeomSvc->Wire(lay, cel)->Id();

          tdc = (aDigi) -> getTimeChannel();
          adc = (aDigi) -> getChargeChannel();
          fgOverFlow = (aDigi) -> getOverflow();

          if ( ((fgOverFlow & 1) !=0 ) || ((fgOverFlow & 12) != 0) ||
               (aDigi->getTimeChannel() == 0x7FFFFFFF) ||
               (aDigi->getChargeChannel() == 0x7FFFFFFF) ) continue;

          traw = tdc * MdcCalTdcCnv;
          qraw = adc * MdcCalAdcCnv;

          m_hLayerHitmapT -> Fill(lay);
          m_hWireHitMapT -> Fill(wir);
          m_hlaymapT[lay] -> Fill(cel);

          m_hLayerHitmapQ -> Fill(lay);
          m_hWireHitMapQ -> Fill(wir);
          m_hlaymapQ[lay] -> Fill(cel);

          if(fgFillTes){
               traw -= tes;
               traw += m_param.timeShift;

               m_htdc[lay] -> Fill(tdc);
               m_htraw[lay] -> Fill(traw);
               m_hqraw[lay] -> Fill(qraw);

               m_htdcTes[lay][nTesFlag] -> Fill(tdc);
               m_htrawTes[lay][nTesFlag] -> Fill(traw);

               m_htrawCel[wir] -> Fill(traw);
               m_hqrawCel[wir] -> Fill(qraw);
          }
     }
     return 1;
}

int IniMdcCalib::updateConst(MdcCalibConst* calconst){
     IMessageSvc* msgSvc;
     Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
     MsgStream log(msgSvc, "IniMdcCalib");
     log << MSG::DEBUG << "IniMdcCalib::updateConst()" << endreq;

     int lay;
     int wir;
     double t0Fit[MdcCalNLayer];
     double t0Cal[MdcCalNLayer];
     double tmax[MdcCalNLayer];
     double initT0 = m_param.initT0;

     int fitTminFg[MdcCalNLayer];
     int fitTmaxFg[MdcCalNLayer];
     double chisq;
     double ndf;
     double chindfTmin[MdcCalNLayer];
     double chindfTmax[MdcCalNLayer];
     char funname[200];

     // fit Tmin
     TF1* ftmin[MdcCalNLayer];
//      sprintf(funname, "ftmin", lay);
//      TF1* ftmin = new TF1(funname, funTmin, 0, 150, 6);
     for(lay=0; lay<m_nLayer; lay++){
          fitTminFg[lay] = 0;
          chindfTmin[lay] = -1;
          sprintf(funname, "ftmin%02d", lay);
          ftmin[lay] = new TF1(funname, funTmin, 0, 150, 6);

          if(1 == m_param.fgCalib[lay]){
               Stat_t nEntryTot = 0;
               for(int ibin=1; ibin<=25; ibin++){
                    Stat_t entry = m_htraw[lay]->GetBinContent(ibin);
                    nEntryTot += entry;
               }
               double c0Ini = (double)nEntryTot / 25.0;
               double c1Ini = (m_htraw[lay]->GetMaximum()) - c0Ini;

               ftmin[lay] -> SetParameter(0, c0Ini);
               ftmin[lay] -> SetParameter(1, c1Ini);
               ftmin[lay] -> SetParameter(2, 0);
               ftmin[lay] -> SetParameter(4, initT0);
               ftmin[lay] -> SetParameter(5, 2);

               m_htraw[lay] -> Fit(funname, "Q", "",
                                   m_param.tminFitRange[lay][0],
                                   m_param.tminFitRange[lay][1]);
               chisq = ftmin[lay]->GetChisquare();
               gStyle -> SetOptFit(11);
               ndf = ftmin[lay]->GetNDF();
               chindfTmin[lay] = chisq / ndf;
               if(chindfTmin[lay] < m_param.tminFitChindf){
                    fitTminFg[lay] = 1;
                    t0Fit[lay] = ftmin[lay]->GetParameter(4);
                    t0Fit[lay] += m_param.t0Shift;
                    t0Cal[lay] = t0Fit[lay] - m_param.timeShift;
               }
          }

          if(0 == fitTminFg[lay]){
               wir = m_mdcGeomSvc->Wire(lay, 0)->Id();
               t0Cal[lay] = calconst->getT0(wir);
               t0Fit[lay] = t0Cal[lay] + m_param.timeShift;
          }
     }

     // fit Tmax
     TF1* ftmax[MdcCalNLayer];
     for(lay=0; lay<m_nLayer; lay++){
          fitTmaxFg[lay] = 0;
          chindfTmax[lay] = -1;
          sprintf(funname, "ftmax%02d", lay);
          ftmax[lay] = new TF1(funname, funTmax, 250, 500, 4);

          if(1 == m_param.fgCalib[lay]){
               ftmax[lay] -> SetParameter(2, m_param.initTm[lay]);
               ftmax[lay] -> SetParameter(3, 10);

               m_htraw[lay] -> Fit(funname, "Q+", "",
                                   m_param.tmaxFitRange[lay][0],
                                   m_param.tmaxFitRange[lay][1]);
               gStyle -> SetOptFit(11);
               chisq = ftmax[lay]->GetChisquare();
               ndf = ftmax[lay]->GetNDF();
               chindfTmax[lay] = chisq / ndf;
               if(chindfTmax[lay] < m_param.tmaxFitChindf){
                    fitTmaxFg[lay] = 1;
                    tmax[lay] = ftmax[lay]->GetParameter(2);
               }
          }

          if(0 == fitTmaxFg[lay]){
               tmax[lay] = (calconst->getXtpar(lay, 0, 0, 6)) + t0Fit[lay];
          }
     }

     // output for check
     ofstream ft0("iniT0.dat");
     for(lay=0; lay<m_nLayer; lay++){
          ft0 << setw(5) << lay << setw(3) << fitTminFg[lay]
              << setw(12) << t0Cal[lay] << setw(12) << t0Fit[lay]
              << setw(12) << chindfTmin[lay] << setw(5) << fitTmaxFg[lay]
              << setw(12) << tmax[lay] << setw(12) << tmax[lay] - t0Fit[lay]
              << setw(12) << chindfTmax[lay] << endl;
     }
     ft0.close();
     cout << "iniT0.dat was written." << endl;

     // set T0
     int i;
     int nwire = m_mdcGeomSvc -> getWireSize();
     for(i=0; i<nwire; i++){
          lay = m_mdcGeomSvc -> Wire(i) -> Layer();
          if(1 == m_param.fgCalib[lay]){
               calconst -> resetT0(i, t0Cal[lay]);
               calconst -> resetDelT0(i, 0.0);
          }
     }

     if(0 == m_param.fgIniCalConst){
          // set X-T
          int lr;
          int iEntr;
          int ord;
          double xtpar;
          double xtini[8] = {0, 0.03, 0, 0, 0, 0, 999.9, 0};
          for(lay=0; lay<m_nLayer; lay++){
               if(1 != m_param.fgCalib[lay]) continue;

               for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
                    for(lr=0; lr<MdcCalLR; lr++){
                         for(ord=0; ord<MdcCalXtNPars; ord++){
                              if(6 == ord){
                                   xtpar = tmax[lay] - t0Fit[lay];
                              } else{
                                   xtpar = xtini[ord];
                              }
                              calconst -> resetXtpar(lay, iEntr, lr, ord, xtpar);
                         }
                    }
               }
          }

          // set Q-T
          for(lay=0; lay<m_nLayer; lay++){
               calconst -> resetQtpar0(lay, 0.0);
               calconst -> resetQtpar1(lay, 0.0);
          }

          // set S-D
          int bin;
          double sdpar = m_param.initSigma;     // mm
          for(lay=0; lay<m_nLayer; lay++){
               for(iEntr=0; iEntr<MdcCalNENTRSD; iEntr++){
                    for(lr=0; lr<2; lr++){
                         for(bin=0; bin<MdcCalSdNBIN; bin++){
                              calconst -> resetSdpar(lay, iEntr, lr, bin, sdpar);
                         }
                    }
               }
          }
     } else if(2 == m_param.fgIniCalConst){
          int lr;
          int iEntr;
          double xtpar;
          for(lay=0; lay<m_nLayer; lay++){
               for(iEntr=0; iEntr<MdcCalNENTRXT; iEntr++){
                    for(lr=0; lr<MdcCalLR; lr++){
                         xtpar = tmax[lay] - t0Fit[lay];
                         calconst -> resetXtpar(lay, iEntr, lr, 6, xtpar);
                    }
               }
          }
     }

//      delete ftmin;
     for(lay=0; lay<MdcCalNLayer; lay++){
          delete ftmin[lay];
          delete ftmax[lay];
     }
     return 1;
}

Double_t IniMdcCalib::funTmin(Double_t* x, Double_t* par){
     Double_t fitval;
     fitval = par[0] + par[1]*exp( -par[2]*(x[0]-par[3]) ) /
          ( 1 + exp( -(x[0]-par[4])/par[5] ));
     return fitval;
}

Double_t IniMdcCalib::funTmax(Double_t* x, Double_t* par){
     Double_t fitval;
     fitval = par[0] + par[1] / (1 + exp((x[0]-par[2])/par[3]));
     return fitval;
}

---