00001 #include "MdcCalibAlg/XtInteMdcCalib.h"
00002
00003 #include "GaudiKernel/MsgStream.h"
00004 #include "GaudiKernel/IMessageSvc.h"
00005 #include "GaudiKernel/StatusCode.h"
00006 #include "GaudiKernel/ISvcLocator.h"
00007 #include "GaudiKernel/Bootstrap.h"
00008
00009 #include <iostream>
00010 #include <fstream>
00011 #include <iomanip>
00012 #include <cstring>
00013 #include <cmath>
00014
00015 #include "TF1.h"
00016
00017 using namespace std;
00018
00019 XtInteMdcCalib::XtInteMdcCalib(){
00020 m_fgIni = false;
00021 m_nMaxGrPoint = 5000;
00022 for(int lay=0; lay<MdcCalNLayer; lay++){
00023 m_tbinWid[lay][0] = 5.0;
00024 m_tbinWid[lay][1] = 10.0;
00025 m_tbinWid[lay][2] = 20.0;
00026
00027 m_tbinLim[lay][0] = -10.0;
00028 m_tbinLim[lay][1] = 30.0;
00029 if(lay < 8) m_tbinLim[lay][2] = 210.0;
00030 else m_tbinLim[lay][2] = 350.0;
00031 m_tbinLim[lay][3] = 990.0;
00032 }
00033 }
00034
00035 XtInteMdcCalib::~XtInteMdcCalib(){
00036 }
00037
00038 void XtInteMdcCalib::clear(){
00039 for(int lay=0; lay<MdcCalNLayer; lay++){
00040 for(int iEntr=0; iEntr<NENTR; iEntr++){
00041 for(int lr=0; lr<2; lr++){
00042 delete m_pfNear[lay][iEntr][lr];
00043 delete m_pfMid[lay][iEntr][lr];
00044 delete m_pfFar[lay][iEntr][lr];
00045 delete m_grXt[lay][iEntr][lr];
00046 }
00047 }
00048 }
00049 delete m_fdPf;
00050 MdcCalib::clear();
00051 }
00052
00053 void XtInteMdcCalib::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00054 IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc) {
00055 IMessageSvc* msgSvc;
00056 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00057 MsgStream log(msgSvc, "XtInteMdcCalib");
00058 log << MSG::INFO << "XtInteMdcCalib::initialize()" << endreq;
00059
00060 m_hlist = hlist;
00061 m_mdcGeomSvc = mdcGeomSvc;
00062 m_mdcFunSvc = mdcFunSvc;
00063 m_mdcUtilitySvc = mdcUtilitySvc;
00064
00065 MdcCalib::initialize(m_hlist, m_mdcGeomSvc, m_mdcFunSvc, m_mdcUtilitySvc);
00066
00067 m_fdPf = new TFolder("fdProfile", "fdProfile");
00068 m_hlist -> Add(m_fdPf);
00069
00070 char hname[200];
00071 for(int lay=0; lay<MdcCalNLayer; lay++){
00072 for(int iEntr=0; iEntr<NENTR; iEntr++){
00073 for(int lr=0; lr<2; lr++){
00074 sprintf(hname, "xt%02d_%02d_%d_gr", lay, iEntr, lr);
00075 m_grXt[lay][iEntr][lr] = new TGraph();
00076 m_grXt[lay][iEntr][lr]->SetName(hname);
00077 m_fdPf->Add(m_grXt[lay][iEntr][lr]);
00078
00079 int xbinN = (int)((m_tbinLim[lay][1] - m_tbinLim[lay][0])/m_tbinWid[lay][0] + 0.5);
00080 sprintf(hname, "xt%02d_%02d_%d_near", lay, iEntr, lr);
00081 m_pfNear[lay][iEntr][lr] = new TProfile(hname, hname, xbinN, m_tbinLim[lay][0], m_tbinLim[lay][1]);
00082 m_fdPf->Add(m_pfNear[lay][iEntr][lr]);
00083
00084 int xbinM = (int)((m_tbinLim[lay][2] - m_tbinLim[lay][1])/m_tbinWid[lay][1] + 0.5);
00085 sprintf(hname, "xt%02d_%02d_%d_mid", lay, iEntr, lr);
00086 m_pfMid[lay][iEntr][lr] = new TProfile(hname, hname, xbinM, m_tbinLim[lay][1], m_tbinLim[lay][2]);
00087 m_fdPf->Add(m_pfMid[lay][iEntr][lr]);
00088
00089 int xbinF = (int)((m_tbinLim[lay][3] - m_tbinLim[lay][2])/m_tbinWid[lay][2] + 0.5);
00090 sprintf(hname, "xt%02d_%02d_%d_far", lay, iEntr, lr);
00091 m_pfFar[lay][iEntr][lr] = new TProfile(hname, hname, xbinF, m_tbinLim[lay][2], m_tbinLim[lay][3]);
00092 m_fdPf->Add(m_pfFar[lay][iEntr][lr]);
00093
00094
00095
00096
00097
00098
00099 }
00100 }
00101 }
00102 }
00103
00104 int XtInteMdcCalib::fillHist(MdcCalEvent* event){
00105 IMessageSvc* msgSvc;
00106 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00107 MsgStream log(msgSvc, "XtInteMdcCalib");
00108 log << MSG::DEBUG << "XtInteMdcCalib::fillHist()" << endreq;
00109
00110 MdcCalib::fillHist(event);
00111
00112
00113 bool esCutFg = event->getEsCutFlag();
00114 if( ! esCutFg ) return -1;
00115
00116 int i;
00117 int k;
00118 int lay;
00119 int iLR;
00120 int iEntr;
00121
00122 double dr;
00123 double dz;
00124 double doca;
00125 double tdr;
00126 double resi;
00127 double entrance;
00128
00129 int nhitlay;
00130 bool fgHitLay[MdcCalNLayer];
00131 bool fgTrk;
00132
00133 if(! m_fgIni){
00134 for(lay=0; lay<MdcCalNLayer; lay++){
00135 if(lay < 8) m_docaMax[lay] = m_param.maxDocaInner;
00136 else m_docaMax[lay] = m_param.maxDocaOuter;
00137 }
00138 m_fgIni = true;
00139 }
00140
00141 MdcCalRecTrk* rectrk;
00142 MdcCalRecHit* rechit;
00143
00144 int nhit;
00145 int ntrk = event -> getNTrk();
00146 for(i=0; i<ntrk; i++){
00147 fgTrk = true;
00148 rectrk = event->getRecTrk(i);
00149 nhit = rectrk -> getNHits();
00150
00151
00152 dr = rectrk->getDr();
00153 if(fabs(dr) > m_param.drCut) continue;
00154
00155
00156 dz = rectrk->getDz();
00157 if(fabs(dz) > m_param.dzCut) continue;
00158
00159 for(lay=0; lay<MdcCalNLayer; lay++){
00160 fgHitLay[lay] = false;
00161 }
00162
00163 for(k=0; k<nhit; k++){
00164 rechit = rectrk -> getRecHit(k);
00165 lay = rechit -> getLayid();
00166 doca = rechit -> getDocaInc();
00167 resi = rechit -> getResiInc();
00168 fgHitLay[lay] = true;
00169
00170
00171
00172
00173
00174 }
00175 if(! fgTrk) continue;
00176
00177 nhitlay = 0;
00178 for(lay=0; lay<MdcCalNLayer; lay++){
00179 if(fgHitLay[lay]) nhitlay++;
00180 }
00181 if(nhitlay < m_param.nHitLayCut) continue;
00182
00183 for(k=0; k<nhit; k++){
00184 rechit = rectrk -> getRecHit(k);
00185 lay = rechit -> getLayid();
00186 doca = rechit -> getDocaInc();
00187 resi = rechit -> getResiInc();
00188 iLR = rechit -> getLR();
00189 entrance = rechit -> getEntra();
00190 tdr = rechit -> getTdrift();
00191
00192 if( (fabs(doca) > m_docaMax[lay]) ||
00193 (fabs(resi) > m_param.resiCut[lay]) ){
00194 continue;
00195 }
00196
00197 if(0 == lay){
00198 if( ! fgHitLay[1] ) continue;
00199 } else if(42 == lay){
00200 if( ! fgHitLay[41] ) continue;
00201 } else{
00202 if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ) continue;
00203 }
00204
00205 iEntr = m_mdcFunSvc -> getXtEntrIndex(entrance);
00206 int npoint = m_grXt[lay][iEntr][iLR]->GetN();
00207 if(npoint < m_nMaxGrPoint) m_grXt[lay][iEntr][iLR]->SetPoint(npoint, tdr, fabs(doca));
00208
00209 if(tdr<m_tbinLim[lay][0]) continue;
00210 else if(tdr<m_tbinLim[lay][1]) m_pfNear[lay][iEntr][iLR]->Fill(tdr, fabs(doca));
00211 else if(tdr<m_tbinLim[lay][2]) m_pfMid[lay][iEntr][iLR]->Fill(tdr, fabs(doca));
00212 else if(tdr<m_tbinLim[lay][3]) m_pfFar[lay][iEntr][iLR]->Fill(tdr, fabs(doca));
00213 }
00214 }
00215 return 1;
00216 }
00217
00218 int XtInteMdcCalib::updateConst(MdcCalibConst* calconst){
00219 IMessageSvc* msgSvc;
00220 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00221 MsgStream log(msgSvc, "XtInteMdcCalib");
00222 log << MSG::INFO << "XtInteMdcCalib::updateConst()" << endreq;
00223
00224 MdcCalib::updateConst(calconst);
00225 }
00226