00001 #include "MdcNavigation/MdcNavigation.h"
00002 #include "GaudiKernel/MsgStream.h"
00003 #include "GaudiKernel/SmartDataPtr.h"
00004 #include "EventNavigator/EventNavigator.h"
00005 #include <cmath>
00006 #include "HepPDT/ParticleDataTable.hh"
00007 #include "HepPDT/ParticleData.hh"
00008 #include "GaudiKernel/IPartPropSvc.h"
00009 #include "CLHEP/Geometry/Point3D.h"
00010 #include "EventModel/EventModel.h"
00011 #include "EventModel/Event.h"
00012 #include "EventModel/EventHeader.h"
00013 #include "RawEvent/RawDataUtil.h"
00014 #include "MdcGeom/Constants.h"
00015
00016 #include "MdcGeom/BesAngle.h"
00017 #include "TrkBase/HelixTraj.h"
00018 #include "CLHEP/Matrix/SymMatrix.h"
00019 #include "TrkBase/TrkPocaBase.h"
00020 #include "TrkBase/TrkPoca.h"
00021 #include "MdcGeom/MdcLayer.h"
00022 #include "MdcGeom/MdcDetector.h"
00023 #include "MdcData/MdcHit.h"
00024
00025 #include "Identifier/MdcID.h"
00026 #include "MdcRawEvent/MdcDigi.h"
00027 #include "EvTimeEvent/RecEsTime.h"
00028
00029 #include "GaudiKernel/INTupleSvc.h"
00030 #include "GaudiKernel/NTuple.h"
00031 using namespace Event;
00032 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00033
00034 typedef HepGeom::Point3D<double> HepPoint3D;
00035 #endif
00036
00037
00038
00039
00040 using namespace std;
00041 using namespace Event;
00042 const double epsilon = 0.000000001;
00043
00044 MdcNavigation::MdcNavigation(const std::string& name, ISvcLocator* pSvcLocator) :
00045 Algorithm(name, pSvcLocator) {
00046 declareProperty("hist", m_hist = 0);
00047 declareProperty("nMcHit", m_nMcHit= 5);
00048 declareProperty("mc", m_mc = 1);
00049
00050 declareProperty("maxMdcDigi", m_maxMdcDigi= 0);
00051 declareProperty("keepBadTdc", m_keepBadTdc= 0);
00052 declareProperty("dropHot", m_dropHot= 0);
00053 declareProperty("keepUnmatch", m_keepUnmatch= 0);
00054
00055 declareProperty("poca", m_poca = false);
00056 declareProperty("doSag", m_doSag = false);
00057
00058 declareProperty("d0Cut", m_d0Cut= 1.);
00059 declareProperty("z0Cut", m_z0Cut= 10.);
00060 declareProperty("debug", m_debug = 0);
00061
00062 }
00063
00064
00065
00066
00067 StatusCode MdcNavigation::initialize(){
00068 MsgStream log(msgSvc(), name());
00069 StatusCode sc = StatusCode::SUCCESS;
00070
00071 t_nTk = 0;
00072
00073 sc = service ("MagneticFieldSvc",m_pIMF);
00074 if(sc!=StatusCode::SUCCESS) {
00075 log << MSG::ERROR << "Unable to open Magnetic field service"<<endreq;
00076 return StatusCode::FAILURE;
00077 }
00078
00079
00080
00081 IPartPropSvc* p_PartPropSvc;
00082 static const bool CREATEIFNOTTHERE(true);
00083 sc = service("PartPropSvc", p_PartPropSvc, CREATEIFNOTTHERE);
00084 if (!sc.isSuccess() || 0 == p_PartPropSvc) {
00085 log << MSG::ERROR << " Could not initialize PartPropSvc" << endreq;
00086 return sc;
00087 }
00088
00089 m_particleTable = p_PartPropSvc->PDT();
00090
00091 IRawDataProviderSvc* irawDataProviderSvc;
00092 sc = service ("RawDataProviderSvc", irawDataProviderSvc);
00093 m_rawDataProviderSvc = dynamic_cast<RawDataProviderSvc*> (irawDataProviderSvc);
00094 if ( sc.isFailure() ){
00095 log << MSG::FATAL << "Could not load RawDataProviderSvc!" << endreq;
00096 return StatusCode::FAILURE;
00097 }
00098
00099
00100 if (m_hist){
00101 sc = bookNTuple();
00102 if (!sc.isSuccess()) {
00103 log << MSG::WARNING << " Could not book NTuple" << endreq;
00104 m_hist = 0;
00105 }
00106 }
00107
00108 keepedParticles = new int(211);
00109
00110
00111 return StatusCode::SUCCESS;
00112 }
00113
00114 StatusCode MdcNavigation::beginRun(){
00115 MsgStream log(msgSvc(), name());
00116 log << MSG::INFO << "in beginRun()" << endreq;
00117
00118 m_gm = MdcDetector::instance(m_doSag);
00119 if(NULL == m_gm) return StatusCode::FAILURE;
00120
00121 return StatusCode::SUCCESS;
00122 }
00123
00124
00125
00126 StatusCode MdcNavigation::execute() {
00127 setFilterPassed(false);
00128 MsgStream log(msgSvc(), name());
00129 StatusCode sc = StatusCode::SUCCESS;
00130
00131
00132 SmartDataPtr<EventNavigator> navigator (eventSvc(),"/Event/Navigator");
00133 if( ! navigator ) {
00134 log << MSG::WARNING<< " Unable to retrieve EventNavigator" << endreq;
00135 m_rawData = true;
00136 }
00137 SmartDataPtr<RecMdcTrackCol> recMdcTrackCol(eventSvc(), "/Event/Recon/RecMdcTrackCol");
00138 SmartDataPtr<RecMdcHitCol> recMdcHitCol(eventSvc(), "/Event/Recon/RecMdcHitCol");
00139
00140
00141 if(m_hist){
00142 sc = fillInit();
00143 if ( sc!=StatusCode::SUCCESS ) {
00144 return StatusCode::FAILURE;
00145 }
00146 }
00147 if (m_mc){
00148
00149 SmartDataPtr<Event::McParticleCol> mcParticles(eventSvc(),"/Event/MC/McParticleCol");
00150 SmartDataPtr<Event::MdcMcHitCol> mcHit(eventSvc(),"/Event/MC/MdcMcHitCol");
00151 if( ! mcParticles ) {
00152 log << MSG::WARNING<< " Unable to retrieve McParticleCol" << endreq;
00153 }else{
00154
00155 t_mcTkNum = 0;
00156 McParticleCol::iterator it= mcParticles->begin();
00157 log <<MSG::INFO << "mcParticles size = "<<mcParticles->size() << endreq;
00158 for(; it!= mcParticles->end(); it++ ) {
00159
00160 t_mcTkNum++;
00161 }
00162 }
00163 }
00164 t_mcTkNum = 0;
00165 t_recTkNum = 0;
00166
00167
00168 if (!recMdcTrackCol){
00169 log << MSG::WARNING<< " Unable to retrieve recMdcTrackCol" << endreq;
00170 return StatusCode::SUCCESS;
00171 }
00172 t_recTkNum = recMdcTrackCol->size();
00173
00174
00175 RecMdcTrackCol::iterator it = recMdcTrackCol->begin();
00176 for(;it != recMdcTrackCol->end(); it++ ) {
00177 if ( m_mc && navigator ) {
00178 McParticleVector particles = navigator->getMcParticles(*it);
00179 t_mcTkNum = particles.size();
00180 RecMdcTrackVector tracks = navigator->getMdcTracks(particles[0]);
00181
00182 if ( sc!=StatusCode::SUCCESS ) return StatusCode::FAILURE;
00183 }
00184 sc = fillRecTrack(*it, t_mcTkNum, t_recTkNum);
00185 t_nTk++;
00186 if ( sc!=StatusCode::SUCCESS ) return StatusCode::FAILURE;
00187 }
00188
00189
00190 fillRecHits(*recMdcHitCol);
00191
00192
00193 if(m_hist){fillEvent();}
00194
00195 return StatusCode::SUCCESS;
00196 }
00197
00198
00199
00200 StatusCode MdcNavigation::finalize() {
00201 MsgStream log(msgSvc(), name());
00202 log << MSG::INFO << "in finalize()" << endreq;
00203 std::cout<< "nTk == "<<t_nTk << std::endl;
00204 delete keepedParticles;
00205 return StatusCode::SUCCESS;
00206 }
00207
00208
00209 Hep3Vector MdcNavigation::momentum(const RecMdcTrack* trk) {
00210
00211 double fi0 = trk->helix(1);
00212 double cpa = trk->helix(2);
00213 double tanl = trk->helix(4);
00214
00215 double pxy = 0.;
00216 if(cpa != 0) pxy = 1/fabs(cpa);
00217
00218 double px = pxy * (-sin(fi0));
00219 double py = pxy * cos(fi0);
00220 double pz = pxy * tanl;
00221
00222 Hep3Vector p;
00223 p.set(px, py, pz);
00224 return p;
00225 }
00226
00227 StatusCode MdcNavigation::fillRecTrack(const RecMdcTrack* tk, int mcTkNum, int recTkNum){
00228
00229
00230
00231 m_na_nEvtNoise = nNoise;
00232
00233 m_na_recTkNum = recTkNum;
00234 CLHEP::Hep3Vector rec_mom = momentum(tk);
00235
00236 m_na_p = rec_mom.mag();
00237 m_na_pt = rec_mom.perp();
00238 m_na_pz = rec_mom.z();
00239
00240
00241
00242 m_na_d0 = tk->helix(0);
00243 m_na_phi0 = tk->helix(1);
00244 m_na_cpa = tk->helix(2);
00245 m_na_z0 = tk->helix(3);
00246 m_na_tanl = tk->helix(4);
00247
00248 if( m_na_d0 > m_d0Cut ) {
00249 std::cout<<__FILE__<<" "<<__LINE__<<" evtNo: "<<t_eventNo<<" d0:"<<std::setw(5)<<m_na_d0<<">"<<m_d0Cut<< std::endl;
00250 setFilterPassed(true);
00251 }
00252 if( m_na_z0 > m_z0Cut ) {
00253 std::cout<<__FILE__<<" "<<__LINE__<<" evtNo: "<<t_eventNo<<" z0:"<<std::setw(5)<<m_na_z0<<">"<<m_z0Cut<< std::endl;
00254 setFilterPassed(true);
00255 }
00256
00257 m_na_d0E = tk->err(0);
00258 m_na_phi0E = tk->err(2);
00259 m_na_cpaE = tk->err(5);
00260 m_na_z0E = tk->err(9);
00261 m_na_tanlE = tk->err(14);
00262 m_na_q = tk->charge();
00263
00264 m_na_chi2 = tk->chi2();
00265
00266 m_na_nDof = tk->ndof();
00267 if ( m_na_nDof > 0 ) {
00268 m_na_chi2Dof = m_na_chi2/(float)m_na_nDof;
00269 } else {
00270 m_na_chi2Dof = 200.;
00271 }
00272 m_na_chi2Prob = probab(m_na_nDof, m_na_chi2);
00273 m_na_nSt = tk->nster();
00274 m_na_fiTerm = tk->getFiTerm();
00275
00276 if (tk->stat()==0){
00277 t_trkRecoTk++;
00278 }else if(tk->stat()==1){
00279 t_curlTk++;
00280 }else if(tk->stat()==2){
00281 t_patRecTk++;
00282 }else if(tk->stat()==3){
00283 t_xRecTk++;
00284 }
00285 m_na_tkStat = tk->stat();
00286
00288
00289
00290
00291
00292
00293 int noiseHit=0;
00294 int matchHit=0;
00295 int nAct = 0;
00296 HitRefVec hl = tk->getVecHits();
00297 HitRefVec::iterator hitIt = hl.begin();
00298 for (;hitIt!=hl.end();++hitIt){
00299 RecMdcHit* h = *hitIt;
00300 if ( !h ) continue;
00301 if (h->getStat()!=0) nAct++;
00302
00303
00304
00305
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315
00316
00317
00318
00319
00320
00321
00322
00323 int tlayer = MdcID::layer(h->getMdcId());
00324 int twire = MdcID::wire(h->getMdcId());
00325
00326
00327
00328
00329
00330
00331
00332
00333 if(havedigi[tlayer][twire]<0){
00334 ++noiseHit;
00335 }
00336
00337
00338
00339
00340
00341 }
00342 m_na_nAct = nAct;
00343 m_na_nNoise = noiseHit;
00344 m_na_nMatch = matchHit;
00345 g_tupleMc->write();
00346
00347
00348
00349
00350 if (m_poca){
00351 uint32_t getDigiFlag = 0;
00352 getDigiFlag += m_maxMdcDigi;
00353 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00354 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00355 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00356 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00357 MdcDigiCol::iterator iter = mdcDigiVec.begin();
00358 for (;iter != mdcDigiVec.end(); iter++) {
00359 poca((*iter),tk->helix(),tk->err());
00360 }
00361 }
00362
00363 return StatusCode::SUCCESS;
00364 }
00365
00366 StatusCode MdcNavigation::bookNTuple(){
00367 MsgStream log(msgSvc(), name());
00368 StatusCode sc = StatusCode::SUCCESS;
00369 g_layerEff = histoSvc()->book( "layerEff", "layerEff",43,-0.5,42.5 );
00370
00371 NTuplePtr nt1(ntupleSvc(), "MdcNavigation/rec");
00372 if ( nt1 ) { g_tupleMc = nt1;}
00373 else {
00374 g_tupleMc = ntupleSvc()->book ("MdcNavigation/rec", CLID_ColumnWiseTuple, "rec and delta with mc truth");
00375 if ( g_tupleMc ) {
00376 sc = g_tupleMc->addItem ("tkStat", m_na_tkStat);
00377 sc = g_tupleMc->addItem ("p", m_na_p);
00378 sc = g_tupleMc->addItem ("pt", m_na_pt);
00379 sc = g_tupleMc->addItem ("pz", m_na_pz);
00380 sc = g_tupleMc->addItem ("d0", m_na_d0);
00381 sc = g_tupleMc->addItem ("phi0", m_na_phi0);
00382 sc = g_tupleMc->addItem ("cpa", m_na_cpa);
00383 sc = g_tupleMc->addItem ("z0", m_na_z0);
00384 sc = g_tupleMc->addItem ("tanl", m_na_tanl);
00385 sc = g_tupleMc->addItem ("d0E", m_na_d0E);
00386 sc = g_tupleMc->addItem ("phi0E", m_na_phi0E);
00387 sc = g_tupleMc->addItem ("cpaE", m_na_cpaE);
00388 sc = g_tupleMc->addItem ("z0E", m_na_z0E);
00389 sc = g_tupleMc->addItem ("tanlE", m_na_tanlE);
00390 sc = g_tupleMc->addItem ("q", m_na_q);
00391
00392 sc = g_tupleMc->addItem ("dP", m_na_dP);
00393 sc = g_tupleMc->addItem ("dPt", m_na_dPt);
00394 sc = g_tupleMc->addItem ("dPz", m_na_dPz);
00395 sc = g_tupleMc->addItem ("dD0", m_na_dD0);
00396 sc = g_tupleMc->addItem ("dPhi0", m_na_dPhi0);
00397 sc = g_tupleMc->addItem ("dCpa", m_na_dCpa);
00398 sc = g_tupleMc->addItem ("dZ0", m_na_dZ0);
00399 sc = g_tupleMc->addItem ("dTanl", m_na_dTanl);
00400
00401 sc = g_tupleMc->addItem ("d0Res", m_na_d0Res);
00402 sc = g_tupleMc->addItem ("phiRes", m_na_phi0Res);
00403 sc = g_tupleMc->addItem ("z0Res", m_na_z0Res);
00404 sc = g_tupleMc->addItem ("tanlRes", m_na_tanlRes);
00405 sc = g_tupleMc->addItem ("cpaRes", m_na_cpaRes);
00406
00407 sc = g_tupleMc->addItem ("recTkNum",m_na_recTkNum);
00408 sc = g_tupleMc->addItem ("mcTkId", m_na_mcTkId);
00409 sc = g_tupleMc->addItem ("mcTkNum", m_na_mcTkNum);
00410 sc = g_tupleMc->addItem ("evtNo", m_na_evtNo);
00411 sc = g_tupleMc->addItem ("nSt", m_na_nSt);
00412 sc = g_tupleMc->addItem ("nDof", m_na_nDof);
00413 sc = g_tupleMc->addItem ("chi2", m_na_chi2);
00414 sc = g_tupleMc->addItem ("chi2Dof", m_na_chi2Dof);
00415 sc = g_tupleMc->addItem ("chi2Porb",m_na_chi2Prob);
00416 sc = g_tupleMc->addItem ("fiTerm", m_na_fiTerm);
00417 sc = g_tupleMc->addItem ("nMatch", m_na_nMatch);
00418 sc = g_tupleMc->addItem ("nAct", m_na_nAct);
00419 sc = g_tupleMc->addItem ("nTkNoise",m_na_nNoise);
00420 sc = g_tupleMc->addItem ("nEvtNoise",m_na_nEvtNoise);
00421 sc = g_tupleMc->addItem ("nHit", m_na_nHit, 0, 10000);
00422 sc = g_tupleMc->addItem ("nDigi", m_na_nDigi, 0, 10000);
00423 sc = g_tupleMc->addIndexedItem ("resid", m_na_nHit, m_na_resid);
00424 sc = g_tupleMc->addIndexedItem ("driftD", m_na_nHit, m_na_driftD);
00425 sc = g_tupleMc->addIndexedItem ("driftT", m_na_nHit, m_na_driftT);
00426 sc = g_tupleMc->addIndexedItem ("doca", m_na_nHit, m_na_doca);
00427 sc = g_tupleMc->addIndexedItem ("entra", m_na_nHit, m_na_entra);
00428 sc = g_tupleMc->addIndexedItem ("zhit", m_na_nHit, m_na_zhit);
00429 sc = g_tupleMc->addIndexedItem ("chi2add", m_na_nHit, m_na_chi2add);
00430 sc = g_tupleMc->addIndexedItem ("flaglr", m_na_nHit, m_na_flaglr);
00431 sc = g_tupleMc->addIndexedItem ("hitStat", m_na_nHit, m_na_hitStat);
00432 sc = g_tupleMc->addIndexedItem ("tdc", m_na_nHit, m_na_Tdc);
00433 sc = g_tupleMc->addIndexedItem ("adc", m_na_nHit, m_na_Adc);
00434 sc = g_tupleMc->addIndexedItem ("act", m_na_nHit, m_na_act);
00435 sc = g_tupleMc->addIndexedItem ("layer", m_na_nHit, m_na_layer);
00436 sc = g_tupleMc->addIndexedItem ("wire", m_na_nHit, m_na_wire);
00437 sc = g_tupleMc->addIndexedItem ("gwire", m_na_nHit, m_na_gwire);
00438 sc = g_tupleMc->addIndexedItem ("hitTkId", m_na_nHit, m_na_hitTkId);
00439 sc = g_tupleMc->addIndexedItem ("digiTkId", m_na_nDigi, m_na_digiTkId);
00440 sc = g_tupleMc->addIndexedItem ("mclayer", m_na_nDigi, m_na_digiLayer);
00441 } else {
00442 log << MSG::ERROR << " Cannot book N-tuple:" << long(g_tupleMc) << endmsg;
00443 return StatusCode::FAILURE;
00444 }
00445 }
00446 NTuplePtr nt3(ntupleSvc(), "MdcNavigation/evt");
00447 if ( nt3 ) { g_tupleEvt = nt3;}
00448 else {
00449 g_tupleEvt = ntupleSvc()->book ("MdcNavigation/evt", CLID_ColumnWiseTuple, "event");
00450 if ( g_tupleEvt ) {
00451 sc = g_tupleEvt->addItem ("nTdsTk", m_na_t3recTk);
00452 sc = g_tupleEvt->addItem ("trkReco", m_na_t3TrkReco);
00453 sc = g_tupleEvt->addItem ("curlFinder", m_na_t3Curl);
00454 sc = g_tupleEvt->addItem ("patRec", m_na_t3PatRec);
00455 sc = g_tupleEvt->addItem ("xRec", m_na_t3XRec);
00456 sc = g_tupleEvt->addItem ("mcTkNum", m_na_t3mcTk);
00457 sc = g_tupleEvt->addItem ("evtNo", m_na_t3evtNo);
00458 sc = g_tupleEvt->addItem ("t0", m_na_t3t0);
00459 sc = g_tupleEvt->addItem ("t0Truth", m_na_t3t0Truth);
00460 sc = g_tupleEvt->addItem ("t0Stat", m_na_t3t0Stat);
00461 sc = g_tupleEvt->addItem ("runNo", m_na_t3runNo);
00462 sc = g_tupleEvt->addItem ("nDigi", m_na_t3nDigi, 0, 10000);
00463 sc = g_tupleEvt->addIndexedItem ("layer", m_na_t3nDigi, m_na_t3layer);
00464 sc = g_tupleEvt->addIndexedItem ("wire", m_na_t3nDigi, m_na_t3wire);
00465 sc = g_tupleEvt->addIndexedItem ("gwire", m_na_t3nDigi, m_na_t3gwire);
00466 sc = g_tupleEvt->addIndexedItem ("rt", m_na_t3nDigi, m_na_t3rt);
00467 sc = g_tupleEvt->addIndexedItem ("rtNot0",m_na_t3nDigi, m_na_t3rtNot0);
00468 sc = g_tupleEvt->addIndexedItem ("rc", m_na_t3nDigi, m_na_t3rc);
00469 sc = g_tupleEvt->addIndexedItem ("phi", m_na_t3nDigi, m_na_t3phi);
00470 sc = g_tupleEvt->addIndexedItem ("ovfl", m_na_t3nDigi, m_na_t3ovfl);
00471 sc = g_tupleEvt->addIndexedItem ("tNum", m_na_t3nDigi, m_na_t3tNum);
00472 }
00473 }
00474 }
00475
00476 StatusCode MdcNavigation::fillInit(){
00477 MsgStream log(msgSvc(), name());
00478 StatusCode sc = StatusCode::SUCCESS;
00479
00480
00481
00482 t_trkRecoTk = 0;
00483 t_curlTk = 0;
00484 t_patRecTk = 0;
00485 t_xRecTk = 0;
00486
00487 SmartDataPtr<Event::EventHeader> evtHead(eventSvc(),"/Event/EventHeader");
00488 if (evtHead) {
00489 t_runNo = evtHead->runNumber();
00490 t_eventNo = evtHead->eventNumber();
00491 }else{
00492 log << MSG::WARNING<< "Could not retrieve event header" << endreq;
00493 return StatusCode::FAILURE;
00494 }
00495 if(m_debug) std::cout<< "evtNo:"<<t_eventNo << std::endl;
00496
00497
00498 t_t0 = -1;
00499 t_t0Stat = -1;
00500 SmartDataPtr<RecEsTimeCol> aevtimeCol(eventSvc(),"/Event/Recon/RecEsTimeCol");
00501
00502 if (!aevtimeCol || aevtimeCol->size()==0) {
00503 log << MSG::WARNING << "Could not find RecEsTimeCol" << endreq;
00504 }else{
00505 RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
00506 t_t0 = (*iter_evt)->getTest();
00507 t_t0Stat = (*iter_evt)->getStat();
00508 }
00509
00510
00511 uint32_t getDigiFlag = 0;
00512 getDigiFlag += m_maxMdcDigi;
00513 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00514 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00515 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00516 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00517 if ((mdcDigiVec.size()==0)) {
00518 log << MSG::WARNING << t_eventNo <<"No digi or could not find event in MdcDigiVec" << endreq;
00519 }
00520
00521 for (int ii=0;ii<43;ii++){
00522 for (int jj=0;jj<288;jj++){
00523 havedigi[ii][jj]= -99;
00524 }
00525 }
00526
00527 for(int imc=0;imc<100;imc++){
00528 nDigiTk[imc]=0;
00529 digiLayer[imc]=0;
00530 }
00531
00532 for(int ii=0;ii<43;ii++) for(int jj=0;jj<288;jj++) multiTdcCount[ii][jj]=0;
00533 MdcDigiCol::iterator iter = mdcDigiVec.begin();
00534 for (;iter != mdcDigiVec.end(); iter++ ) {
00535 double t = RawDataUtil::MdcTime((*iter)->getTimeChannel());
00536 double c = (*iter)->getChargeChannel();
00537 int l = MdcID::layer((*iter)->identify());
00538 int w = MdcID::wire((*iter)->identify());
00539 multiTdcCount[l][w]++;
00540 }
00541
00542 int t_i=0;
00543 iter = mdcDigiVec.begin();
00544 for (;iter != mdcDigiVec.end(); iter++,++t_i ) {
00545 double t = RawDataUtil::MdcTime((*iter)->getTimeChannel());
00546 double c = (*iter)->getChargeChannel();
00547 int l = MdcID::layer((*iter)->identify());
00548 int w = MdcID::wire((*iter)->identify());
00549 if(!m_rawData){
00550 int tkid = (*iter)->getTrackIndex();
00551 havedigi[l][w]= tkid;
00552 if (g_tupleMc){
00553
00554 }
00555 if(tkid>-1){
00556 ++nDigiTk[tkid];
00557 }else{
00558 nNoise++;
00559 }
00560 }
00561
00562 }
00563 return sc;
00564 }
00565
00566 StatusCode MdcNavigation::skipMcParticle(const McParticle* it, int nKindKeeped, int* pid){
00567
00568 MsgStream log(msgSvc(), name());
00569
00570 int pdg_code = (*it).particleProperty();
00571 if (fabs(pdg_code)>=8) {
00572 const HepPDT::ParticleData* particles = m_particleTable->particle(abs(pdg_code));
00573 if( ! particles ){
00574 log<<MSG::WARNING<<"Exotic particle found with PDG code "<<pdg_code<<endreq;
00575 }else{
00576
00577 if( particles->charge() == 0 ){
00578 log<<MSG::INFO<<"Skip charge == 0 mc particle "<<pdg_code<<endreq;
00579 return StatusCode::FAILURE;
00580 }
00581 }
00582 }
00583
00584 int mcTkId = (*it).trackIndex();
00585 int nMcHit=0;
00586 SmartDataPtr<Event::MdcMcHitCol> mcMdcMcHitCol(eventSvc(),"/Event/MC/MdcMcHitCol");
00587 if (!mcMdcMcHitCol) {
00588 log << MSG::INFO << "Could not find MdcMcHit" << endreq;
00589 }else{
00590 Event::MdcMcHitCol::iterator iter_mchit = mcMdcMcHitCol->begin();
00591 for (;iter_mchit != mcMdcMcHitCol->end(); iter_mchit++ ) {
00592 int iMcTk = (*iter_mchit)->getTrackIndex();
00593 if (mcTkId == iMcTk) nMcHit++;
00594 }
00595 }
00596 if(nMcHit<m_nMcHit) return StatusCode::FAILURE;
00597
00598 bool keeped = false;
00599 for (int i=0; i<nKindKeeped; i++){
00600 if (abs(pdg_code) == pid[i]) keeped = true;
00601 }
00602
00603 if (!keeped) return StatusCode::FAILURE;
00604
00605 return StatusCode::SUCCESS;
00606 }
00607
00608 StatusCode MdcNavigation::fillEvent(){
00609 if (!g_tupleEvt) return StatusCode::FAILURE;
00610 uint32_t getDigiFlag = 0;
00611 getDigiFlag += m_maxMdcDigi;
00612 if(m_dropHot) getDigiFlag |= MdcRawDataProvider::b_dropHot;
00613 if(m_keepBadTdc) getDigiFlag |= MdcRawDataProvider::b_keepBadTdc;
00614 if(m_keepUnmatch) getDigiFlag |= MdcRawDataProvider::b_keepUnmatch;
00615 MdcDigiVec mdcDigiVec = m_rawDataProviderSvc->getMdcDigiVec(getDigiFlag);
00616
00617 MdcDigiCol::iterator iter = mdcDigiVec.begin();
00618 int t_i =0;
00619 for (;iter != mdcDigiVec.end(); iter++,++t_i ) {
00620 double t = RawDataUtil::MdcTime((*iter)->getTimeChannel());
00621 double c = (*iter)->getChargeChannel();
00622 int l = MdcID::layer((*iter)->identify());
00623 int w = MdcID::wire((*iter)->identify());
00624 if(g_tupleEvt){
00625 m_na_t3layer[t_i] = l;
00626 m_na_t3wire[t_i] = w;
00627 m_na_t3gwire[t_i] = Constants::nWireBeforeLayer[l] + w;
00628 m_na_t3rt[t_i] = t;
00629 m_na_t3rtNot0[t_i] = t - t_t0;
00630 m_na_t3rc[t_i] = c;
00631 const MdcDigi* digi =const_cast<const MdcDigi*> (*iter);
00632 m_na_t3ovfl[t_i] = (*iter)->getOverflow();
00633 m_na_t3tNum[t_i] = ((*iter)->getOverflow()&0xF0)>>4;
00634 }
00635 }
00636 if(g_tupleEvt) m_na_t3nDigi = t_i;
00637
00638 m_na_t3TrkReco = t_trkRecoTk;
00639 m_na_t3Curl = t_curlTk;
00640 m_na_t3PatRec= t_patRecTk;
00641 m_na_t3XRec= t_xRecTk;
00642
00643 m_na_t3t0 = t_t0;
00644 m_na_t3t0Stat = t_t0Stat;
00645
00646 m_na_t3recTk = t_recTkNum;
00647 m_na_t3mcTk = t_mcTkNum;
00648
00649 m_na_t3evtNo = t_eventNo;
00650 m_na_t3runNo = t_runNo;
00651 g_tupleEvt->write();
00652
00653 return StatusCode::SUCCESS;
00654 }
00655
00656 double MdcNavigation::poca(const MdcDigi* aDigi,const HepVector helixPar,const HepSymMatrix errMat){
00657 int lay = MdcID::layer(aDigi->identify());
00658 int wire = MdcID::wire(aDigi->identify());
00659
00660 double ALPHA_loc,rho,pt,kappa,phiIn;
00661 int charge;
00662 double Bz = m_pIMF->getReferField()*1000.;
00663 ALPHA_loc = 333.567*Bz;
00664 charge = ( kappa >=0 ) ? 1 : -1 ;
00665 rho = ALPHA_loc/kappa ;
00666 pt = fabs( 1.0 / kappa);
00667 HepVector helix(helix);
00668 helix[0] = -helixPar[0];
00669 helix[1] = helixPar[1]+ pi/2;
00670 helix[2] = -1./helixPar[2];
00671 helix[3] = helixPar[3];
00672 helix[4] = helixPar[4];
00673 HelixTraj* m_helixTraj;
00674 MdcSagTraj* m_wireTraj_I;
00675 MdcSWire* m_mdcSWire_I;
00676 TrkPoca* m_trkPoca;
00677 TrkPoca* m_trkPoca_I;
00678 TrkPoca* m_trkPoca_O;
00679 MdcSagTraj* m_wireTraj_O;
00680 MdcSWire* m_mdcSWire_O;
00681 m_helixTraj = new HelixTraj(helix,errMat);
00682
00683 const MdcLayer* layPtr = m_gm->Layer(lay);
00684 double fltLenIn = layPtr->rMid();
00685 phiIn = helix[1] + helix[2]*fltLenIn;
00686 BesAngle tmp(phiIn - layPtr->phiEPOffset());
00687 int wlow = (int)floor(layPtr->nWires() * tmp.rad() / twopi );
00688 int wbig = (int)ceil(layPtr->nWires() * tmp.rad() / twopi );
00689 if (tmp == 0 ){ wlow = -1; wbig = 1; }
00690 int wireIn,wireOut;
00691 wireIn = wbig;
00692 wireOut = wlow;
00693 std::cout<<" in MdcNavigation lay/4 "<<lay/4<<" phi "<<tmp<<" wire "<<wireIn<<" "<<wireOut<<std::endl;
00694 }
00695
00696 StatusCode MdcNavigation::fillRecHits(RecMdcHitCol& hitCol){
00697 int ihit=0;
00698 RecMdcHitCol::iterator itHit = hitCol.begin();
00699 for(;itHit != hitCol.end(); itHit++ ) {
00700 RecMdcHit* h = *itHit;
00701 if ( !h ) continue;
00702 double m_na_lr = h->getFlagLR();
00703 double ddrift = -999;
00704 double ddoca = -999;
00705 ddoca = h->getDoca();
00706 if( 0 == m_na_lr ){ ddrift = h->getDriftDistLeft();
00707 }else{ ddrift = h->getDriftDistRight();}
00708 m_na_resid[ihit] = fabs(ddrift) - fabs(ddoca);
00709 if( 0 == m_na_lr ){ m_na_resid[ihit] *= -1.0;}
00710 m_na_driftD[ihit] = ddrift;
00711 m_na_driftT[ihit] = h->getDriftT();
00712 m_na_doca[ihit] = ddoca;
00713 m_na_entra[ihit] = h->getEntra();
00714 m_na_zhit[ihit] = h->getZhit();
00715 m_na_chi2add[ihit] = h->getChisqAdd();
00716 m_na_flaglr[ihit] = h->getFlagLR();
00717 m_na_hitStat[ihit] = h->getStat();
00718 m_na_Tdc[ihit] = h->getTdc();
00719 m_na_Adc[ihit] = h->getAdc();
00720 m_na_act[ihit] = h->getStat();
00721 int tlayer = MdcID::layer(h->getMdcId());
00722 int twire = MdcID::wire(h->getMdcId());
00723 m_na_layer[ihit] = tlayer;
00724 m_na_wire[ihit] = twire;
00725 m_na_gwire[ihit] = Constants::nWireBeforeLayer[tlayer] + twire;
00726 ++ihit;
00727 }
00728 m_na_nHit = ihit;
00729 return StatusCode::SUCCESS;
00730 }
00731
00732 double MdcNavigation::probab(const int& ndof, const double& chisq){
00733
00734
00735 static double srtopi=2.0/sqrt(2.0*M_PI);
00736 static double upl=100.0;
00737
00738 double prob=0.0;
00739 if(ndof<=0) {return prob;}
00740 if(chisq<0.0) {return prob;}
00741 if(ndof<=60) {
00742
00743 if(chisq>upl) {return prob;}
00744 double sum=exp(-0.5*chisq);
00745 double term=sum;
00746
00747 int m=ndof/2;
00748 if(2*m==ndof){
00749 if(m==1){return sum;}
00750 for(int i=2; i<=m;i++){
00751 term=0.5*term*chisq/(i-1);
00752 sum+=term;
00753 }
00754 return sum;
00755
00756
00757 }else{
00758
00759 double srty=sqrt(chisq);
00760 double temp=srty/M_SQRT2;
00761 prob=erfc(temp);
00762 if(ndof==1) {return prob;}
00763 if(ndof==3) {return (srtopi*srty*sum+prob);}
00764 m=m-1;
00765 for(int i=1; i<=m; i++){
00766 term=term*chisq/(2*i+1);
00767 sum+=term;
00768 }
00769 return (srtopi*srty*sum+prob);
00770
00771 }
00772
00773 }else{
00774
00775 double srty=sqrt(chisq)-sqrt(ndof-0.5);
00776 if(srty<12.0) {prob=0.5*erfc(srty);};
00777
00778 }
00779
00780 return prob;
00781 }