#include <MucRecRoadFinder.h>
Public Member Functions | |
| StatusCode | execute () |
| StatusCode | execute () |
| StatusCode | finalize () |
| StatusCode | finalize () |
| StatusCode | initialize () |
| StatusCode | initialize () |
| MucRecRoadFinder (const std::string &name, ISvcLocator *pSvcLocator) | |
| MucRecRoadFinder (const std::string &name, ISvcLocator *pSvcLocator) | |
| void | TrackFinding (RecMucTrack *aTrack) |
| void | TrackFinding (RecMucTrack *aTrack) |
| ~MucRecRoadFinder () | |
| ~MucRecRoadFinder () | |
Private Attributes | |
| MucRecHitContainer * | aMucRecHitContainer |
| MucRecHitContainer * | aMucRecHitContainer |
| NTuple::Item< double > | m_angle_cosmic |
| NTuple::Item< double > | m_angle_cosmic |
| NTuple::Item< double > | m_angle_updown |
| NTuple::Item< double > | m_angle_updown |
| std::string | m_configFile |
| NTuple::Item< double > | m_diff |
| NTuple::Item< double > | m_diff |
| NTuple::Item< double > | m_dist |
| NTuple::Item< double > | m_dist |
| NTuple::Item< double > | m_emcDown |
| NTuple::Item< double > | m_emcDown |
| NTuple::Item< double > | m_emcUp |
| NTuple::Item< double > | m_emcUp |
| NTuple::Item< double > | m_event |
| NTuple::Item< double > | m_event |
| int | m_fittingMethod |
| NTuple::Item< double > | m_gap |
| NTuple::Item< double > | m_gap |
| NTuple::Item< double > | m_maxhit |
| NTuple::Item< double > | m_maxhit |
| int | m_maxHitsRec |
| int | m_mccosmic |
| bool | m_MsOutput |
| NTuple::Item< double > | m_mucDown |
| NTuple::Item< double > | m_mucDown |
| NTuple::Item< double > | m_mucUp |
| NTuple::Item< double > | m_mucUp |
| NTuple::Item< double > | m_multihit |
| NTuple::Item< double > | m_multihit |
| int | m_NEvent |
| int | m_NEventReco |
| int | m_NEventWithHit |
| NTuple::Item< double > | m_ngapwithhits |
| NTuple::Item< double > | m_ngapwithhits |
| NTuple::Item< double > | m_nhit |
| NTuple::Item< double > | m_nhit |
| std::vector< int > | m_NHitsLost |
| std::vector< int > | m_NHitsLost |
| std::vector< int > | m_NHitsLostInGap |
| std::vector< int > | m_NHitsLostInGap |
| int | m_NHitsLostTotal |
| int | m_NHitsTotal |
| int | m_NtOutput |
| int | m_onlyseedfit |
| NTuple::Item< double > | m_part |
| NTuple::Item< double > | m_part |
| NTuple::Item< double > | m_phi |
| NTuple::Item< double > | m_phi |
| NTuple::Item< double > | m_phi_mc |
| NTuple::Item< double > | m_phi_mc |
| NTuple::Item< double > | m_phi_mc_pipe |
| NTuple::Item< double > | m_phi_mc_pipe |
| NTuple::Item< double > | m_phi_pipe |
| NTuple::Item< double > | m_phi_pipe |
| NTuple::Item< double > | m_projx |
| NTuple::Item< double > | m_projx |
| NTuple::Item< double > | m_projz |
| NTuple::Item< double > | m_projz |
| NTuple::Item< double > | m_px |
| NTuple::Item< double > | m_px |
| NTuple::Item< double > | m_px_mc |
| NTuple::Item< double > | m_px_mc |
| NTuple::Item< double > | m_py |
| NTuple::Item< double > | m_py |
| NTuple::Item< double > | m_py_mc |
| NTuple::Item< double > | m_py_mc |
| NTuple::Item< double > | m_pz |
| NTuple::Item< double > | m_pz |
| NTuple::Item< double > | m_pz_mc |
| NTuple::Item< double > | m_pz_mc |
| NTuple::Item< double > | m_run |
| NTuple::Item< double > | m_run |
| int | m_seedtype |
| NTuple::Item< double > | m_seg |
| NTuple::Item< double > | m_seg |
| NTuple::Item< double > | m_strip |
| NTuple::Item< double > | m_strip |
| NTuple::Item< double > | m_theta |
| NTuple::Item< double > | m_theta |
| NTuple::Item< double > | m_theta_mc |
| NTuple::Item< double > | m_theta_mc |
| NTuple::Item< double > | m_theta_mc_pipe |
| NTuple::Item< double > | m_theta_mc_pipe |
| NTuple::Item< double > | m_theta_pipe |
| NTuple::Item< double > | m_theta_pipe |
| NTuple::Tuple * | m_tuple |
| NTuple::Tuple * | m_tuple |
| int | m_united |
|
||||||||||||
|
00049 :
00050 Algorithm(name, pSvcLocator)
00051 {
00052 // Declare the properties
00053 declareProperty("FittingMethod", m_fittingMethod = 2);
00054 declareProperty("ConfigFile", m_configFile = "MucConfig.xml");
00055 declareProperty("McCosmic", m_mccosmic = 0);
00056 declareProperty("OnlySeedFit", m_onlyseedfit = 0);
00057 declareProperty("NtOutput", m_NtOutput = 0);
00058 declareProperty("MaxHitsRec", m_maxHitsRec = 200);
00059 declareProperty("United", m_united = 0);
00060 declareProperty("SeedType", m_seedtype = 0);
00061 declareProperty("MsOutput", m_MsOutput = false);
00062 }
|
|
|
00032 {};
|
|
||||||||||||
|
|
|
|
00032 {};
|
|
|
|
|
|
00147 {
00148
00149 MsgStream log(msgSvc(), name());
00150 log << MSG::INFO << "in execute()" << endreq;
00151
00152 // Part 1: Get the event header, print out event and run number
00153 int event, run;
00154
00155 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00156 if (!eventHeader) {
00157 log << MSG::FATAL << "Could not find Event Header" << endreq;
00158 return( StatusCode::FAILURE);
00159 }
00160 log << MSG::INFO << "Run: " << eventHeader->runNumber() << " Event: " << eventHeader->eventNumber() << endreq;
00161
00162 int digiId;
00163
00164 //Part 2: Retrieve MC truth
00165
00166 Hep3Vector cosmicMom;
00167 if(m_mccosmic==1){
00168 SmartDataPtr<McParticleCol> mcParticleCol(eventSvc(),"/Event/MC/McParticleCol");
00169 if (!mcParticleCol) {
00170 log << MSG::FATAL << "Could not find McParticle" << endreq;
00171 return( StatusCode::FAILURE);
00172 }
00173
00174 McParticleCol::iterator iter_mc = mcParticleCol->begin();
00175 for (;iter_mc != mcParticleCol->end(); iter_mc++, digiId++) {
00176 log << MSG::DEBUG
00177 << " particleId = " << (*iter_mc)->particleProperty()
00178 << endreq;
00179 int pid = (*iter_mc)->particleProperty();
00180 //cout<<"in mucroadfinder pid = "<<pid<<endl;
00181 if(fabs(pid)!=13) continue;
00182
00183 HepLorentzVector initialMomentum = (*iter_mc)->initialFourMomentum();
00184 Hep3Vector startP(initialMomentum.px(), initialMomentum.py(), initialMomentum.pz());
00185 Hep3Vector rotate(-initialMomentum.px(),initialMomentum.pz(),initialMomentum.py());
00186
00187 if(m_NtOutput>=1){
00188 m_px_mc = initialMomentum.px();
00189 m_py_mc = initialMomentum.py();
00190 m_pz_mc = initialMomentum.pz();
00191
00192 m_theta_mc = rotate.theta();
00193 m_phi_mc = rotate.phi();
00194
00195 m_theta_mc_pipe = startP.theta();
00196 m_phi_mc_pipe = startP.phi();
00197
00198 //m_tuple->write();
00199 }
00200
00201 //if(initialMomentum.py()>0)cout<<"py>0????? "<<pid<<endl;
00202 cosmicMom = startP;
00203
00204 }
00205 }
00206 /*
00207 //Part 3: Retrieve MDC digi
00208 SmartDataPtr<MdcDigiCol> mdcDigiCol(eventSvc(),"/Event/Digi/MdcDigiCol");
00209 if (!mdcDigiCol) {
00210 log << MSG::FATAL << "Could not find MDC digi" << endreq;
00211 return( StatusCode::FAILURE);
00212 }
00213
00214 MdcDigiCol::iterator iter1 = mdcDigiCol->begin();
00215 digiId = 0;
00216
00217 for (;iter1 != mdcDigiCol->end(); iter1++, digiId++) {
00218 log << MSG::INFO << "MDC digit No: " << digiId << endreq;
00219
00220 log << MSG::INFO
00221 << " time_channel = " << (*iter1)->getTimeChannel()
00222 << " charge_channel = " << (*iter1)->getChargeChannel()
00223 << endreq;
00224 }
00225
00226 //Part 4: Retrieve TOF digi
00227 SmartDataPtr<TofDigiCol> tofDigiCol(eventSvc(),"/Event/Digi/TofDigiCol");
00228 if (!tofDigiCol) {
00229 log << MSG::FATAL << "Could not find TOF digi" << endreq;
00230 return( StatusCode::FAILURE);
00231 }
00232
00233 TofDigiCol::iterator iter2 = tofDigiCol->begin();
00234 digiId = 0;
00235
00236 for (;iter2 != tofDigiCol->end(); iter2++, digiId++) {
00237 log << MSG::INFO << "TOF digit No: " << digiId << endreq;
00238
00239 }
00240
00241 //Part 5: Retrieve EMC digi
00242 SmartDataPtr<EmcDigiCol> emcDigiCol(eventSvc(),"/Event/Digi/EmcDigiCol");
00243 if (!emcDigiCol) {
00244 log << MSG::FATAL << "Could not find EMC digi" << endreq;
00245 return( StatusCode::FAILURE);
00246 }
00247
00248 EmcDigiCol::iterator iter3 = emcDigiCol->begin();
00249 digiId = 0;
00250
00251 for (;iter3 != emcDigiCol->end(); iter3++, digiId++) {
00252 log << MSG::INFO << "Emc digit No: " << digiId << endreq;
00253
00254 log << MSG::INFO
00255 << " time_channel = " << (*iter3)->getTimeChannel()
00256 << " charge_channel = " << (*iter3)->getChargeChannel()
00257 << endreq;
00258 }
00259 */
00260
00261 //Part 6: Retrieve MUC digi
00262 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00263 if (!mucDigiCol) {
00264 log << MSG::FATAL << "Could not find MUC digi" << endreq;
00265 return( StatusCode::FAILURE);
00266 }
00267
00268 MucDigiCol::iterator iter4 = mucDigiCol->begin();
00269 digiId = 0;
00270 for (;iter4 != mucDigiCol->end(); iter4++, digiId++) {
00271 }
00272 log << MSG::INFO << "MUC digis:: " << digiId << endreq;
00273 if( digiId == 0) return( StatusCode::SUCCESS );
00274
00275 //********************************//
00276 RecMucTrackCol *aMucTrackCol;
00277 int trackId = -1;
00278 int muctrackId = -1;
00279
00280 if(m_united != 1) // if this algorithm run individualy, we need create hitcollection and trackcollection...
00281 {
00282 Identifier mucID;
00283 if (!aMucRecHitContainer) {
00284 aMucRecHitContainer = new MucRecHitContainer();
00285 }
00286 aMucRecHitContainer->Clear();
00287
00288 MucRecHitCol *aMucRecHitCol = new MucRecHitCol();
00289 aMucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00290
00291 IDataManagerSvc *dataManSvc = dynamic_cast<IDataManagerSvc*> (eventSvc());
00292 DataObject *mucRecHitCol;
00293 eventSvc()->findObject("/Event/Recon/MucRecHitCol",mucRecHitCol);
00294 if(mucRecHitCol != NULL) {
00295 dataManSvc->clearSubTree("/Event/Recon/MucRecHitCol");
00296 eventSvc()->unregisterObject("/Event/Recon/MucRecHitCol");
00297 }
00298
00299 StatusCode sc = eventSvc()->registerObject("/Event/Recon/MucRecHitCol", aMucRecHitContainer->GetMucRecHitCol());
00300
00301 MucDigiCol::iterator iterMuc = mucDigiCol->begin();
00302 int mucDigiId = 0;
00303 for (;iterMuc != mucDigiCol->end(); iterMuc++, mucDigiId++) {
00304 mucID = (*iterMuc)->identify();
00305 int part = MucID::barrel_ec(mucID);
00306 int seg = MucID::segment(mucID);
00307 int gap = MucID::layer(mucID);
00308 int strip = MucID::channel(mucID);
00309 //aMucRecHitContainer->AddHit(mucID);
00310 //if(part==1 && seg==2 && gap==8 && (strip==19||strip==16));//cout<<"noise hit!!!=========="<<endl;
00311 //else /////this problem has been solved!
00312 aMucRecHitContainer->AddHit(part, seg, gap, strip);
00313 log << MSG::INFO << " digit" << mucDigiId << " : "
00314 << " part " << part
00315 << " seg " << seg
00316 << " gap " << gap
00317 << " strip " << strip
00318 << endreq;
00319 }
00320
00321 DataObject *aReconEvent ;
00322 eventSvc()->findObject("/Event/Recon",aReconEvent);
00323 if(aReconEvent==NULL) {
00324 //then register Recon
00325 aReconEvent = new ReconEvent();
00326 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00327 if(sc!=StatusCode::SUCCESS) {
00328 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00329 return( StatusCode::FAILURE);
00330 }
00331 }
00332 StatusCode fr = eventSvc()->findObject("/Event/Recon",aReconEvent);
00333 if(fr!=StatusCode::SUCCESS) {
00334 log << MSG::WARNING << "Could not find register ReconEvent, will create it" <<endreq;
00335 StatusCode sc = eventSvc()->registerObject("/Event/Recon",aReconEvent);
00336 if(sc!=StatusCode::SUCCESS) {
00337 log << MSG::FATAL << "Could not register ReconEvent" <<endreq;
00338 return( StatusCode::FAILURE);
00339 }
00340 }
00341
00342 //********************************//
00343 // Create track Container;
00344 DataObject *mucTrackCol;
00345 eventSvc()->findObject("/Event/Recon/RecMucTrackCol",mucTrackCol);
00346 if(mucTrackCol != NULL) {
00347 dataManSvc->clearSubTree("/Event/Recon/RecMucTrackCol");
00348 eventSvc()->unregisterObject("/Event/Recon/RecMucTrackCol");
00349 }
00350
00351 aMucTrackCol = new RecMucTrackCol();
00352 sc = eventSvc()->registerObject("/Event/Recon/RecMucTrackCol", aMucTrackCol);
00353 aMucTrackCol->clear();
00354
00355 // check MucTrackCol registered
00356 SmartDataPtr<RecMucTrackCol> findMucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
00357 if (!findMucTrackCol) {
00358 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
00359 return( StatusCode::FAILURE);
00360 }
00361 aMucTrackCol->clear();
00362
00363
00364 log << MSG::INFO <<"MaxHitsRec : "<<m_maxHitsRec<< endreq;
00365 if(digiId> m_maxHitsRec) return StatusCode::SUCCESS; //too many digit! abnormal------2008-04-17
00366 //********************************//
00367 // Create track Container;
00368
00369 trackId = 0;
00370 muctrackId = 0;
00371 }
00372 else if(m_united == 1){
00373
00374 Identifier mucID;
00375 if (!aMucRecHitContainer) {
00376 aMucRecHitContainer = new MucRecHitContainer();
00377 }
00378 aMucRecHitContainer->Clear();
00379
00380 SmartDataPtr<MucRecHitCol> aMucRecHitCol (eventSvc(),"/Event/Recon/MucRecHitCol");
00381 if(aMucRecHitCol == NULL) {
00382 log << MSG::WARNING << "MucRecHitCol is NULL" << endreq;
00383 }
00384
00385 aMucRecHitContainer->SetMucRecHitCol(aMucRecHitCol);
00386
00387 SmartDataPtr<RecMucTrackCol> mucTrackCol (eventSvc(),"/Event/Recon/RecMucTrackCol");
00388 if(mucTrackCol == NULL) {
00389 log << MSG::WARNING << "aMucTrackCol is NULL" << endreq;
00390 }
00391
00392 log << MSG::INFO <<"mucTrackCol size: "<<mucTrackCol->size()<<" hitCol size: "<<aMucRecHitCol->size()<<endreq;
00393 aMucTrackCol = mucTrackCol;
00394
00395 // Retrieve Ext track Col from TDS for setting trackId
00396 SmartDataPtr<RecExtTrackCol> aExtTrackCol(eventSvc(),"/Event/Recon/RecExtTrackCol");
00397 if (!aExtTrackCol) {
00398 log << MSG::WARNING << "Can't find ExtTrackCol in TDS!" << endreq;
00399 }
00400
00401 SmartDataPtr<RecMdcTrackCol> aMdcTrackCol(eventSvc(),"/Event/Recon/RecMdcTrackCol");
00402 if (!aMdcTrackCol) {
00403 log << MSG::WARNING << "Can't find MdcTrackCol in TDS!" << endreq;
00404 }
00405
00406 if( aExtTrackCol && aMdcTrackCol ) trackId = aMdcTrackCol->size();
00407 muctrackId = aMucTrackCol->size();
00408 }
00409
00410 int hasEmcUp = 0;
00411 int hasEmcDown = 0;
00412 SmartDataPtr<RecEmcShowerCol> aShowerCol(eventSvc(),"/Event/Recon/RecEmcShowerCol");
00413 if (!aShowerCol) {
00414 log << MSG::WARNING << "Could not find RecEmcShowerCol" << endreq;
00415 //return( StatusCode::FAILURE);
00416 }
00417 else{
00418 RecEmcShowerCol::iterator iShowerCol;
00419 for(iShowerCol=aShowerCol->begin();
00420 iShowerCol!=aShowerCol->end();
00421 iShowerCol++){
00422 /*
00423 cout<<"check EmcRecShowerCol registered:"<<endl;
00424 cout<<"shower id: "<<(*iShowerCol)->ShowerId()<<"\t"
00425 <<"shower energy: "<<(*iShowerCol)->Energy()<<"\n"
00426 <<"shower position: "<<(*iShowerCol)->Position()<<endl;
00427 */
00428 if(((*iShowerCol)->position()).phi()>0&&((*iShowerCol)->position()).phi()<3.1415926) hasEmcUp++;
00429 else hasEmcDown++;
00430 }
00431 }
00432 if(m_NtOutput >= 1){
00433 m_emcUp = hasEmcUp; m_emcDown = hasEmcDown;
00434 }
00435
00436 //********************************//
00437 m_NEvent++;
00438 MucRec3DRoadContainer *p3DRoadC = new MucRec3DRoadContainer();
00439 // Find or create the 3D road container ...
00440 if (!p3DRoadC) {
00441 p3DRoadC = new MucRec3DRoadContainer();
00442 }
00443 p3DRoadC->clear(); // make sure that the container is empty
00444
00445 MucRec2DRoadContainer *p2DRoad0C = new MucRec2DRoadContainer();
00446 // Find or create the 2D road0 container ...
00447 if (!p2DRoad0C) {
00448 p2DRoad0C = new MucRec2DRoadContainer();
00449 }
00450 p2DRoad0C->clear(); // make sure that the container is empty
00451
00452 MucRec2DRoadContainer *p2DRoad1C = new MucRec2DRoadContainer();
00453 // Find or create the 2D road1 container ...
00454 if (!p2DRoad1C) {
00455 p2DRoad1C = new MucRec2DRoadContainer();
00456 }
00457 p2DRoad1C->clear(); // make sure that the container is empty
00458 log << MSG::INFO << "Muc 2D 3D Road Container created" << endreq;
00459
00460 // We have all of the input and output parameters now, so do
00461 // something useful with them!
00462 //static bool first_call = true;
00463 //
00464 // Now find some roads!
00465 //
00466
00467 MucRecHit *pHit0 = 0, *pHit1 = 0;
00468 int count0, count1, count, iGap0, iGap1;
00469 int orient;
00470
00471 for (int iPart = 0; iPart < (int)MucID::getPartNum(); iPart++) {
00472 for (int iSeg = 0; iSeg < (int)MucID::getSegNum(iPart); iSeg++) {
00473 for (int iOrient = 0; iOrient < 2; iOrient++) {
00474 int nLoops = 1;
00475 if(m_seedtype == 1) nLoops = kNSeedLoops;
00476 for (int iLoop = 0; iLoop < nLoops; iLoop++) { //now only 1 loop
00477 // checkout the seed gap(s) from search order
00478 int length = kSearchLength[iLoop][iPart][iOrient];
00479 if(m_seedtype == 1){
00480 iGap0 = kSearchOrder[iLoop][iPart][iOrient][0];
00481 iGap1 = kSearchOrder[iLoop][iPart][iOrient][1];
00482 }
00483 else { //new method to calc seed gaps------//
00484 int setgap0 = 0;
00485 MucDigiCol::iterator iterMuc = mucDigiCol->begin();
00486 int mucDigiId = 0;
00487 Identifier mucID;
00488 iGap0 = 0; iGap1 = 0;
00489 for (;iterMuc != mucDigiCol->end(); iterMuc++, mucDigiId++) {
00490 mucID = (*iterMuc)->identify();
00491 int part = MucID::barrel_ec(mucID);
00492 int seg = MucID::segment(mucID);
00493 int gap = MucID::layer(mucID);
00494 int strip = MucID::channel(mucID);
00495 int orient = 0;
00496 if ( (part == 1 && gap % 2 == 0) || (part != 1 && gap % 2 == 1) ) orient = 1;
00497
00498 if(part == iPart && seg == iSeg && orient == iOrient && setgap0 == 0){
00499 iGap0 = gap;
00500 setgap0 = 1;
00501 }//make sure that iGap0 record the 1st gap in this orient
00502 if(part == iPart && seg == iSeg && orient == iOrient && setgap0 == 1 && gap != iGap0 ){
00503 iGap1 = gap;
00504 }//make sure that iGap1 record the last gap in this orient
00505 }
00506 }
00507 //
00508 if(m_MsOutput) cout <<"Find seed gap in ( "<<iPart<<" "<<iSeg<<" ) gap0: "<<iGap0<<" gap1: "<<iGap1<<endl;
00509
00510 //new method to calc seed gaps------//
00511 if(iGap0 > iGap1){
00512 int tempGap = iGap0;
00513 iGap0 = iGap1;
00514 iGap1 = tempGap;
00515 }
00516 if(iGap1 == iGap0) continue;
00517
00518 count0 = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap0);
00519 for (int iHit0 = 0; iHit0 < count0; iHit0++) {
00520 //cout << "iHit0 = " << iHit0 << endl;
00521 pHit0 = aMucRecHitContainer->GetHit(iPart, iSeg, iGap0, iHit0);
00522 if (!pHit0) {
00523 log << MSG::FATAL << "MucRecRoadFinder-E10 " << " part " << iPart
00524 << " seg " << iSeg << " gap " << iGap0 << " null pointer"
00525 << endl;
00526 }
00527 else {
00528
00529 //this algo use with ext and this hit has been used before;
00530 if(m_united == 1 && pHit0->GetHitMode() != -1) continue;
00531
00532 count1 = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap1);
00533 if(m_MsOutput) cout << "part " << iPart << " seg " << iSeg << " orient " << iOrient
00534 << " gap0 " << iGap0 << " count0 " << count0
00535 << " gap1 " << iGap1 << " count1 " << count1 << endl;
00536 for (int iHit1 = 0; iHit1 < count1; iHit1++) {
00537 //cout << "iHit1 = " << iHit1 << endl;
00538 pHit1 = aMucRecHitContainer->GetHit(iPart, iSeg, iGap1, iHit1);
00539 if (!pHit1) {
00540 log << MSG::ERROR << "MucRecRoadFinder-E10 " << " part " << iPart
00541 << " seg " << iSeg << " gap " << iGap1 << " null pointer"
00542 << endl;
00543 } else {
00544
00545 //this algo use with ext and this hit has been used before;
00546 if(m_united == 1 && pHit1->GetHitMode() != -1) continue;
00547
00548 // Found seed hit(s), create a new 2D road object,
00549 // and attach hit pHit1 and pHit0
00550 MucRec2DRoad *road2D = new MucRec2DRoad(iPart, iSeg, iOrient, 0.0, 0.0, 3000.0);
00551 if (!road2D) {
00552 log << MSG::FATAL << "MucRecRoadFinder-E20 " << "failed to create 2D road!" << endreq;
00553 continue;
00554 }
00555 road2D->SetMaxNSkippedGaps(kMaxSkippedGaps);
00556
00557 if (!pHit0->GetGap()) log << MSG::ERROR << "pHit0->GetGap(), null pointer" << endreq;
00558 if (pHit0->GetGap()->Orient() == iOrient) {
00559 //float xStrip, yStrip, zStrip;
00560 //MucRecGeometry::GetPointer()->GetStripPointer(pHit0->GetSoftID())->GetCenterPos(xStrip, yStrip, zStrip);
00561 //cout << " x = " << xStrip << " y = " << yStrip << " z = " << zStrip << endl;
00562
00563 //set hit as seed then
00564 bool isseed = true;
00565 pHit0->SetHitSeed(true);
00566 pHit1->SetHitSeed(true);
00567
00568 road2D->AttachHit(pHit0);
00569 if(m_MsOutput) cout << "pHit0 attached, " << road2D->GetTotalHits()
00570 << ", hitId= "<<pHit0->Part()<<" "<<pHit0->Seg()<<" "<<pHit0->Gap()<<" "<<pHit0->Strip()<<endl;
00571 }
00572
00573 if (pHit1->GetGap()->Orient() == iOrient) {
00574 //cout << pHit1->GetSoftID() << endl;
00575 //float xStrip, yStrip, zStrip;
00576 //MucRecGeometry::GetPointer()->GetStripPointer(pHit1->GetSoftID())->GetCenterPos(xStrip, yStrip, zStrip);
00577 //cout << " x = " << xStrip << " y = " << yStrip << " z = " << zStrip << endl;
00578 road2D->AttachHit(pHit1);
00579 if(m_MsOutput) cout << "pHit1 attached, " << road2D->GetTotalHits()
00580 << ", hitId= "<<pHit1->Part()<<" "<<pHit1->Seg()<<" "<<pHit1->Gap()<<" "<<pHit1->Strip()<<endl;
00581 }
00582 if(m_MsOutput) cout << "Hit cout " << road2D->GetTotalHits() << ", slope " << road2D->GetSlope() << endl;
00583
00584 // After first (two) hit(s) is(are) attached,
00585 // calculate reference pos, direction.
00586 // Project to the other planes;
00587 // attach clusters within search window.
00588
00589 for (int i = 0; i < length; i++) {
00590 int iGap = kSearchOrder[iLoop][iPart][iOrient][i];
00591
00592 float dXmin = kInfinity;
00593 MucRecHit *pHit = 0;
00594
00595 float Win = kWindowSize[iPart][iGap];
00596 //Win = road2D->GetSearchWindowSize(iGap);
00597
00598 // following line should be removed after
00599 // we have a good function classes for search window
00600 //Win = WindowParGamma[iGap];
00601
00602 //search hit in iGap
00603 count = aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap);
00604 for (int iHit = 0; iHit < count; iHit++) {
00605 pHit = aMucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
00606 if (!pHit) {
00607 log << MSG::FATAL << "MucRecRoadFinder-E20 null cluster pointer pHit" << endreq;
00608 continue;
00609 }
00610 else {
00611
00612 //this algo use with ext and this hit has been used before;
00613 if(m_united == 1 && pHit->GetHitMode() != -1) continue;
00614
00615 // Get the displacement of hit pHit to road2D
00616 float dX = road2D->GetHitDistance(pHit);
00617 if(m_MsOutput) cout<<"dX = "<<dX<<" Win = "<<Win<<", hit: "<<pHit->Part()<<" "<<pHit->Seg()<<" "<<pHit->Gap()<<" "<<pHit->Strip()<<endl;
00618
00619 if (dX < Win) {
00620 // Attach the hit if it exists
00621 // cout << "meet window " << pHit->GetSoftID() << endl;
00622 if(iGap == iGap0 || iGap == iGap1) { //hits in these gap is used to be seed, so unvailable for calibrate effi;
00623 if(pHit->GetHitMode() == -1){
00624 pHit->SetHitMode(3); //self road finder mode
00625 }
00626 else if(pHit->GetHitMode() == 3){
00627 pHit->SetHitMode(4); //self road finder mode used!!!
00628 }
00629 }
00630 // attach hits if we need fit or hit. OR, only attach hits not in seed layers
00631 if(m_onlyseedfit == 0)road2D->AttachHit(pHit);
00632 else {
00633 if(iGap == iGap0 || iGap == iGap1) road2D->AttachHit(pHit);
00634 else road2D->AttachHitNoFit(pHit);
00635 }
00636 }
00637 }
00638 }
00639
00640
00641 } // for (int i = 0; i < length; ...), Search all gaps
00642
00643
00644 // Ok we found a road already, we need to know
00645 // whether we should save it or not
00646 // check road quality...
00647 //
00648
00649 // 1. last gap of the road
00650 bool lastGapOK = false;
00651 if ( road2D->GetLastGap() >= kMinLastGap2D ) {
00652 lastGapOK = true;
00653 }
00654 else {
00655 log<<MSG::INFO << " 2D lastGap " << road2D->GetLastGap() << " < " << kMinLastGap2D << endreq;
00656 }
00657
00658 // 2. positon at reference plane
00659
00660 // 3. number of gaps contain hits
00661 bool firedGapsOK = false;
00662 if (road2D->GetNGapsWithHits() >= kMinFiredGaps) {
00663 firedGapsOK = true;
00664 }
00665 else {
00666 log<<MSG::INFO << " 2D firedGaps " << road2D->GetNGapsWithHits() << " < " << kMinFiredGaps << endreq;
00667 }
00668
00669 // 4. duplicated road check
00670 bool duplicateRoad = false;
00671 int maxSharedHits = 0, sharedHits = 0;
00672 if (iOrient == 0) {
00673 for (int index = 0; index < (int)p2DRoad0C->size(); index++) {
00674 sharedHits =(*p2DRoad0C)[index]->GetNSharedHits(road2D);
00675 }
00676 if (maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00677 }
00678 else {
00679 for (int index = 0; index < (int)p2DRoad1C->size(); index++) {
00680 sharedHits = (*p2DRoad1C)[index]->GetNSharedHits(road2D);
00681 }
00682 if (maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00683 }
00684
00685 if(road2D->GetTotalHits() == maxSharedHits
00686 || maxSharedHits >= kMinSharedHits2D ) {
00687 duplicateRoad = true;
00688 log<<MSG::INFO << " maxSharedHits " << maxSharedHits << " > " << kMinSharedHits2D << endreq;
00689 }
00690
00691 // Here is our criteria for a candidate road
00692 // i.e. 1) There are at least two gaps containing hits
00693 // 2) LastGap of the road passes last plane cut
00694 // 3) Reference position passes vertex cut
00695 // 4) No existing road share same hits with the road
00696
00697 if (lastGapOK && firedGapsOK && !duplicateRoad) {
00698 if (iOrient == 0) {
00699 log<<MSG::INFO << "Add a road0" << endreq;
00700 p2DRoad0C->add(road2D);
00701 }
00702 else {
00703 log<<MSG::INFO << "Add a road1" << endreq;
00704 p2DRoad1C->add(road2D);
00705 }
00706 }
00707 else {
00708
00709 // reset hit mode to be -1 if this road useless!
00710 vector<MucRecHit*> road2DHits = road2D->GetHits();
00711 for(int i=0; i< road2DHits.size(); i++)
00712 {
00713 MucRecHit *ihit = road2DHits[i];
00714 if(ihit->Gap() == iGap0 || ihit->Gap() == iGap1){
00715 if(ihit->GetHitMode() == 3) ihit->SetHitMode(-1);
00716 if(ihit->GetHitMode() == 4) ihit->SetHitMode(3);
00717 }
00718 }
00719 delete road2D;
00720 }
00721 } // else {
00722 } // for ( int iHit = 0 ...)
00723 } // else {
00724 } // for ( int iHit0 ...)
00725 } // for (int iLoop ...)
00726 } // for (int iOrient ...)
00727 } // for (int iSeg ...)
00728 } // for(iPartArm ...)
00729
00730 int print2DRoadInfo = 0;
00731 if (print2DRoadInfo == 1) {
00732 // print 2DRoad container info.
00733 cout << "In 2DRoad container : " << endl
00734 << " Num of 2DRoad0 = " << p2DRoad0C->size() << endl
00735 << endl;
00736
00737 for (int iRoad = 0; iRoad < (int)p2DRoad0C->size(); iRoad++) {
00738
00739 MucRec2DRoad *road = (*p2DRoad0C)[iRoad];
00740 cout << " " << iRoad << "th 2DRoad0 :" << endl
00741 << " Part = " << road->GetPart() << endl
00742 << " Seg = " << road->GetSeg() << endl
00743 << " Orient = " << road->GetOrient() << endl
00744 << " LastGap = " << road->GetLastGap() << endl
00745 << " TotalHits = " << road->GetTotalHits() << endl
00746 << " DOF = " << road->GetDegreesOfFreedom() << endl
00747 << " Slope = " << road->GetSlope() << endl
00748 << " Intercept = " << road->GetIntercept() << endl
00749 << endl;
00750 //road->PrintHitsInfo();
00751 }
00752
00753 cout << " Num of 2DRoad1 = " << p2DRoad1C->size() << endl
00754 << endl;
00755
00756 for ( int iRoad = 0; iRoad < (int)p2DRoad1C->size(); iRoad++) {
00757
00758 MucRec2DRoad *road = (*p2DRoad1C)[iRoad];
00759 cout << " " << iRoad << "th 2DRoad1 :" << endl
00760 << " Part = " << road->GetPart() << endl
00761 << " Seg = " << road->GetSeg() << endl
00762 << " Orient = " << road->GetOrient() << endl
00763 << " LastGap = " << road->GetLastGap() << endl
00764 << " TotalHits = " << road->GetTotalHits() << endl
00765 << " DOF = " << road->GetDegreesOfFreedom() << endl
00766 << " Slope = " << road->GetSlope() << endl
00767 << " Intercept = " << road->GetIntercept() << endl
00768 << endl;
00769 //road->PrintHitsInfo();
00770 }
00771 }
00772
00773 // We found all possible 2D roads in one part and seg
00774 // and now it is time to construct 3D roads base on those 2D roads
00775 for ( int iRoad0 = 0; iRoad0 < (int)p2DRoad0C->size(); iRoad0++) {
00776 MucRec2DRoad *road0 = (*p2DRoad0C)[iRoad0];
00777 for ( int iRoad1 = 0; iRoad1 < (int)p2DRoad1C->size(); iRoad1++){
00778 MucRec2DRoad *road1 = (*p2DRoad1C)[iRoad1];
00779
00780 // Create a new road object with road0 and road1
00781 if ( !(road0 &&road1) ) {
00782 cout << "Null pointer to road0 or road1: "
00783 << "road0 = " << road0
00784 << "road1 = " << road1
00785 << endl;
00786 continue;
00787 }
00788
00789 // Check that both 1D roads are in the same part and segment.
00790 if ( (road0->GetPart() != road1->GetPart())
00791 || (road0->GetSeg() != road1->GetSeg()) ) {
00792 continue;
00793 }
00794
00795 MucRec3DRoad *road3D = new MucRec3DRoad(road0, road1);
00796
00797 // What are our criteria for accepting a road as valid?
00798 // Check them here.
00799 // 1. difference of total number clusters in road0 and in road1
00800 // 2. difference of last gaps in road0 and road1
00801 // 3. InterSection OK or not (if X Y clusters in each gap are
00802 // in same panel or overlap region)
00803 // 4. Reference position check
00804 // 5. Chisquare check
00805 // 6. Last Gap check
00806 // 7. Duplicate road check
00807
00808 bool lastGapDeltaOK = false;
00809 if ( road3D->GetLastGapDelta() <= kMaxDeltaLastGap ) {
00810 lastGapDeltaOK = true;
00811 }
00812 else {
00813 cout << "LastGapDelta = " << road3D->GetLastGapDelta() << endl;
00814 }
00815
00816 bool totalHitsDeltaOK = false;
00817 if ( road3D->GetTotalHitsDelta() <= kMaxDeltaTotalHits ) {
00818 totalHitsDeltaOK = true;
00819 }
00820 else {
00821 //cout << "TotalHitsDelta = " << road3D->GetTotalHitsDelta() << endl;
00822 }
00823
00824 bool chiSquareCutOK = false;
00825 float xChiSquare = road0->GetReducedChiSquare();
00826 float yChiSquare = road1->GetReducedChiSquare();
00827 if ( xChiSquare <= kMaxXChisq
00828 && yChiSquare <= kMaxYChisq ) {
00829 chiSquareCutOK = true;
00830 }
00831 else {
00832 cout << "xChiSquare = " << xChiSquare
00833 << "yChiSquare = " << yChiSquare
00834 << endl;
00835 }
00836
00837 bool minLastGapOK = false;
00838 if ( road3D->GetLastGap() >= kMinLastGap3D ) {
00839 minLastGapOK = true;
00840 }
00841 else {
00842 log<<MSG::INFO << " minLastGap = " << road3D->GetLastGap() << endreq;
00843 }
00844
00845 bool duplicateRoad = false;
00846 int maxSharedHits = 0, sharedHits = 0;
00847 for ( int i = 0; i < (int)p3DRoadC->size(); i++){
00848 sharedHits =(*p3DRoadC)[i]->GetNSharedHits(road3D);
00849 if ( maxSharedHits < sharedHits) maxSharedHits = sharedHits;
00850
00851 //cout<<"judge shared hits: road["<<i<<"] max = "<<maxSharedHits<<endl;
00852 }
00853 if(road3D->GetTotalHits() == maxSharedHits
00854 || maxSharedHits >= kMinSharedHits2D ) {
00855 duplicateRoad = true;
00856 log<<MSG::INFO << " MaxShareHits = " << maxSharedHits << endreq;
00857 }
00858
00859 if ( lastGapDeltaOK &&
00860 totalHitsDeltaOK &&
00861 chiSquareCutOK &&
00862 minLastGapOK &&
00863 !duplicateRoad ) {
00864 float vx, vy, x0, y0;
00865 float slope1 = 100, slope0 = 100;
00866 if(fabs(road1->GetSlope())<100) slope1 = road1->GetSlope();
00867 if(fabs(road0->GetSlope())<100) slope0 = road0->GetSlope();
00868
00869 if ( road3D->GetPart() == 1 ){
00870 road3D->TransformPhiRToXY( slope1, slope0,
00871 road1->GetIntercept(), road0->GetIntercept(),
00872 vx, vy, x0, y0);
00873 }
00874 else {
00875 vx = slope1;
00876 x0 = road1->GetIntercept();
00877 vy = slope0;
00878 y0 = road0->GetIntercept();
00879 }
00880 road3D->SetSimpleFitParams(vx, x0, vy, y0);
00881
00882 log<<MSG::INFO << "Add a 3D Road ... " << endreq;
00883
00884 float startx = 0.0, starty = 0.0, startz = 0.0;
00885 float sigmax = 0.0, sigmay = 0.0, sigmaz = 0.0;
00886 road3D->ProjectWithSigma(0, startx, starty, startz, sigmax, sigmay, sigmaz);//gap0
00887
00888 //cout<<"slope1,0 = "<<slope1<<" "<<slope0<<" vx,y = "<<vx<<" "<<vy<<endl;
00889 //cout<<"startxyz= "<<startx<<" "<<starty<<" "<<startz<<endl;
00890 //mom(vx,vy,1)
00891 float vz = 1;
00892 float sign = vy/fabs(vy);
00893 float signx = vx/fabs(vx);
00894 //cout<<"vxyz = "<<vx<<" "<<vy<<" "<<vz<<endl;
00895 if(road3D->GetPart() == 1){
00896 if(road3D->GetSeg() > 4){ //down segment
00897
00898 vx *= -sign;
00899 vy *= -sign;
00900 vz *= -sign;
00901 }
00902 else if(road3D->GetSeg()<2){
00903 vx *= signx;
00904 vy *= signx;
00905 vz *= signx;
00906 }
00907 else if(road3D->GetSeg()>=3 && road3D->GetSeg()<=4){
00908 vx *= -signx;
00909 vy *= -signx;
00910 vz *= -signx;
00911 }
00912 else{
00913 vx *= sign;
00914 vy *= sign;
00915 vz *= sign;
00916 }
00917 }
00918 else if(road3D->GetPart() == 0){
00919 //fix me
00920
00921 //cout<<"startxyz= "<<startx<<" "<<starty<<" "<<startz<<endl;
00922 //cout<<"vx,y,z = "<<vx<<" "<<vy<<" "<<vz<<endl;
00923 //cout<<"in road finder a endcap finded!!! -------------"<<endl;
00924 }
00925 else if(road3D->GetPart() == 2){
00926 //fix me
00927
00928 vx *= -1;
00929 vy *= -1;
00930 vz *= -1;
00931 //cout<<"startxyz= "<<startx<<" "<<starty<<" "<<startz<<endl;
00932 //cout<<"vx,y,z = "<<vx<<" "<<vy<<" "<<vz<<endl;
00933 //cout<<"in road finder a endcap finded!!! ------------2-"<<endl;
00934
00935 }
00936
00937
00938 Hep3Vector mom(vx, vy, vz);
00939
00941 //MucTrack *aTrack = new MucTrack(road3D);
00942 //cout<<"startxyz = "<<startx<<" "<<starty<<" "<<startz<<endl;
00943 //cout<<"vxyz = "<<vx<<" "<<vy<<" "<<vz<<endl;
00944
00945 startx /= 10; starty /= 10; startz /= 10; //mm->cm
00946 startx -= vx/mom.mag(); starty -= vy/mom.mag(); startz -= vz/mom.mag(); // decrease a little
00947
00948 //cout<<"startxyz = "<<startx<<" "<<starty<<" "<<startz<<endl;
00949 RecMucTrack *aTrack = new RecMucTrack();
00950 aTrack->SetExtMucPos(startx, starty, startz); //mm->cm
00951 aTrack->SetExtMucMomentum(vx, vy, vz);
00952 aTrack->SetMucPos(startx, starty, startz);
00953 aTrack->SetMucMomentum(vx, vy, vz);
00954 aTrack->SetCurrentPos( startx, starty, startz);
00955 aTrack->SetCurrentDir( vx, vy, vz);
00956 aTrack->SetRecMode(3);
00957 //aTrack->LineFit(1);
00958 //aTrack->ComputeTrackInfo(1);
00959
00960 aTrack->setTrackId(trackId);
00961 aTrack->setId(muctrackId);
00962 trackId++;
00963 muctrackId++;
00964
00965 //cout<<"find a track!!!"<<endl;
00966 aMucTrackCol->add(aTrack);
00967 TrackFinding(aTrack);
00968 p3DRoadC->add(road3D);
00969
00971
00972 vector<MucRecHit*> attachedHits = aTrack->GetHits();
00973 vector<MucRecHit*> expectedHits = aTrack->GetExpectedHits();
00974 vector<float> distanceHits = aTrack->getDistHits();
00975
00976 for(int i=0; i< expectedHits.size(); i++)
00977 {
00978 MucRecHit *ihit = expectedHits[i];
00979 //cout<<"expected Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<" "<<ihit->Strip()<<endl;
00980 }
00981
00982 vector<int> multiHit;
00983 for(int i=0; i< attachedHits.size(); i++)
00984 {
00985 MucRecHit *ihit = attachedHits[i];
00986 //cout<<"attached Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<" "<<ihit->Strip()<<" hitmode: "<<ihit->GetHitMode()<<endl;
00987 //calc multiplicity hits;
00988 int hitnum = 0;
00989 for(int k=0; k < attachedHits.size(); k++){
00990 MucRecHit *khit = attachedHits[k];
00991 if((ihit->Part()==khit->Part())&&(ihit->Seg()==khit->Seg())&&(ihit->Gap()==khit->Gap()))
00992 hitnum++;
00993 }
00994 multiHit.push_back(hitnum);
00995 //cout<<"multi hit: "<<hitnum<<" "<<multiHit[i]<<endl;
00996
00997 }
00998
00999 for(int i=0; i< expectedHits.size(); i++)
01000 { //calc distance between attached hits and expected hits
01001
01002 MucRecHit *ihit = expectedHits[i];
01003 //cout<<"attached Hits: "<<ihit->Part()<<" "<<ihit->Seg()<<" "<<ihit->Gap()<<" "<<ihit->Strip()<<endl;
01004
01005 int diff = -999;
01006
01007 for(int j=0; j< attachedHits.size(); j++){
01008 MucRecHit *jhit = attachedHits[j];
01009
01010 // if(attachedHits.size()!=distanceHits.size())cout<<"attached hits size no same as disthits!!!"<<endl;
01011
01012 if((jhit->Part()==ihit->Part())&&(jhit->Seg()==ihit->Seg())&&(jhit->Gap()==ihit->Gap())&&attachedHits.size()==distanceHits.size())
01013 { // same gap, cale distance
01014 diff = ihit->Strip() - jhit->Strip();
01015 //cout<<"diff = "<<diff<<endl;
01016
01017 if(m_NtOutput>=2){
01018
01019 m_part = ihit->Part(); m_seg = ihit->Seg(); m_gap = ihit->Gap();
01020 m_strip = jhit->Strip();
01021 m_diff = diff;
01022 m_dist = distanceHits[j];
01023 //cout<<"distance = "<<m_dist<<endl;
01024
01025 m_angle_cosmic = -999;
01026 m_angle_updown = -999;
01027 //m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01028 //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01029 //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01030 m_ngapwithhits = aTrack->numLayers();
01031 m_nhit = aTrack->numHits();
01032 m_maxhit = aTrack->maxHitsInLayer();
01033 m_multihit = multiHit[j];
01034 m_run = eventHeader->runNumber();
01035 m_event = eventHeader->eventNumber();
01036
01037 m_tuple->write();
01038 }
01039 }
01040
01041
01042 }
01043
01044
01045
01046 if(diff == -999) { // to calc effi of this strip
01047 if(m_NtOutput>=2){
01048 m_part = ihit->Part(); m_seg = ihit->Seg(); m_gap = ihit->Gap();
01049 m_strip = ihit->Strip();
01050 m_diff = diff;
01051 m_dist = -999;
01052 m_angle_updown = -999;
01053 m_angle_cosmic = -999;
01054 //m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01055 //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01056 //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01057 m_ngapwithhits = aTrack->numLayers();
01058 m_run = eventHeader->runNumber();
01059 m_event = eventHeader->eventNumber();
01060 //m_tuple->write();
01061 }
01062 }
01063 //if(diff != -999) cout<< "has hit in this layer"<<endl;
01064
01065 }
01067 /*
01068 m_part = -999;
01069 m_seg = -999;
01070 m_gap = -999;
01071 m_strip= -999;
01072 m_diff = -999;
01073
01074 m_angle_updown = -999;
01075 m_angle_cosmic = cosmicMom.angle(aTrack->getMucMomentum());
01076 if(m_angle_cosmic>1.57) m_angle_cosmic = 3.14159265 - m_angle_cosmic;
01077 m_run = eventHeader->runNumber();
01078 m_event = eventHeader->eventNumber();
01079
01080 double px,py,pz,p,theta,phi;
01081 px = (aTrack->getMucMomentum()).x();
01082 py = (aTrack->getMucMomentum()).y();
01083 pz = (aTrack->getMucMomentum()).z();
01084
01085 Hep3Vector rotate(-px,pz,py);
01086 theta = rotate.theta();
01087 phi = rotate.phi();
01088
01089
01090 m_px = px; m_py = py; m_pz = pz;
01091 m_theta = theta; m_phi = phi;
01092 m_theta_pipe = (aTrack->getMucMomentum()).theta();
01093 m_phi_pipe = (aTrack->getMucMomentum()).phi();
01094 //m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01095 //m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01097
01098 */
01099
01100 }
01101 else {
01102 //cout << "Delete a 3D Road ... " << endl;
01103 delete road3D;
01104 // don't keep it if it's not a good condidate
01105 }
01106 } // for ( int iRoad1 ...
01107 } // for ( int iRoad0 ..
01108
01109
01110 //for cosmic ray, to combine 2 track to 1
01111 RecMucTrack *aaTrack = 0;
01112 RecMucTrack *bbTrack = 0;
01113
01114 int hasMucUp = 0;
01115 int hasMucDown = 0;
01116 for(int iTrack = 0; iTrack < (int)aMucTrackCol->size(); iTrack++) {
01117 aaTrack = (*aMucTrackCol)[iTrack];
01118 if((aaTrack->getMucMomentum()).phi()<3.1415926&&(aaTrack->getMucMomentum()).phi()>0) hasMucUp++;
01119 else hasMucDown++;
01120
01121
01122 double px,py,pz,p,theta,phi;
01123 px = (aaTrack->getMucMomentum()).x();
01124 py = (aaTrack->getMucMomentum()).y();
01125 pz = (aaTrack->getMucMomentum()).z();
01126
01127 if(py<0)continue;
01128
01129 if(m_NtOutput>=1){
01130
01131 m_angle_updown = -999;
01132 m_angle_cosmic = cosmicMom.angle(aaTrack->getMucMomentum());
01133 if(m_angle_cosmic>1.57) m_angle_cosmic = 3.14159265 - m_angle_cosmic;
01134 m_run = eventHeader->runNumber();
01135 m_event = eventHeader->eventNumber();
01136
01137 Hep3Vector rotate(-px,pz,py);
01138 theta = rotate.theta();
01139 phi = rotate.phi();
01140
01141 m_px = px; m_py = py; m_pz = pz;
01142 m_theta = theta; m_phi = phi;
01143 m_theta_pipe = (aaTrack->getMucMomentum()).theta();
01144 m_phi_pipe = (aaTrack->getMucMomentum()).phi();
01145
01146 //if(fabs(m_phi_pipe*57.2958-90)>3)cout<<"px,y,z = "<<m_px<<" "<<m_py<<" "<<m_pz<<endl;
01147
01148 //calc proj point in y=0 plane
01149 Hep3Vector mucPos = aaTrack->getMucPos();
01150 double posx, posy, posz;
01151 posx = mucPos.x(); posy = mucPos.y(); posz = mucPos.z();
01152 //double projx, projz;
01153 m_projx = -999; m_projz = -999;
01154 if(py!=0){
01155 m_projx = posx - px/py*posy;
01156 m_projz = posz - pz/py*posy;
01157 }
01158 //cout<<"projection: "<<projx<<" "<<projz<<endl;
01159 }
01160
01161 }
01162 if(m_NtOutput>=1){
01163 m_mucUp = hasMucUp; m_mucDown = hasMucDown;
01164 m_tuple->write();
01165 }
01166
01167 int forCosmic = 0;
01168 if(aMucTrackCol->size()>=2 && forCosmic == 1){
01169 for(int iTrack = 0; iTrack < (int)aMucTrackCol->size(); iTrack++) {
01170 log << MSG::DEBUG << "iTrack " << iTrack << " / " <<(int)aMucTrackCol->size()<<endreq;
01171 aaTrack = (*aMucTrackCol)[iTrack];
01172 if(aaTrack->trackId()>=0) continue; // this track has been used
01173 aaTrack->setTrackId(iTrack);
01174 //cout<<"atrack set index "<<iTrack<<endl;
01175 for(int jTrack = iTrack+1; jTrack < (int)aMucTrackCol->size(); jTrack++){
01176 bbTrack = (*aMucTrackCol)[jTrack];
01177
01178 Hep3Vector mom_a = aaTrack->getMucMomentum();
01179 Hep3Vector mom_b = bbTrack->getMucMomentum();
01180
01181 //cout<<"angle is : "<<mom_a.angle(mom_b)<<endl;
01182 if(fabs(mom_a.angle(mom_b) - 3.1415926)<0.2)
01183 {
01184 bbTrack->setTrackId(iTrack);
01185 //cout<<"btrack set index "<<iTrack<<endl;
01186 //cout<<"find one cosmic ray"<<endl;
01187
01188 //cout<<"angle = "<<cosmicMom.angle(mom_a)<<" "<<cosmicMom.angle(mom_b)<<endl;
01189
01190 }
01191
01192 if(m_NtOutput>=1){
01193 m_angle_cosmic = cosmicMom.angle(mom_a);
01194 if(cosmicMom.angle(mom_a)>1.57) m_angle_cosmic = 3.14159265 - cosmicMom.angle(mom_a);
01195
01196 m_angle_updown = fabs(mom_a.angle(mom_b) - 3.1415926);
01197 m_px = -999; m_py = -999; m_pz = -999; m_theta = -999; m_phi = -999; m_theta_pipe = -999; m_phi_pipe = -999;
01198 m_px_mc = -999; m_py_mc = -999; m_pz_mc = -999; m_theta_mc = -999; m_phi_mc = -999;
01199 m_theta_mc_pipe = -999; m_phi_mc_pipe = -999;
01200
01201 //m_tuple->write();
01202 }
01203 }
01204
01205 }
01206
01207
01208 }
01209
01210 if( p3DRoadC->size() !=0 ) {
01211 log<<MSG::INFO << "In 3DRoad container : "
01212 << " Num of 3DRoad = " << p3DRoadC->size() <<endreq;
01213
01214 int print2DRoadInfo = 0;
01215 if (print2DRoadInfo == 1) {
01216 for ( int iRoad = 0; iRoad < (int)p3DRoadC->size(); iRoad++) {
01217 MucRec3DRoad *road = (*p3DRoadC)[iRoad];
01218 cout << endl << " " << iRoad << "th 3DRoad :" << endl
01219 << " Part = " << road->GetPart() << endl
01220 << " Seg = " << road->GetSeg() << endl
01221 << " NumGapsWithHits = " << road->GetNGapsWithHits() << endl
01222 << " TotalHits = " << road->GetTotalHits() << endl
01223 << " MaxHitsPerGap = " << road->GetMaxHitsPerGap() << endl
01224 << " LastGapDelta = " << road->GetLastGapDelta() << endl
01225 << " TotalHitsDelta = " << road->GetTotalHitsDelta() << endl
01226 << " DegreesOfFreedom = " << road->GetDegreesOfFreedom() << endl
01227 << " ReducedChiSquare = " << road->GetReducedChiSquare() << endl
01228 << " SlopeZX = " << road->GetSlopeZX() << endl
01229 << " InterceptZX = " << road->GetInterceptZX() << endl
01230 << " SlopeZY = " << road->GetSlopeZY() << endl
01231 << " InterceptZY = " << road->GetInterceptZY() << endl
01232 << " Hits Info : " << endl;
01233 //road->PrintHitsInfo();
01234 }
01235 }
01236
01237 m_NEventReco ++;
01238 }
01239
01240 //delete p3DRoadC
01241 for(int i = 0 ; i < p3DRoadC->size(); i++)
01242 delete (*p3DRoadC)[i];
01243 for(int i = 0 ; i < p2DRoad0C->size(); i++)
01244 delete (*p2DRoad0C)[i];
01245 for(int i = 0 ; i < p2DRoad1C->size(); i++)
01246 delete (*p2DRoad1C)[i];
01247
01248 p3DRoadC->clear();
01249 p2DRoad0C->clear();
01250 p2DRoad1C->clear();
01251 delete p3DRoadC;
01252 delete p2DRoad0C;
01253 delete p2DRoad1C;
01254 return StatusCode::SUCCESS;
01255 }
|
|
|
|
|
|
01259 {
01260 MsgStream log(msgSvc(), name());
01261 log << MSG::INFO << "in finalize()" << endreq;
01262
01263 //cout << m_NHitsLostTotal << " of " << m_NHitsTotal << " total hits" << endl;
01264 //for(int i = 0; i < 20; i++) cout << "lost " << i << " hits event " << m_NHitsLost[i] << endl;
01265 //for(int i = 0; i < 9; i++) cout << "lost on gap " << i << " is " << m_NHitsLostInGap[i] << endl;
01266
01267 return StatusCode::SUCCESS;
01268 }
|
|
|
|
|
|
00066 {
00067 MsgStream log(msgSvc(), name());
00068 log << MSG::INFO << "in initialize()" << endreq;
00069
00070 m_NHitsTotal = 0;
00071 m_NHitsLostTotal = 0;
00072 for(int i = 0; i < 20; i++) m_NHitsLost.push_back(0);
00073 for(int i = 0; i < 10; i++) m_NHitsLostInGap.push_back(0);
00074
00075 m_NEvent = 0;
00076 m_NEventWithHit = 0;
00077 m_NEventReco = 0;
00078
00079 IMucGeomSvc* mucGeomSvc;
00080 StatusCode sc = service("MucGeomSvc", mucGeomSvc);
00081 if (sc == StatusCode::SUCCESS) {
00082 mucGeomSvc->Dump();
00083 //cout<<"1st wire id:"<<mucGeomSvc->Wire(0)->Id()<<endl;
00084 //cout<<"2nd wire lyr id:"<<mucGeomSvc->Wire(1)->Lyr()->Id()<<endl;
00085 //cout<<"6860th wire lyr id:"<<mucGeomSvc->Wire(6859)->Lyr()->Id()<<endl;
00086 } else {
00087 return StatusCode::FAILURE;
00088 }
00089
00090 aMucRecHitContainer = 0;
00091
00092 if(m_NtOutput>=1){
00093 NTuplePtr nt1(ntupleSvc(), "FILE401/T");
00094
00095 if ( nt1 ) { m_tuple = nt1;}
00096 else {
00097 // m_tuple = ntupleSvc()->book ("MyTuples/1", CLID_RowWiseTuple, "MdcTrkRecon N-Tuple");
00098 m_tuple = ntupleSvc()->book ("FILE401/T", CLID_RowWiseTuple, "MucTrkRecon N-Tuple");
00099 if ( m_tuple ) {
00100 sc = m_tuple->addItem ("part", m_part);
00101 sc = m_tuple->addItem ("seg", m_seg);
00102 sc = m_tuple->addItem ("gap", m_gap);
00103 sc = m_tuple->addItem ("strip", m_strip);
00104 sc = m_tuple->addItem ("diff", m_diff);
00105 sc = m_tuple->addItem ("dist", m_dist);
00106 sc = m_tuple->addItem ("run", m_run);
00107 sc = m_tuple->addItem ("event", m_event);
00108 sc = m_tuple->addItem ("ngap", m_ngapwithhits);
00109 sc = m_tuple->addItem ("nhit", m_nhit);
00110 sc = m_tuple->addItem ("maxhit", m_maxhit);
00111 sc = m_tuple->addItem ("multihit",m_multihit);
00112 sc = m_tuple->addItem ("angleCosmic",m_angle_cosmic);
00113 sc = m_tuple->addItem ("angleUpdown",m_angle_updown);
00114 sc = m_tuple->addItem ("px",m_px);
00115 sc = m_tuple->addItem ("py",m_py);
00116 sc = m_tuple->addItem ("pz",m_pz);
00117 sc = m_tuple->addItem ("theta",m_theta);
00118 sc = m_tuple->addItem ("phi",m_phi);
00119 sc = m_tuple->addItem ("theta_pipe",m_theta_pipe);
00120 sc = m_tuple->addItem ("phi_pipe",m_phi_pipe);
00121 sc = m_tuple->addItem ("pxmc",m_px_mc);
00122 sc = m_tuple->addItem ("pymc",m_py_mc);
00123 sc = m_tuple->addItem ("pzmc",m_pz_mc);
00124 sc = m_tuple->addItem ("thetamc",m_theta_mc);
00125 sc = m_tuple->addItem ("phimc",m_phi_mc);
00126 sc = m_tuple->addItem ("thetamc_pipe",m_theta_mc_pipe);
00127 sc = m_tuple->addItem ("phimc_pipe",m_phi_mc_pipe);
00128 sc = m_tuple->addItem ("emcUp",m_emcUp);
00129 sc = m_tuple->addItem ("emcDown",m_emcDown);
00130 sc = m_tuple->addItem ("mucUp",m_mucUp);
00131 sc = m_tuple->addItem ("mucDown",m_mucDown);
00132 sc = m_tuple->addItem ("projx",m_projx);
00133 sc = m_tuple->addItem ("projz",m_projz);
00134
00135 }
00136 else { // did not manage to book the N tuple....
00137 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
00138 //return StatusCode::FAILURE;
00139 }
00140 }
00141 }
00142
00143 return StatusCode::SUCCESS;
00144 }
|
|
|
|
|
|
01272 {
01273 MsgStream log(msgSvc(), name());
01274
01275 Hep3Vector currentPos = aTrack->GetCurrentPos();
01276 Hep3Vector currentDir = aTrack->GetCurrentDir();
01277 // if(currentPos.mag() < kMinor) {
01278 // log << MSG::WARNING << "No MUC intersection in track " << endreq;
01279 // continue;
01280 // }
01281
01282 int firstHitFound[2] = { 0, 0}; // Has the fist position in this orient determined? if so, could CorrectDirection()
01283 int firstHitGap[2] = {-1, -1}; // When the fist position in this orient determined, the gap it is on
01284 for(int partSeq = 0; partSeq < (int)MucID::getPartNum(); partSeq++) {
01285 int iPart = kPartSeq[partSeq];
01286 for(int iGap = 0; iGap < (int)MucID::getGapNum(iPart); iGap++) {
01287 int seg = -1;
01288 int orient = MucGeoGeneral::Instance()->GetGap(iPart, 0, iGap)->Orient();;
01289
01290 float xInsct, yInsct, zInsct;
01291 aTrack->Project(iPart, iGap, xInsct, yInsct, zInsct, seg);
01292 if(m_MsOutput) cout<<"part "<<iPart<<" gap "<<iGap<<" x "<<xInsct<<" y "<<yInsct<<" z "<<zInsct<<" seg "<<seg<<endl;
01293
01294 if(seg == -1) {
01295 //log << MSG::DEBUG << "no intersection with part " << iPart
01296 // << " gap " << iGap << " in this track !" << endl;
01297 continue;
01298 }
01299
01300 aTrack->SetCurrentInsct(xInsct, yInsct, zInsct);
01301
01302 for(int iDeltaSeg = 0; iDeltaSeg < kNSegSearch; iDeltaSeg++) {
01303 int iSeg = seg + kDeltaSeg[iDeltaSeg]; // also find on neighbor seg;
01304 if(iSeg < 0) iSeg += MucID::getSegNum(iPart);
01305 if(iSeg > (int)MucID::getSegNum(iPart) - 1) iSeg -= MucID::getSegNum(iPart);
01306
01307 float window = kWindowSize[iPart][iGap];
01308
01309 for (int iHit = 0; iHit < aMucRecHitContainer->GetGapHitCount(iPart, iSeg, iGap); iHit++) {
01310 log << MSG::DEBUG << "iSeg " << iSeg << " iHit " << iHit << endreq;
01311 MucRecHit* pHit = aMucRecHitContainer->GetHit(iPart, iSeg, iGap, iHit);
01312 //cout<< "strip " << pHit->Strip() << " center " << pHit->GetCenterPos() << endl;
01313
01314 if (!pHit) {
01315 log << MSG::WARNING << "MucRecTrkExt: null pointer to pHit" << endreq;
01316 continue;
01317 }
01318 else {
01319
01320 //this algo use with ext and this hit has been used before;
01321 if(pHit->GetHitMode() != -1 && pHit->GetHitMode() != 3) continue;
01322
01323 // Get the displacement of hit pHit to intersection
01324 float dX = aTrack->GetHitDistance(pHit);
01325 log << MSG::DEBUG << "distance = " << setw(8) << dX << " size " << setw(4) << window << endreq;
01326
01327 //cout<<"dX= "<<dX<<" window="<<window<<endl;
01328
01329 if (dX < window) {
01330 // Attach the hit if it exists
01331 //cout << "meet window " << pHit->GetSoftID() << endl;
01332 //****************if this if emc track, we abondon used hit in mdc*******************
01333 //if(m_emcrec == 1 )
01334 if(aTrack->GetRecMode() == 0) {
01335 pHit->SetHitMode(1); //mdc ext
01336 aTrack->AttachHit(pHit);
01337 //cout<<"in MucRecTrkExt: trackmode==0 so mdc ext "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01338 }
01339 if(aTrack->GetRecMode() == 1) {
01340 //cout<<"in MucRecTrkExt: HitMode = "<<pHit->GetHitMode()<<" part: "<<iPart<<" "<<iSeg<<" "<<iGap<<" "<<iHit<<endl;
01341 if(pHit->GetHitMode()!=1) {
01342 aTrack->AttachHit(pHit); //this hit has not been used by mdc ext
01343 pHit->SetHitMode(2); //emc ext
01344 }
01345 }
01346 if(aTrack->GetRecMode() == 2) {
01347 aTrack->AttachHit(pHit); //this hit has not been used by mdc ext
01348 //pHit->SetHitMode(2); //emc ext
01349 }
01350 if(aTrack->GetRecMode() == 3) {
01351 if(pHit->GetHitMode()!=1) {
01352 aTrack->AttachHit(pHit); //this hit has not been used by mdc ext
01353 pHit->SetHitMode(3); //emc ext
01354 }
01355 }
01356
01357 if (firstHitGap[orient] == -1) firstHitGap[orient] = iGap;
01358 firstHitFound[orient] = 1;
01359 //cout << "You could correct directon in orient " << orient << endl;
01360
01361 //cout<< " part " << iPart << " seg " << iSeg << " gap " << iGap
01362 // << " strip " << setw(2) << pHit->Strip() << " attatched" << endl;
01363
01364 log << MSG::DEBUG << " part " << iPart << " seg " << iSeg << " gap " << iGap
01365 << " strip " << setw(2) << pHit->Strip() << " attatched" << endreq;
01366 log << MSG::DEBUG << "current total hits " << aTrack->GetTotalHits() << endreq;
01367 }
01368 else {
01369 m_NHitsLostInGap[iGap]++;
01370 //log << MSG::DEBUG cout << " part " << iPart << " seg " << iSeg << " gap " << iGap
01371 // << " strip " << setw(2) << pHit->GetSoftID().GetStrip()
01372 // << " not attached !" << endreq;
01373 }
01374 }
01375 }
01376 aTrack->CalculateInsct(iPart, iSeg, iGap);
01377 }
01378
01379 if(m_onlyseedfit == 0) {//
01380 // When correct dir in the orient is permitted and this gap is not the gap first hit locates.
01381 if (firstHitFound[orient] && firstHitGap[orient] != iGap) aTrack->CorrectDir();
01382 aTrack->CorrectPos();
01383 //cout << "Current Intersection " << aTrack->GetCurrentInsct() << endl;
01384 //cout << "Current Direction " << aTrack->GetCurrentDir() << endl;
01385 aTrack->AttachInsct(aTrack->GetCurrentInsct());
01386 }// else; not correct pos & dir.
01387
01388 }
01389 }
01390 aTrack->LineFit(1);
01391 aTrack->ComputeTrackInfo(1);
01392 log << MSG::INFO << "Fit track done! trackIndex: " << aTrack->trackId() << ", mucId: " << aTrack->id() << ", RecMode: " << aTrack->GetRecMode() << endreq;
01393 //cout<<" in TrackFinding: depth= "<<aTrack->GetDepth()<<endl;
01394 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1