00001 #include "ValidRecMdcKalTrack/ValidRecMdcKalTrackAlg.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "GaudiKernel/IHistogramSvc.h"
00005 #include "MdcRecEvent/RecMdcKalTrack.h"
00006 #include "Identifier/MdcID.h"
00007 using namespace std;
00008
00009
00010 ValidRecMdcKalTrackAlg::ValidRecMdcKalTrackAlg(const std::string& name, ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator),_ntuple(1),m_nt1(0) {
00011
00012 declareProperty("ntuple",_ntuple=1);
00013
00014 }
00015
00016
00017 void ValidRecMdcKalTrackAlg::histDef(void){
00018 if(_ntuple) {
00019 NTuplePtr nt1(ntupleSvc(),"FILE_KalValid/n101");
00020 StatusCode status;
00021 if ( nt1 ) m_nt1 = nt1;
00022 else {
00023 m_nt1= ntupleSvc()->book("FILE_KalValid/n101",CLID_ColumnWiseTuple,"KalFitValid");
00024 if ( m_nt1 ) {
00025 status = m_nt1->addItem("stat",5,2,m_stat);
00026 status = m_nt1->addItem("ndf",5,2,m_ndf);
00027 status = m_nt1->addItem("chisq",5,2,m_chisq);
00028 status = m_nt1->addItem("length",5,m_length);
00029 status = m_nt1->addItem("tof",5,m_tof);
00030 status = m_nt1->addItem("nhits",5,m_nhits);
00031 status = m_nt1->addItem("zhelix",5,m_zhelix);
00032 status = m_nt1->addItem("zhelixe",5,m_zhelixe);
00033 status = m_nt1->addItem("zhelixmu",5,m_zhelixmu);
00034 status = m_nt1->addItem("zhelixk",5,m_zhelixk);
00035 status = m_nt1->addItem("zhelixp",5,m_zhelixp);
00036 status = m_nt1->addItem("zptot",m_zptot);
00037 status = m_nt1->addItem("zptote",m_zptote);
00038 status = m_nt1->addItem("zptotmu",m_zptotmu);
00039 status = m_nt1->addItem("zptotk",m_zptotk);
00040 status = m_nt1->addItem("zptotp",m_zptotp);
00041
00042 status = m_nt1->addItem("zpt",m_zpt);
00043 status = m_nt1->addItem("zpte",m_zpte);
00044 status = m_nt1->addItem("zptmu",m_zptmu);
00045 status = m_nt1->addItem("zptk",m_zptk);
00046 status = m_nt1->addItem("zptp",m_zptp);
00047
00048 status = m_nt1->addItem("fptot",m_fptot);
00049 status = m_nt1->addItem("fptote",m_fptote);
00050 status = m_nt1->addItem("fptotmu",m_fptotmu);
00051 status = m_nt1->addItem("fptotk",m_fptotk);
00052 status = m_nt1->addItem("fptotp",m_fptotp);
00053 status = m_nt1->addItem("fpt",m_fpt);
00054 status = m_nt1->addItem("fpte",m_fpte);
00055 status = m_nt1->addItem("fptmu",m_fptmu);
00056
00057 status = m_nt1->addItem("fptk",m_fptk);
00058 status = m_nt1->addItem("fptp",m_fptp);
00059 status = m_nt1->addItem("lptot",m_lptot);
00060 status = m_nt1->addItem("lptote",m_lptote);
00061 status = m_nt1->addItem("lptotmu",m_lptotmu);
00062 status = m_nt1->addItem("lptotk",m_lptotk);
00063 status = m_nt1->addItem("lptotp",m_lptotp);
00064 status = m_nt1->addItem("lpt",m_lpt);
00065 status = m_nt1->addItem("lpte",m_lpte);
00066 status = m_nt1->addItem("lptmu",m_lptmu);
00067 status = m_nt1->addItem("lptk",m_lptk);
00068 status = m_nt1->addItem("lptp",m_lptp);
00069
00070 status = m_nt1->addItem("zsigp",m_zsigp);
00071 status = m_nt1->addItem("zsigpe",m_zsigpe);
00072 status = m_nt1->addItem("zsigpmu",m_zsigpmu);
00073 status = m_nt1->addItem("zsigpk",m_zsigpk);
00074 status = m_nt1->addItem("zsigpp",m_zsigpp);
00075 status = m_nt1->addItem("fhelix",5,m_fhelix);
00076 status = m_nt1->addItem("fhelixe",5,m_fhelixe);
00077 status = m_nt1->addItem("fhelixmu",5,m_fhelixmu);
00078 status = m_nt1->addItem("fhelixk",5,m_fhelixk);
00079 status = m_nt1->addItem("fhelixp",5,m_fhelixp);
00080 status = m_nt1->addItem("lhelix",5,m_lhelix);
00081 status = m_nt1->addItem("lhelixe",5,m_lhelixe);
00082 status = m_nt1->addItem("lhelixmu",5,m_lhelixmu);
00083 status = m_nt1->addItem("lhelixk",5,m_lhelixk);
00084 status = m_nt1->addItem("lhelixp",5,m_lhelixp);
00085
00086 status = m_nt1->addItem("zerror",15,m_zerror);
00087 status = m_nt1->addItem("zerrore",15,m_zerrore);
00088 status = m_nt1->addItem("zerrormu",15,m_zerrormu);
00089 status = m_nt1->addItem("zerrork",15,m_zerrork);
00090 status = m_nt1->addItem("zerrorp",15,m_zerrorp);
00091 status = m_nt1->addItem("ferror",15,m_ferror);
00092 status = m_nt1->addItem("ferrore",15,m_ferrore);
00093 status = m_nt1->addItem("ferrormu",15,m_ferrormu);
00094 status = m_nt1->addItem("ferrork",15,m_ferrork);
00095 status = m_nt1->addItem("ferrorp",15,m_ferrorp);
00096 status = m_nt1->addItem("lerror",15,m_lerror);
00097 status = m_nt1->addItem("lerrore",15,m_lerrore);
00098 status = m_nt1->addItem("lerrormu",15,m_lerrormu);
00099 status = m_nt1->addItem("lerrork",15,m_lerrork);
00100 status = m_nt1->addItem("lerrorp",15,m_lerrorp);
00101
00102 if( status.isFailure() ) cout<<"Ntuple m_nt1 add item failed!"<<endl;
00103 }
00104 }
00105
00106
00107 NTuplePtr nt2(ntupleSvc(),"FILE_KalValid/n102");
00108 StatusCode status2;
00109 if ( nt2 ) m_nt2 = nt2;
00110 else {
00111 m_nt2= ntupleSvc()->book("FILE_KalValid/n102",CLID_ColumnWiseTuple,"KalFitValid");
00112 if ( m_nt2 ) {
00113 if( status2.isFailure() ) cout<<"Ntuple2 add item failed!"<<endl;
00114 }
00115 }
00116 }
00117 }
00118
00119
00120 StatusCode ValidRecMdcKalTrackAlg::initialize(){
00121
00122 MsgStream log( msgSvc(), name() );
00123 log << MSG::INFO << "Booking Histograms" << endreq;
00124
00125 histDef();
00126
00127 log << MSG::INFO << "Finished booking Histograms" << endreq;
00128 return StatusCode::SUCCESS;
00129 }
00130
00131
00132 StatusCode ValidRecMdcKalTrackAlg::execute() {
00133 MsgStream log(msgSvc(), name());
00134 SmartDataPtr<RecMdcKalTrackCol> mdcKalTrkCol(eventSvc(),"/Event/Recon/RecMdcKalTrackCol");
00135 if( ! mdcKalTrkCol )
00136 {
00137 log << MSG::ERROR << "Unable to retrieve mdcKalTrackCol" << endreq;
00138 return StatusCode::FAILURE;
00139 } else {
00140 for(RecMdcKalTrackCol::iterator iter_trk=mdcKalTrkCol->begin(); iter_trk!=mdcKalTrkCol->end(); iter_trk++) {
00141 for( int jj=0, iii=0; jj<5; jj++){
00142 m_nhits[jj] = (*iter_trk)->getNhits(jj);
00143 m_zhelix[jj] = (*iter_trk)->getZHelix()[jj];
00144 m_zhelixe[jj] = (*iter_trk)->getZHelixE()[jj];
00145 m_zhelixmu[jj] = (*iter_trk)->getZHelixMu()[jj];
00146 m_zhelixk[jj] = (*iter_trk)->getZHelixK()[jj];
00147 m_zhelixp[jj] = (*iter_trk)->getZHelixP()[jj];
00148 m_fhelix[jj] = (*iter_trk)->getFHelix()[jj];
00149 m_fhelixe[jj] = (*iter_trk)->getFHelixE()[jj];
00150 m_fhelixmu[jj] = (*iter_trk)->getFHelixMu()[jj];
00151 m_fhelixk[jj] = (*iter_trk)->getFHelixK()[jj];
00152 m_fhelixp[jj] = (*iter_trk)->getFHelixP()[jj];
00153 m_lhelix[jj] = (*iter_trk)->getLHelix()[jj];
00154 m_lhelixe[jj] = (*iter_trk)->getLHelixE()[jj];
00155 m_lhelixmu[jj] = (*iter_trk)->getLHelixMu()[jj];
00156 m_lhelixk[jj] = (*iter_trk)->getLHelixK()[jj];
00157 m_lhelixp[jj] = (*iter_trk)->getLHelixP()[jj];
00158
00159
00160 for(int kk=0; kk<=jj; kk++,iii++) {
00161 m_zerror[iii] = (*iter_trk)->getZError()[jj][kk];
00162 m_zerrore[iii] = (*iter_trk)->getZErrorE()[jj][kk];
00163 m_zerrormu[iii] = (*iter_trk)->getZErrorMu()[jj][kk];
00164 m_zerrork[iii] = (*iter_trk)->getZErrorK()[jj][kk];
00165 m_zerrorp[iii] = (*iter_trk)->getZErrorP()[jj][kk];
00166 m_ferror[iii] = (*iter_trk)->getFError()[jj][kk];
00167 m_ferrore[iii] = (*iter_trk)->getFErrorE()[jj][kk];
00168 m_ferrormu[iii] = (*iter_trk)->getFErrorMu()[jj][kk];
00169 m_ferrork[iii] = (*iter_trk)->getFErrorK()[jj][kk];
00170 m_ferrorp[iii] = (*iter_trk)->getFErrorP()[jj][kk];
00171 m_lerror[iii] = (*iter_trk)->getLError()[jj][kk];
00172 m_lerrore[iii] = (*iter_trk)->getLErrorE()[jj][kk];
00173 m_lerrormu[iii] = (*iter_trk)->getLErrorMu()[jj][kk];
00174 m_lerrork[iii] = (*iter_trk)->getLErrorK()[jj][kk];
00175 m_lerrorp[iii] = (*iter_trk)->getLErrorP()[jj][kk];
00176 }
00177 }
00178
00179
00180 m_ndf[0][0] = (*iter_trk)->getNdf(0,0);
00181 m_ndf[1][0] = (*iter_trk)->getNdf(0,1);
00182 m_ndf[2][0] = (*iter_trk)->getNdf(0,2);
00183 m_ndf[3][0] = (*iter_trk)->getNdf(0,3);
00184 m_ndf[4][0] = (*iter_trk)->getNdf(0,4);
00185 m_ndf[0][1] = (*iter_trk)->getNdf(1,0);
00186 m_ndf[1][1] = (*iter_trk)->getNdf(1,1);
00187 m_ndf[2][1] = (*iter_trk)->getNdf(1,2);
00188 m_ndf[3][1] = (*iter_trk)->getNdf(1,3);
00189 m_ndf[4][1] = (*iter_trk)->getNdf(1,4);
00190
00191 m_stat[0][0] = (*iter_trk)->getStat(0,0);
00192 m_stat[1][0] = (*iter_trk)->getStat(0,1);
00193 m_stat[2][0] = (*iter_trk)->getStat(0,2);
00194 m_stat[3][0] = (*iter_trk)->getStat(0,3);
00195 m_stat[4][0] = (*iter_trk)->getStat(0,4);
00196 m_stat[0][1] = (*iter_trk)->getStat(1,0);
00197 m_stat[1][1] = (*iter_trk)->getStat(1,1);
00198 m_stat[2][1] = (*iter_trk)->getStat(1,2);
00199 m_stat[3][1] = (*iter_trk)->getStat(1,3);
00200 m_stat[4][1] = (*iter_trk)->getStat(1,4);
00201
00202
00203 m_fptot = sqrt(1+pow(m_fhelix[4],2))/m_fhelix[2];
00204 m_fptote = sqrt(1+pow(m_fhelixe[4],2))/m_fhelixe[2];
00205 m_fptotmu = sqrt(1+pow(m_fhelixmu[4],2))/m_fhelixmu[2];
00206 m_fptotk = sqrt(1+pow(m_fhelixk[4],2))/m_fhelixk[2];
00207 m_fptotp = sqrt(1+pow(m_fhelixp[4],2))/m_fhelixp[2];
00208
00209 m_zpt = 1/m_zhelix[2];
00210 m_zpte = 1/m_zhelixe[2];
00211 m_zptmu = 1/m_zhelixmu[2];
00212 m_zptk = 1/m_zhelixk[2];
00213 m_zptp = 1/m_zhelixp[2];
00214
00215 m_fpt = 1/m_fhelix[2];
00216 m_fpte = 1/m_fhelixe[2];
00217 m_fptmu = 1/m_fhelixmu[2];
00218 m_fptk = 1/m_fhelixk[2];
00219 m_fptp = 1/m_fhelixp[2];
00220
00221 m_lpt = 1/m_lhelix[2];
00222 m_lpte = 1/m_lhelixe[2];
00223 m_lptmu = 1/m_lhelixmu[2];
00224 m_lptk = 1/m_lhelixk[2];
00225 m_lptp = 1/m_lhelixp[2];
00226
00227 m_lptot = sqrt(1+pow(m_lhelix[4],2))/m_lhelix[2];
00228 m_lptote = sqrt(1+pow(m_lhelixe[4],2))/m_lhelixe[2];
00229 m_lptotmu = sqrt(1+pow(m_lhelixmu[4],2))/m_lhelixmu[2];
00230 m_lptotk = sqrt(1+pow(m_lhelixk[4],2))/m_lhelixk[2];
00231 m_lptotp = sqrt(1+pow(m_lhelixp[4],2))/m_lhelixp[2];
00232
00233 m_zptot = sqrt(1+pow(m_zhelix[4],2))/m_zhelix[2];
00234 m_zptote = sqrt(1+pow(m_zhelixe[4],2))/m_zhelixe[2];
00235 m_zptotmu = sqrt(1+pow(m_zhelixmu[4],2))/m_zhelixmu[2];
00236 m_zptotk = sqrt(1+pow(m_zhelixk[4],2))/m_zhelixk[2];
00237 m_zptotp = sqrt(1+pow(m_zhelixp[4],2))/m_zhelixp[2];
00238
00239 StatusCode sc1 = m_nt1->write();
00240 if( sc1.isFailure() ) cout<<"Ntuple1 filling failed!"<<endl;
00241 }
00242 }
00243 return StatusCode::SUCCESS;
00244 }
00245
00246
00247 StatusCode ValidRecMdcKalTrackAlg::finalize() {
00248
00249 MsgStream log(msgSvc(), name());
00250 log << MSG::INFO << "Finalizing..." << endreq;
00251 return StatusCode::SUCCESS;
00252 }
00253
00254
00255
