#include <DQA_MUC.h>
|
||||||||||||
|
00036 :
00037 Algorithm(name, pSvcLocator) {
00038
00039 //Declare the properties
00040 declareProperty("EffWindow", m_effWindow = 4);
00041 }
|
|
||||||||||||
|
|
|
|
|
|
|
00227 {
00228
00229 MsgStream log(msgSvc(), name());
00230 log << MSG::INFO << "in execute()" << endreq;
00231
00232 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
00233 m_run = eventHeader->runNumber();
00234 m_event = eventHeader->eventNumber();
00235 log << MSG::DEBUG <<"Run " << m_run << "\tEvent " << m_event <<endreq;
00236
00237 SmartDataPtr<DQAEvent::DQAEvent> dqaevt(eventSvc(), "/Event/DQATag");
00238 if ( dqaevt ) {
00239 log << MSG::INFO << "success get DQAEvent" << endreq;
00240 } else {
00241 log << MSG::ERROR << "Error accessing DQAEvent" << endreq;
00242 return StatusCode::FAILURE;
00243 }
00244 log << MSG::DEBUG << "DQA event tag = " << dqaevt->EventTag() << endreq;
00245
00246 int part, segment, layer, strip;
00247 char name[100];
00248 bool isDimu = false;
00249 if ( dqaevt->Dimu() ) isDimu = true;
00250 log << MSG::INFO << "DQADimuTag:\t" << dqaevt->Dimu() << endreq;
00251
00252 //---> Retrieve MUC digi
00253 log << MSG::INFO << "Retrieve digis" << endreq;
00254 SmartDataPtr<MucDigiCol> mucDigiCol(eventSvc(),"/Event/Digi/MucDigiCol");
00255 if(!mucDigiCol) {
00256 log << MSG::FATAL << "Could not find MUC digi" << endreq;
00257 return( StatusCode::FAILURE);
00258 }
00259
00260 Identifier mucId;
00261 MucDigiCol::iterator digiIter = mucDigiCol->begin();
00262 int eventDigi = 0;
00263 for ( int digiId =0; digiIter != mucDigiCol->end(); digiIter++, digiId++ )
00264 {
00265 mucId = (*digiIter)->identify();
00266 part = MucID::barrel_ec(mucId); segment = MucID::segment(mucId);
00267 layer = MucID::layer(mucId); strip = MucID::channel(mucId);
00268 eventDigi ++;
00269
00270 // Add digi
00271 MucMark *aMark = new MucMark( part, segment, layer, strip );
00272 m_digiCol.push_back( aMark );
00273 m_segDigiCol[part][segment].push_back( aMark );
00274 m_recordAll[part][segment][layer][strip][0] ++;
00275 if( isDimu ) m_recordDimu[part][segment][layer][strip][0] ++;
00276 }
00277 log << MSG::INFO << "Total digis in this event: " << eventDigi << endreq;
00278 if( eventDigi > 500) cout <<"Run:\t"<< m_run <<"\tEvent:\t"<< m_event << "\tdigits inflation:\t" << eventDigi << endreq;
00279
00280 // Search cluster in digis
00281 m_clusterCol = (*m_ptrMucMark).CreateClusterCol(1, m_digiCol );
00282 log << MSG::INFO << "Total clusters in this event: " << m_clusterCol.size() << endreq;
00283 //<--- End retrieve digis
00284
00285 //---> Retrieve rec tracks
00286 SmartDataPtr<EvtRecEvent> evtRecEvent(eventSvc(), EventModel::EvtRec::EvtRecEvent);
00287 SmartDataPtr<RecMucTrackCol> mucTrackCol(eventSvc(),"/Event/Recon/RecMucTrackCol");
00288 if (!mucTrackCol) {
00289 log << MSG::FATAL << "Could not find RecMucTrackCol" << endreq;
00290 return( StatusCode::FAILURE);
00291 }
00292
00293 RecMucTrackCol *aRecMucTrackCol = mucTrackCol;
00294 if (aRecMucTrackCol->size() < 1) {
00295 log << MSG::INFO << "No MUC tracks in this event" << endreq;
00296 return StatusCode::SUCCESS;
00297 }
00298 log << MSG::INFO << "Total tracks of this event: " << aRecMucTrackCol->size() << endreq;
00299
00300 int trackHitNum, expectedHitNum, segNum, lastLayerBR, lastLayerEC;
00301 int layerPassNum[3], passMax[TRACK_SEG_MAX][2];
00302 bool firedLay[TRACK_SEG_MAX][LAYER_MAX];
00303 double costheta, phi;
00304 TH1* htmp(0);
00305 trackHitNum = expectedHitNum = segNum = lastLayerBR = lastLayerEC = 0;
00306 layerPassNum[0] = layerPassNum[1] = layerPassNum[2] = 0;
00307 for( int segi=0; segi<TRACK_SEG_MAX; segi++ ) {
00308 passMax[segi][0] = passMax[segi][1] = 0;
00309 for( int layi=0; layi<LAYER_MAX; layi++ ) firedLay[segi][layi] = 0;
00310 }
00311
00312 vector<MucRecHit*> mucRawHitCol;
00313 vector<MucRecHit*> mucExpHitCol;
00314 RecMucTrackCol::iterator trackIter = mucTrackCol->begin();
00315 for (int trackId = 0; trackIter != mucTrackCol->end(); trackIter++, trackId++)
00316 {
00317 trackHitNum = (*trackIter)->numHits();
00318
00319 if( trackHitNum == 0 ) {
00320 log << MSG::INFO << "Track " << trackId << " no hits" << endreq;
00321 continue;
00322 }
00323
00324 lastLayerBR = (*trackIter)->brLastLayer();
00325 lastLayerEC = (*trackIter)->ecLastLayer();
00326 // First fit position in MUC
00327 CLHEP::Hep3Vector a3Vector((*trackIter)->xPos(),(*trackIter)->yPos(),(*trackIter)->zPos());
00328 costheta = a3Vector.cosTheta();
00329 phi = a3Vector.phi();
00330 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_All", htmp).isSuccess()) {
00331 htmp->Fill( costheta );
00332 } else {
00333 log << MSG::ERROR << "Fail to retrieve Costheta_All" << endreq;
00334 }
00335 if(m_thsvc->getHist("/DQAHist/MUC/Phi_All", htmp).isSuccess()) {
00336 htmp->Fill( phi );
00337 } else {
00338 log << MSG::ERROR << "Fail to retrieve Phi_All" << endreq;
00339 }
00340
00341 if( isDimu )
00342 {
00343 if(m_thsvc->getHist("/DQAHist/MUC/Costheta_Dimu", htmp).isSuccess()) {
00344 htmp->Fill( costheta );
00345 } else {
00346 log << MSG::ERROR << "Fail to retrieve Costheta_Dimu" << endreq;
00347 }
00348 if(m_thsvc->getHist("/DQAHist/MUC/Phi_Dimu", htmp).isSuccess()) {
00349 htmp->Fill( phi );
00350 } else {
00351 log << MSG::ERROR << "Fail to retrieve Phi_Dimu" << endreq;
00352 }
00353 }
00354 log << MSG::INFO << "Fill costheta and phi:\t" << costheta << "\t" << phi << endreq;
00355
00356 MucRecHit* pMucRawHit;
00357 MucRecHit* pMucExpHit;
00358
00359
00360 // Expected hits in this rec track
00361 mucExpHitCol = (*trackIter)->GetExpectedHits();
00362 //mucExpHitCol = (*trackIter)->getExpHits();
00363 expectedHitNum += mucExpHitCol.size();
00364 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
00365 {
00366 pMucRawHit = mucExpHitCol[ hitId ];
00367 part = pMucRawHit->Part(); segment = pMucRawHit->Seg();
00368 layer = pMucRawHit->Gap(); strip = pMucRawHit->Strip();
00369
00370 MucMark* currentMark = new MucMark( part, segment, layer, strip );
00371 m_expHitCol.push_back( currentMark );
00372
00373 // Judge efficiency hit
00374 int isInPos = -1;
00375 bool isInEffWindow = false;
00376 isInPos = currentMark->IsInCol( m_segDigiCol[part][segment] );
00377
00378 // Avoid bias in outer layers caused by low momentum tracks
00379 if( part == BRID && (layer-lastLayerBR>1) ) continue;
00380 if( part != BRID && (layer-lastLayerEC>1) ) continue;
00381
00382 // Avoid bias in both sides of the innermost layer of Barrel
00383 if( part==BRID && layer==0 && (strip<2 || strip>45) )
00384 {
00385 if( isInPos != -1) // expHit is fired
00386 {
00387 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
00388 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
00389 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
00390
00391 if(isDimu) {
00392 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
00393 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
00394 }
00395 }
00396 else {
00397 m_recordAll[part][segment][layer][strip][1] ++;
00398 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
00399 }
00400 continue; // Judge next hit
00401 }
00402 // Eff calibration
00403 if( isInPos != -1 ) // expHit is fired
00404 {
00405 m_recordAll[part][segment][layer][strip][2] ++; // Efficiency hit number
00406 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
00407 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
00408
00409 if(isDimu) {
00410 m_recordDimu[part][segment][layer][strip][2] ++; // Efficiency hit number
00411 m_recordDimu[part][segment][layer][strip][1] ++; // Rec track number
00412 }
00413
00414 continue; // Judge next hit
00415 }
00416 else for(int tempStrip=0, hiti=-m_effWindow; hiti<=m_effWindow; hiti++ )
00417 {
00418 if( hiti == 0 ) continue;
00419 tempStrip = strip + hiti;
00420 if( tempStrip < 0 || tempStrip > m_ptrIdTr->GetStripMax(part,segment,layer) ) continue;
00421
00422 isInPos = m_ptrMucMark->IsInCol( part, segment, layer, tempStrip, m_segDigiCol[part][segment] );
00423 if( isInPos != -1 )
00424 {
00425 m_recordAll[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
00426 m_recordAll[part][segment][layer][tempStrip][1] ++; // Rec track number
00427 m_effHitCol.push_back( m_segDigiCol[part][segment][isInPos] );
00428
00429 if(isDimu) {
00430 m_recordDimu[part][segment][layer][tempStrip][2] ++; // Efficiency hit number
00431 m_recordDimu[part][segment][layer][tempStrip][1] ++; // Rec track number
00432 }
00433
00434 isInEffWindow = true;
00435 }
00436 } // End else
00437
00438 if( isInEffWindow ) { continue; } // Judge next hit
00439 else { // A hit should be fired but not fired and not in the EffWindow
00440 m_recordAll[part][segment][layer][strip][1] ++; // Rec track number
00441 if(isDimu) m_recordDimu[part][segment][layer][strip][1] ++;
00442 }
00443
00444 } // End expected hits
00445
00446 // Fill residual, and for the other way of eff calculation
00447 log << MSG::INFO << "Fill residual" << endreq;
00448 vector<float> m_lineResCol = (*trackIter)->getDistHits();
00449 // Digis belong to this rec track
00450 mucRawHitCol = (*trackIter)->GetHits(); // Get hit collection of a track
00451 //mucRawHitCol = (*trackIter)->getVecHits(); // Get hit collection of a track
00452
00453 if( trackHitNum > 4 && m_lineResCol[0] != -99) // track is good for res
00454 {
00455 // Fill res histograms
00456 bool firedFlag[PART_MAX][LAYER_MAX][2];
00457 for(int iprt=0; iprt<PART_MAX; iprt++)
00458 for(int jlay=0; jlay<LAYER_MAX; jlay++)
00459 firedFlag[iprt][jlay][0] = firedFlag[iprt][jlay][1] = false;
00460
00461 for(unsigned int hitId = 0; hitId < mucExpHitCol.size(); hitId++)
00462 {
00463 pMucExpHit = mucExpHitCol[ hitId ];
00464 part = pMucExpHit->Part(); segment = pMucExpHit->Seg(); layer = pMucExpHit->Gap();
00465 firedFlag[part][layer][0] = true;
00466 }
00467
00468 for(unsigned int hitId = 0; hitId < mucRawHitCol.size(); hitId++)
00469 {
00470 pMucRawHit = mucRawHitCol[ hitId ];
00471 part = pMucRawHit->Part(); segment = pMucRawHit->Seg(); layer = pMucRawHit->Gap();
00472
00473 // if exp is true and fired is true, and not filled yet
00474 if( firedFlag[part][layer][0] == true && firedFlag[part][layer][1] == false )
00475 {
00476 if( part == BRID )
00477 {
00478 sprintf( name, "/DQAHist/MUC/BrResDist_All_L%d", layer );
00479 if(m_thsvc->getHist(name, htmp).isSuccess()) {
00480 htmp->Fill( m_lineResCol[hitId] );
00481 } else {
00482 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00483 }
00484 //m_hBrResDist[layer][0]->Fill( m_lineResCol[hitId] );
00485
00486 if( isDimu )
00487 {
00488 sprintf( name, "/DQAHist/MUC/BrResDist_Dimu_L%d", layer );
00489 if(m_thsvc->getHist(name, htmp).isSuccess()) {
00490 htmp->Fill( m_lineResCol[hitId] );
00491 } else {
00492 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00493 }
00494 //m_hBrResDist[layer][1]->Fill( m_lineResCol[hitId] );
00495 }
00496 } else {
00497 sprintf( name, "/DQAHist/MUC/EcResDist_All_L%d", layer );
00498 if(m_thsvc->getHist(name, htmp).isSuccess()) {
00499 htmp->Fill( m_lineResCol[hitId] );
00500 } else {
00501 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00502 }
00503 //m_hEcResDist[layer][0]->Fill( m_lineResCol[hitId] );
00504
00505 if( isDimu )
00506 {
00507 sprintf( name, "/DQAHist/MUC/EcResDist_Dimu_L%d", layer );
00508 if(m_thsvc->getHist(name, htmp).isSuccess()) {
00509 htmp->Fill( m_lineResCol[hitId] );
00510 } else {
00511 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00512 }
00513 //m_hEcResDist[layer][1]->Fill( m_lineResCol[hitId] );
00514 }
00515 }
00516 }
00517
00518 firedFlag[part][layer][1] = true;
00519 }
00520 } // End fill residual, if track is good for res
00521
00522 mucRawHitCol.clear();
00523 mucExpHitCol.clear();
00524
00525 } // End read all tracks
00526 //<--- End retrieve rec tracks
00527
00528 //---> Searching inc/noise hits
00529 log << MSG::INFO << "Searching inc/noise hits" << endreq;
00530 bool isNosHit;
00531 bool hasEffHit;
00532 for( unsigned int i=0; i < m_digiCol.size(); i++ )
00533 {
00534 isNosHit = true;
00535
00536 if( m_digiCol[i]->IsInCol( m_effHitCol )!=-1 ) continue; // digi in effHitCol
00537 else
00538 {
00539 for( unsigned int j=0; j < m_clusterCol.size(); j++ )
00540 {
00541 hasEffHit = false;
00542 for( unsigned int k=0; k<m_clusterCol[j].size(); k++)
00543 {
00544 if( m_clusterCol[j][k]->IsInCol(m_effHitCol) != -1)
00545 { // Clusters have efficiency hit
00546 hasEffHit = true;
00547 break; // Out a cluster
00548 }
00549 }
00550
00551 if ( hasEffHit && (m_digiCol[i]->IsInCol( m_clusterCol[j] ) != -1) ) {
00552 isNosHit = false;
00553 break; // Out cluster collection
00554 }
00555 } // End cluster col
00556
00557 if( isNosHit ) {
00558 part = (*m_digiCol[i]).Part(); segment = (*m_digiCol[i]).Segment();
00559 layer = (*m_digiCol[i]).Layer(); strip = (*m_digiCol[i]).Strip();
00560
00561 //log << MSG::INFO << "Noise hit:\t" << part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip << endreq;
00562 m_recordAll[part][segment][layer][strip][3] ++;
00563 if( isDimu ) m_recordDimu[part][segment][layer][strip][3] ++;
00564 m_nosHitCol.push_back(m_digiCol[i]);
00565 }
00566 }// End else
00567 } // End digi collection
00568 //<--- End searching inc/noise hits
00569
00570 DQA_MUC::FillHistograms( isDimu );
00571
00572 // Reset mark collection
00573 for( unsigned int i=0; i<m_digiCol.size(); i++ ) {
00574 if( m_digiCol[i] != NULL ) delete m_digiCol[i];
00575 }
00576
00577 for( unsigned int i=0; i<m_expHitCol.size(); i++ ) {
00578 if( m_expHitCol[i] != NULL ) delete m_expHitCol[i];
00579 }
00580
00581 if( m_effHitCol.size() != 0 ) m_effHitCol.clear();
00582 if( m_expHitCol.size() != 0 ) m_expHitCol.clear();
00583 if( m_nosHitCol.size() != 0 ) m_nosHitCol.clear();
00584 if( m_digiCol.size() != 0 ) m_digiCol.clear();
00585 if( m_clusterCol.size() != 0 ) m_clusterCol.clear();
00586
00587 for( int i=0; i<PART_MAX; i++ ) {
00588 for( int j=0; j<SEGMENT_MAX; j++ ) {
00589 if( m_segDigiCol[i][j].size() != 0 ) m_segDigiCol[i][j].clear();
00590 }
00591 }
00592
00593 return StatusCode::SUCCESS;
00594 }
|
|
|
|
|
|
00834 {
00835
00836 MsgStream log(msgSvc(), name());
00837 log << MSG::INFO << "Filling histograms" << endmsg;
00838 char name[100];
00839 int part, segment, layer, strip, boxId, strId;
00840 part = segment = layer = strip = boxId = strId = 0;
00841 TH1 *hEff[2][2][3];
00842
00843 for(int i=0; i<TAGN; i++) {
00844 for(int j=0; j<2; j++) {
00845 for(int k=0; k<LVLN; k++)
00846 {
00847 sprintf( name, "/DQAHist/MUC/%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
00848 if(!m_thsvc->getHist(name, hEff[i][j][k]).isSuccess())
00849 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00850 }
00851 }
00852 }
00853
00854 // Numberator of eff
00855 for(int i=0; i<m_effHitCol.size(); i++)
00856 {
00857 part = m_effHitCol[i]->Part();
00858 segment = m_effHitCol[i]->Segment();
00859 layer = m_effHitCol[i]->Layer();
00860 strip = m_effHitCol[i]->Strip();
00861 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
00862 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
00863 hEff[0][0][0]->Fill(layer);
00864 hEff[0][0][1]->Fill(boxId);
00865 hEff[0][0][2]->Fill(strId);
00866
00867 if( isDimu ) {
00868 hEff[1][0][0]->Fill(layer);
00869 hEff[1][0][1]->Fill(boxId);
00870 hEff[1][0][2]->Fill(strId);
00871 }
00872 }
00873
00874 // Denominator of eff
00875 for(int i=0; i<m_expHitCol.size(); i++)
00876 {
00877 part = m_expHitCol[i]->Part();
00878 segment = m_expHitCol[i]->Segment();
00879 layer = m_expHitCol[i]->Layer();
00880 strip = m_expHitCol[i]->Strip();
00881 //cout<<"ExpHit:\t"<<part<<"\t"<<segment<<"\t"<<layer<<"\t"<<strip<<endl;
00882 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
00883 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
00884 hEff[0][1][0]->Fill(layer);
00885 hEff[0][1][1]->Fill(boxId);
00886 hEff[0][1][2]->Fill(strId);
00887
00888 if( isDimu ) {
00889 hEff[1][1][0]->Fill(layer);
00890 hEff[1][1][1]->Fill(boxId);
00891 hEff[1][1][2]->Fill(strId);
00892 }
00893
00894 }
00895
00896 TH1 *hNosRatio[2][2];
00897 for(int i=0; i<TAGN; i++) {
00898 for(int j=0; j<2; j++) {
00899 sprintf( name, "/DQAHist/MUC/%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
00900 if(!m_thsvc->getHist(name, hNosRatio[i][j]).isSuccess())
00901 log << MSG::ERROR << "Fail to retrieve " << name << endreq;
00902 }
00903 }
00904
00905 // Numberator of box noise ratio
00906 for(int i=0; i<m_nosHitCol.size(); i++)
00907 {
00908 part = m_nosHitCol[i]->Part();
00909 segment = m_nosHitCol[i]->Segment();
00910 layer = m_nosHitCol[i]->Layer();
00911 strip = m_nosHitCol[i]->Strip();
00912 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
00913 strId = m_ptrIdTr->GetStripId(part, segment, layer, strip);
00914 hNosRatio[0][0]->Fill(boxId);
00915
00916 if( isDimu ) hNosRatio[1][0]->Fill(boxId);
00917 }
00918
00919 // Denominator of box noise ratio
00920 for(int i=0; i<m_digiCol.size(); i++)
00921 {
00922 part = m_digiCol[i]->Part();
00923 segment = m_digiCol[i]->Segment();
00924 layer = m_digiCol[i]->Layer();
00925 strip = m_digiCol[i]->Strip();
00926 boxId = m_ptrIdTr->GetBoxId(part, segment, layer);
00927 hNosRatio[0][1]->Fill(boxId);
00928
00929 if( isDimu ) hNosRatio[1][1]->Fill(boxId);
00930 }
00931
00932 return StatusCode::SUCCESS;
00933 }
|
|
|
|
|
|
00597 {
00598
00599 MsgStream log(msgSvc(), name());
00600 log << MSG::INFO << "in finalize()" << endmsg;
00601
00602 /*
00603 TH1 *htmp(0);
00604
00605 // Fit spacial resolution
00606 log << MSG::INFO << "Spacial resolution fitting" << endreq;
00607 double resSigma, resSigmaErr;
00608 resSigma = resSigmaErr = 0.;
00609
00610 for( int i=0; i<B_LAY_NUM; i++ )
00611 {
00612 m_hBrResDist[i][0]->Fit("gaus");
00613 resSigma = m_hBrResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
00614 resSigmaErr = m_hBrResDist[i][0]->GetFunction("gaus")->GetParError(2);
00615 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_All", htmp).isSuccess()) {
00616 htmp->Fill( i, resSigma );
00617 htmp->SetBinError( i+1, resSigmaErr );
00618 } else {
00619 log << MSG::ERROR << "Fail to retrieve BrRes_All" << endreq;
00620 }
00621 //log << MSG::INFO << "Barrle layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
00622
00623 m_hBrResDist[i][1]->Fit("gaus");
00624 resSigma = m_hBrResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
00625 resSigmaErr = m_hBrResDist[i][1]->GetFunction("gaus")->GetParError(2);
00626 if(m_thsvc->getHist("/DQAHist/MUC/BrRes_Dimu", htmp).isSuccess()) {
00627 htmp->Fill( i, resSigma );
00628 htmp->SetBinError( i+1, resSigmaErr );
00629 } else {
00630 log << MSG::ERROR << "Fail to retrieve BrRes_Dimu" << endreq;
00631 }
00632
00633 }
00634
00635 for( int i=0; i<E_LAY_NUM; i++ )
00636 {
00637 m_hEcResDist[i][0]->Fit("gaus");
00638 resSigma = m_hEcResDist[i][0]->GetFunction("gaus")->GetParameter("Sigma");
00639 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
00640 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_All", htmp).isSuccess()) {
00641 htmp->Fill( i, resSigma );
00642 htmp->SetBinError( i+1, resSigmaErr );
00643 } else {
00644 log << MSG::ERROR << "Fail to retrieve EcRes_All" << endreq;
00645 }
00646 //log << MSG::INFO << "Endcap layer:\t" << i << "\t" << resSigma << "\t" << resSigmaErr << endreq;
00647
00648 m_hEcResDist[i][1]->Fit("gaus");
00649 resSigma = m_hEcResDist[i][1]->GetFunction("gaus")->GetParameter("Sigma");
00650 resSigmaErr = m_hEcResDist[i][0]->GetFunction("gaus")->GetParError(2);
00651 if(m_thsvc->getHist("/DQAHist/MUC/EcRes_Dimu", htmp).isSuccess()) {
00652 htmp->Fill( i, resSigma );
00653 htmp->SetBinError( i+1, resSigmaErr );
00654 } else {
00655 log << MSG::ERROR << "Fail to retrieve EcRes_Dimu" << endreq;
00656 }
00657 }
00658
00659
00660 // Efficiency and noise ratio
00661 int part, segment, layer, stripMax;
00662 part = segment = layer = stripMax = 0;
00663 double digi[TAGN][3], effHit[TAGN][3], trackNum[TAGN][3], nosHit[TAGN][3];
00664 double eff[TAGN], effErr[TAGN], nos[TAGN], nosErr[TAGN];
00665 for(int i=0; i<TAGN; i++) {
00666 for( int j=0; j<3; j++) {
00667 digi[i][j] = effHit[i][j] = trackNum[i][j] = nosHit[i][j] = 0.;
00668 }
00669
00670 eff[i] = effErr[i] = nos[i] = nosErr[i] = 0.;
00671 }
00672
00673 // Layer efficiency
00674 for( int i=0; i<LAYER_MAX; i++ )
00675 {
00676 effHit[0][0] = trackNum[0][0] = 0;
00677 effHit[1][0] = trackNum[1][0] = 0;
00678
00679 for( int j=0; j<PART_MAX; j++ ) {
00680 for( int k=0; k<SEGMENT_MAX; k++) {
00681 stripMax = m_ptrIdTr->GetStripMax( j, k, i );
00682 for( int n=0; n<stripMax; n++ ) {
00683 trackNum[0][0] += m_recordAll[j][k][i][n][1];
00684 effHit[0][0] += m_recordAll[j][k][i][n][2];
00685
00686 trackNum[1][0] += m_recordDimu[j][k][i][n][1];
00687 effHit[1][0] += m_recordDimu[j][k][i][n][2];
00688 }
00689 }
00690 }
00691
00692 if( trackNum[0][0] == 0 ) {
00693 eff[0] = effErr[0] = 0.0;
00694 } else {
00695 eff[0] = ( (double)effHit[0][0] )/trackNum[0][0];
00696 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][0] );
00697 }
00698 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_All", htmp).isSuccess()) {
00699 htmp->Fill( i, eff[0] );
00700 htmp->SetBinError( i+1, effErr[0] );
00701 } else {
00702 log << MSG::ERROR << "Fail to retrieve LayerEff_All" << endreq;
00703 }
00704 log << MSG::INFO << "LayerEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
00705
00706 if( trackNum[1][0] == 0 ) {
00707 eff[1] = effErr[1] = 0.0;
00708 } else {
00709 eff[1] = ( (double)effHit[1][0] )/trackNum[1][0];
00710 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][0] );
00711 }
00712 if(m_thsvc->getHist("/DQAHist/MUC/LayerEff_Dimu", htmp).isSuccess()) {
00713 htmp->Fill( i, eff[1] );
00714 htmp->SetBinError( i+1, effErr[1] );
00715 } else {
00716 log << MSG::ERROR << "Fail to retrieve LayerEff_Dimu" << endreq;
00717 }
00718
00719 } // End layer
00720
00721 // Box efficiency and noise ratio, strip efficiency
00722 for( int i=0; i<BOX_MAX; i++ )
00723 {
00724 m_ptrIdTr->SetBoxPos( i, &part, &segment, &layer );
00725 stripMax = m_ptrIdTr->GetStripMax( part, segment, layer );
00726
00727 digi[0][1] = effHit[0][1] = trackNum[0][1] = nosHit[0][1] = 0;
00728 digi[1][1] = effHit[1][1] = trackNum[1][1] = nosHit[1][1] = 0;
00729 for( int j=0; j<stripMax; j++ )
00730 {
00731 // For box
00732 digi[0][1] += m_recordAll[part][segment][layer][j][0];
00733 trackNum[0][1] += m_recordAll[part][segment][layer][j][1];
00734 effHit[0][1] += m_recordAll[part][segment][layer][j][2];
00735 nosHit[0][1] += m_recordAll[part][segment][layer][j][3];
00736
00737 digi[1][1] += m_recordDimu[part][segment][layer][j][0];
00738 trackNum[1][1] += m_recordDimu[part][segment][layer][j][1];
00739 effHit[1][1] += m_recordDimu[part][segment][layer][j][2];
00740 nosHit[1][1] += m_recordDimu[part][segment][layer][j][3];
00741
00742 // For strip
00743 trackNum[0][2] = m_recordAll[part][segment][layer][j][1];
00744 effHit[0][2] = m_recordAll[part][segment][layer][j][2];
00745 trackNum[1][2] = m_recordDimu[part][segment][layer][j][1];
00746 effHit[1][2] = m_recordDimu[part][segment][layer][j][2];
00747
00748 if( trackNum[0][2] == 0 ) {
00749 eff[0] = 0.0;
00750 } else {
00751 eff[0] = ( (double)effHit[0][2] )/trackNum[0][2];
00752 }
00753 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_All", htmp).isSuccess()) {
00754 htmp->Fill( eff[0] );
00755 } else {
00756 log << MSG::ERROR << "Fail to retrieve StripEff_All" << endreq;
00757 }
00758
00759 if( trackNum[1][2] == 0 ) {
00760 eff[1] = 0.0;
00761 } else {
00762 eff[1] = ( (double)effHit[1][2] )/trackNum[1][2];
00763 }
00764 if(m_thsvc->getHist("/DQAHist/MUC/StripEff_Dimu", htmp).isSuccess()) {
00765 htmp->Fill( eff[1] );
00766 } else {
00767 log << MSG::ERROR << "Fail to retrieve StripEff_Dimu" << endreq;
00768 }
00769 } // End StripMax
00770
00771 // Box efficiency
00772 if( trackNum[0][1] == 0 ) {
00773 eff[0] = effErr[0] = 0.0;
00774 } else {
00775 eff[0] = ( (double)effHit[0][1] )/trackNum[0][1];
00776 effErr[0] = sqrt( eff[0]*(1-eff[0])/trackNum[0][1] );
00777 }
00778 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_All", htmp).isSuccess()) {
00779 htmp->Fill( i, eff[0] );
00780 htmp->SetBinError( i+1, effErr[0] );
00781 } else {
00782 log << MSG::ERROR << "Fail to retrieve BoxEff_All" << endreq;
00783 }
00784 log << MSG::INFO << "BoxEff:\t" << i << "\t" << eff[0] << "\t" << effErr[0] << endreq;
00785
00786 if( trackNum[1][1] == 0 ) {
00787 eff[1] = effErr[1] = 0.0;
00788 } else {
00789 eff[1] = ( (double)effHit[1][1] )/trackNum[1][1];
00790 effErr[1] = sqrt( eff[1]*(1-eff[1])/trackNum[1][1] );
00791 }
00792 if(m_thsvc->getHist("/DQAHist/MUC/BoxEff_Dimu", htmp).isSuccess()) {
00793 htmp->Fill( i, eff[1] );
00794 htmp->SetBinError( i+1, effErr[1] );
00795 } else {
00796 log << MSG::ERROR << "Fail to retrieve BoxEff_Dimu" << endreq;
00797 }
00798
00799 // Box noise ratio
00800 if( digi[0][1] == 0 ) {
00801 nos[0] = nosErr[0] = 0.0;
00802 } else {
00803 nos[0] = ( (double)nosHit[0][1] )/digi[0][1];
00804 nosErr[0] = sqrt( nos[0]*(1-nos[0])/digi[0][1] );
00805 }
00806 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_All", htmp).isSuccess()) {
00807 htmp->Fill( i, nos[0] );
00808 htmp->SetBinError( i+1, nosErr[0] );
00809 } else {
00810 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_All" << endreq;
00811 }
00812 log << MSG::INFO << "BoxNosRatio:\t" << i << "\t" << nos[0] << "\t" << nosErr[0] << endreq;
00813
00814 if( digi[1][1] == 0 ) {
00815 nos[1] = nosErr[1] = 0.0;
00816 } else {
00817 nos[1] = ( (double)nosHit[1][1] )/digi[1][1];
00818 nosErr[1] = sqrt( nos[1]*(1-nos[1])/digi[1][1] );
00819 }
00820 if(m_thsvc->getHist("/DQAHist/MUC/BoxNosRatio_Dimu", htmp).isSuccess()) {
00821 htmp->Fill( i, nos[1] );
00822 htmp->SetBinError( i+1, nosErr[1] );
00823 } else {
00824 log << MSG::ERROR << "Fail to retrieve BoxNosRatio_Dimu" << endreq;
00825 }
00826 }// End BOX_MAX
00827 */
00828
00829 log << MSG::INFO << "Finalize successfully! " << endmsg;
00830
00831 return StatusCode::SUCCESS;
00832 }
|
|
|
|
|
|
00044 {
00045 MsgStream log(msgSvc(), name());
00046
00047 log << MSG::INFO << "in initialize()" << endmsg;
00048 StatusCode status;
00049
00050 // Call Histogram service
00051 if(service("THistSvc", m_thsvc).isFailure()) {
00052 log << MSG::ERROR << "Couldn't get THistSvc" << endreq;
00053 return StatusCode::FAILURE;
00054 }
00055
00056 char name[60];
00057 char title[60];
00058
00059 // Resolution histograms
00060 for( int i=0; i<B_LAY_NUM; i++ )
00061 {
00062 sprintf( name, "BrResDist_All_L%d", i );
00063 sprintf( title, "Barrel spacial resolution in L%d",i );
00064 m_hBrResDist[i][0] = new TH1F(name,title, 200, -100, 100 );
00065 sprintf( name, "/DQAHist/MUC/BrResDist_All_L%d", i );
00066 if(m_thsvc->regHist(name, m_hBrResDist[i][0]).isFailure())
00067 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00068
00069 sprintf( name, "BrResDist_Dimu_L%d", i );
00070 m_hBrResDist[i][1] = new TH1F(name,title, 200, -100, 100 );
00071 sprintf( name, "/DQAHist/MUC/BrResDist_Dimu_L%d", i );
00072 if(m_thsvc->regHist(name, m_hBrResDist[i][1]).isFailure())
00073 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00074 }
00075
00076 for( int i=0; i<E_LAY_NUM; i++ )
00077 {
00078 sprintf( name, "EcResDist_All_L%d", i );
00079 sprintf( title, "Endcap spacial resolution in L%d",i );
00080 m_hEcResDist[i][0] = new TH1F(name,title, 200, -100, 100 );
00081 sprintf( name, "/DQAHist/MUC/EcResDist_All_L%d", i );
00082 if(m_thsvc->regHist(name, m_hEcResDist[i][0]).isFailure())
00083 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00084
00085 sprintf( name, "EcResDist_Dimu_L%d", i );
00086 m_hEcResDist[i][1] = new TH1F(name,title, 200, -100, 100 );
00087 sprintf( name, "/DQAHist/MUC/EcResDist_Dimu_L%d", i );
00088 if(m_thsvc->regHist(name, m_hEcResDist[i][1]).isFailure())
00089 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00090 }
00091
00092 //m_hBrRes[0] = new TH1F("BrRes_All", "BrResAll", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
00093 //m_hBrRes[0]->GetXaxis()->SetTitle("Layer id");
00094 //m_hBrRes[0]->GetXaxis()->CenterTitle();
00095 //m_hBrRes[0]->SetStats(0);
00096 //m_hBrRes[1] = new TH1F("BrRes_Dimu", "BrResDimu", B_LAY_NUM+1, -0.5, B_LAY_NUM+0.5 );
00097 //m_hBrRes[1]->GetXaxis()->SetTitle("Layer id");
00098 //m_hBrRes[1]->GetXaxis()->CenterTitle();
00099 //m_hBrRes[1]->SetStats(0);
00100 //
00101 //m_hEcRes[0] = new TH1F("EcRes_All", "EcResAll", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
00102 //m_hEcRes[0]->GetXaxis()->SetTitle("Layer id");
00103 //m_hEcRes[0]->GetXaxis()->CenterTitle();
00104 //m_hEcRes[0]->SetStats(0);
00105 //m_hEcRes[1] = new TH1F("EcRes_Dimu", "EcResDimu", E_LAY_NUM+1, -0.5, E_LAY_NUM+0.5 );
00106 //m_hEcRes[1]->GetXaxis()->SetTitle("Layer id");
00107 //m_hEcRes[1]->GetXaxis()->CenterTitle();
00108 //m_hEcRes[1]->SetStats(0);
00109
00110 //if(m_thsvc->regHist("/DQAHist/MUC/BrRes_All", m_hBrRes[0]).isFailure())
00111 //{ log << MSG::ERROR << "Couldn't register BrRes_All" << endreq; }
00112 //if(m_thsvc->regHist("/DQAHist/MUC/BrRes_Dimu", m_hBrRes[1]).isFailure())
00113 //{ log << MSG::ERROR << "Couldn't register BrRes_Dimu" << endreq; }
00114 //if(m_thsvc->regHist("/DQAHist/MUC/EcRes_All", m_hEcRes[0]).isFailure())
00115 //{ log << MSG::ERROR << "Couldn't register EcRes_All" << endreq; }
00116 //if(m_thsvc->regHist("/DQAHist/MUC/EcRes_Dimu", m_hEcRes[1]).isFailure())
00117 //{ log << MSG::ERROR << "Couldn't register BrRes_Dimu" << endreq; }
00118
00119 // Efficiency histograms, j=0: numberator, j=1; denominator
00120 for(int i=0; i<TAGN; i++)
00121 {
00122 for(int j=0; j<2; j++)
00123 {
00124 for(int k=0; k<LVLN; k++)
00125 {
00126 sprintf( name, "%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
00127 sprintf( title, "%s efficiency",LNAME[k] );
00128
00129 if(k==0) {
00130 m_hEff[k][i][j] = new TH1F(name,title, LAYER_MAX+1, -0.5, LAYER_MAX+0.5);
00131 m_hEff[k][i][j]->GetXaxis()->SetTitle("Layer id");
00132 }else if(k==1) {
00133 m_hEff[k][i][j] = new TH1F(name,title, BOX_MAX+1, -0.5, BOX_MAX+0.5);
00134 m_hEff[k][i][j]->GetXaxis()->SetTitle("Box id [EE:0~31, BR:32~103, WE:104~135]");
00135 }else {
00136 m_hEff[k][i][j] = new TH1F(name,title, STRIP_MAX+1, -0.5, STRIP_MAX+0.5);
00137 m_hEff[k][i][j]->GetXaxis()->SetTitle("Strip id");
00138 }
00139
00140 m_hEff[k][i][j]->GetXaxis()->CenterTitle();
00141 //m_hEff[k][i][j]->SetStats(0);
00142
00143 sprintf( name, "/DQAHist/MUC/%sEff_%s_%s", LNAME[k], ENAME[i], HTYPE[j] );
00144 if(m_thsvc->regHist(name, m_hEff[k][i][j]).isFailure())
00145 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00146 }
00147
00148 sprintf( name, "%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
00149 sprintf( title, "%s noise ratio",LNAME[1] );
00150 m_hNosRatio[i][j] = new TH1F(name, title, BOX_MAX+1, -0.5, BOX_MAX+0.5);
00151 m_hNosRatio[i][j]->GetXaxis()->SetTitle("Box id [EE:0~31, BR:32~103, WE:104~135]");
00152 m_hNosRatio[i][j]->GetXaxis()->CenterTitle();
00153 //m_hNosRatio[i][j]->GetYaxis()->SetRangeUser(0., 0.3);
00154 //m_hNosRatio[i][j]->SetStats(0);
00155
00156 sprintf( name, "/DQAHist/MUC/%sNosRatio_%s_%s", LNAME[1], ENAME[i], HTYPE[j] );
00157 if(m_thsvc->regHist(name, m_hNosRatio[i][j]).isFailure())
00158 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00159 }
00160
00161 // Costheta and phi
00162 sprintf( name, "Costheta_%s", ENAME[i]);
00163 sprintf( title, "Costheta");
00164 m_hCostheta[i] = new TH1F(name, title, 360, -1.0, 1.0);
00165 m_hCostheta[i]->GetXaxis()->SetTitle("cos#theta");
00166 m_hCostheta[i]->GetXaxis()->CenterTitle();
00167
00168 sprintf( name, "Phi_%s", ENAME[i]);
00169 sprintf( title, "Phi");
00170 m_hPhi[i] = new TH1F(name, title, 360, -PI, PI);
00171 m_hPhi[i]->GetXaxis()->SetTitle("#phi");
00172 m_hPhi[i]->GetXaxis()->CenterTitle();
00173
00174 sprintf( name, "/DQAHist/MUC/Costheta_%s", ENAME[i] );
00175 if(m_thsvc->regHist(name, m_hCostheta[i]).isFailure())
00176 { log << MSG::ERROR << "Couldn't register " << name << endreq; }
00177 sprintf( name, "/DQAHist/MUC/Phi_%s", ENAME[i] );
00178 if(m_thsvc->regHist(name, m_hPhi[i]).isFailure())
00179 { log << MSG::ERROR << "Couldn't register Phi_All" << endreq; }
00180 }
00181
00182 //m_hStripEff[0] = new TH1F("StripEff_All", "Strip efficiency", 101, -0.005, 1.005 );
00183 //m_hStripEff[0]->GetXaxis()->SetTitle("Efficiency");
00184 //m_hStripEff[0]->GetXaxis()->CenterTitle();
00185 //m_hStripEff[1] = new TH1F("StripEff_Dimu", "Strip efficiency", 101, -0.005, 1.005 );
00186 //m_hStripEff[1]->GetXaxis()->SetTitle("Efficiency");
00187 //m_hStripEff[1]->GetXaxis()->CenterTitle();
00188
00189 // Parameters
00190 for( int i=0; i<PART_MAX; i++ )
00191 for( int j=0; j<SEGMENT_MAX; j++ )
00192 for( int k=0; k<LAYER_MAX; k++ )
00193 for( int n=0; n<STRIP_INBOX_MAX; n++ )
00194 {
00195 m_recordAll[i][j][k][n][0] = 0;
00196 m_recordAll[i][j][k][n][1] = 0;
00197 m_recordAll[i][j][k][n][2] = 0;
00198 m_recordAll[i][j][k][n][3] = 0;
00199 m_recordAll[i][j][k][n][4] = 0;
00200
00201 m_recordDimu[i][j][k][n][0] = 0;
00202 m_recordDimu[i][j][k][n][1] = 0;
00203 m_recordDimu[i][j][k][n][2] = 0;
00204 m_recordDimu[i][j][k][n][3] = 0;
00205 m_recordDimu[i][j][k][n][4] = 0;
00206 }
00207
00208 m_ptrMucMark = new MucMark(0,0,0,0);
00209 m_ptrIdTr = new MucIdTransform();
00210
00211 m_digiCol.resize(0);
00212 m_expHitCol.resize(0);
00213 m_effHitCol.resize(0);
00214 m_nosHitCol.resize(0);
00215 m_clusterCol.resize(0);
00216
00217 for( int i=0; i<PART_MAX; i++ )
00218 for( int j=0; j<SEGMENT_MAX; j++ ) {
00219 m_segDigiCol[i][j].resize(0);
00220 }
00221
00222 log << MSG::INFO << "Initialize done!" <<endmsg;
00223 return StatusCode::SUCCESS;
00224 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1