00001 #include "GaudiKernel/Bootstrap.h"
00002 #include "GaudiKernel/PropertyMgr.h"
00003 #include "GaudiKernel/StatusCode.h"
00004
00005 #include <fstream>
00006 #include <string>
00007 #include "TFile.h"
00008 #include "TTree.h"
00009
00010 #include "DedxCalibAlg/DedxCalib.h"
00011
00012 DedxCalib::DedxCalib(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator)
00013 {
00014 declareProperty("ParticleType",ParticleFlag=-1);
00015 declareProperty("EventType", m_eventType = "isBhabha");
00016 declareProperty("CalibFlag",m_calibflag=63);
00017 declareProperty("PhShapeFlag",m_phShapeFlag=0);
00018 declareProperty("TruncRate",truncate = 0.7);
00019 declareProperty("RecFileList", m_recFileList = "no file");
00020 declareProperty("RootFile",m_rootfile = "no file");
00021 declareProperty("CurveFile",m_curvefile = std::string("no rootfile"));
00022 }
00023
00024 StatusCode DedxCalib::initialize()
00025 {
00026 MsgStream log(msgSvc(), name());
00027 log << MSG::INFO << "DedxCalib initialze()" << endreq;
00028
00029 StatusCode gesc = service("MdcGeomSvc", geosvc);
00030 if (gesc == StatusCode::SUCCESS)
00031 {
00032 log << MSG::INFO <<"MdcGeomSvc is running"<<endreq;
00033 }
00034 else
00035 {
00036 log << MSG::ERROR <<"MdcGeomSvc is failed"<<endreq;
00037 return StatusCode::SUCCESS;
00038 }
00039
00040 StatusCode scex = service("DedxCorrecSvc", exsvc);
00041 if (scex == StatusCode::SUCCESS) {
00042 log << MSG::INFO <<"Hi, DedxCorrectSvc is running"<<endl;
00043 } else {
00044 log << MSG::ERROR <<"DedxCorrectSvc is failed"<<endl;
00045 return StatusCode::SUCCESS;
00046 }
00047 exsvc->set_flag( m_calibflag );
00048
00049 initializing();
00050 BookHists();
00051 getCurvePar();
00052
00053 return StatusCode::SUCCESS;
00054 }
00055
00056 StatusCode DedxCalib::execute()
00057 {
00058 MsgStream log(msgSvc(), name());
00059 log << MSG::INFO << "DedxCalib execute()" << endreq;
00060
00061 genNtuple();
00062 FillHists();
00063 AnalyseHists();
00064 WriteHists();
00065
00066 return StatusCode::SUCCESS;
00067 }
00068
00069 StatusCode DedxCalib::finalize()
00070 {
00071 MsgStream log(msgSvc(), name());
00072 log << MSG::INFO << "DedxCalib finalize()" << endreq;
00073
00074 return StatusCode::SUCCESS;
00075 }
00076
00077 void DedxCalib::initializing() {}
00078
00079 void DedxCalib::genNtuple() {}
00080
00081 void DedxCalib::BookHists() {}
00082
00083 void DedxCalib::FillHists() {}
00084
00085 void DedxCalib::AnalyseHists() {}
00086
00087 void DedxCalib::WriteHists() {}
00088
00089 void DedxCalib::ReadRecFileList()
00090 {
00091
00092 MsgStream log(msgSvc(), name());
00093 log<<MSG::INFO<<"DedxCalib::ReadRecFileList()"<<endreq;
00094
00095 ifstream filea(m_recFileList.c_str(),ifstream::in);
00096 cout<<m_recFileList.c_str()<<endl;
00097
00098 string runlist;
00099 filea>>runlist;
00100 do
00101 {
00102 m_recFileVector.push_back(runlist);
00103 cout<<runlist.c_str()<<endl;
00104 filea>>runlist;
00105 }while(filea);
00106 }
00107
00108 double DedxCalib::cal_dedx_bitrunc(float truncate, std::vector<double> phlist)
00109 {
00110 sort(phlist.begin(),phlist.end());
00111 int smpl = (int)(phlist.size()*(truncate+0.05));
00112 int min_cut = (int)( phlist.size()*0.05 + 0.5 );
00113 double qSum = 0;
00114 unsigned i = 0;
00115 for(vector<double>::iterator ql= phlist.begin();ql!=phlist.end();ql++)
00116 {
00117 i++;
00118 if(i<= smpl && i>=min_cut ) qSum += (*ql);
00119 }
00120 double trunc=-99;
00121 int usedhit = smpl-min_cut+1;
00122 if(usedhit>0) trunc=qSum/usedhit;
00123
00124 return trunc;
00125 }
00126
00127 double DedxCalib::cal_dedx(float truncate, std::vector<double> phlist)
00128 {
00129 sort(phlist.begin(),phlist.end());
00130 int smpl = (int)(phlist.size()*(truncate+0.05));
00131 int min_cut = 0;
00132 double qSum = 0;
00133 unsigned i = 0;
00134 for(vector<double>::iterator ql= phlist.begin();ql!=phlist.end();ql++)
00135 {
00136 i++;
00137 if(i<= smpl && i>=min_cut ) qSum += (*ql);
00138 }
00139 double trunc=-99;
00140 int usedhit = smpl-min_cut+1;
00141 if(usedhit>0) trunc=qSum/usedhit;
00142
00143 return trunc;
00144 }
00145
00146 void DedxCalib::getCurvePar()
00147
00148 {
00149 if(m_curvefile=="no rootfile")
00150 {
00151 cout<<"no curve file! can not calculate chi!!! "<<endl;
00152 return;
00153 }
00154
00155 double Curve[100];
00156 double Sigma[100];
00157 int CurveSize,SigmaSize;
00158 TFile* f = new TFile(m_curvefile.c_str());
00159 TTree *curve = (TTree*) f->Get("curve");
00160 TTree *sigma = (TTree*) f->Get("sigma");
00161 curve->SetBranchAddress("CurveSize",&CurveSize);
00162 curve->SetBranchAddress("curve",Curve);
00163 sigma->SetBranchAddress("SigmaSize",&SigmaSize);
00164 sigma->SetBranchAddress("sigma",Sigma);
00165 curve->GetEntry(0);
00166 sigma->GetEntry(0);
00167
00168 for(int i=0; i<CurveSize; i++)
00169 {
00170 cout<<Curve[i]<<endl;
00171 if(i==0) vFlag[0] = (int) Curve[i];
00172 else Curve_Para.push_back(Curve[i]);
00173 }
00174 for(int i=0; i<SigmaSize; i++)
00175 {
00176 cout<<Sigma[i]<<endl;
00177 if(i==0) vFlag[1] = (int) Sigma[i];
00178 else Sigma_Para.push_back(Sigma[i]);
00179 }
00180
00181
00182 }
00183
00184 void DedxCalib::set_dEdx( int landau, float dEdx, int trkalg, int runflag, int dedxhit_use, float ptrk, float theta, float pl_rp, int vflag[3] , double t0)
00185 {
00186
00187
00188
00189 if(runflag>0) vFlag[2]=0;
00190 else vFlag[2]=1;
00191
00192
00193 if(runflag ==1 || runflag ==2 )
00194 dedx_pid_exp_old( landau, runflag, (float)dEdx, (int)dedxhit_use,
00195 (float)ptrk, (float)theta, (float)t0,(float)pl_rp,
00196 m_dedx_exp, m_ex_sigma, m_pid_prob, m_chi_ex);
00197
00198 else
00199 dedx_pid_exp( vflag, (float)dEdx, trkalg,(int)dedxhit_use,
00200 (float)ptrk, (float)theta, (float)t0,(float)pl_rp,
00201 m_dedx_exp, m_ex_sigma, m_pid_prob, m_chi_ex, Curve_Para, Sigma_Para);
00202 }