BOSS_7.0.2: /home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Mdc/MdcAlignAlg/MdcAlignAlg-00-01-04/src/ResiAlign.cxx Source File

00001 #include "MdcAlignAlg/ResiAlign.h"
00002 #include "MdcAlignAlg/MdcAlignAlg.h"
00003 #include "MdcAlignAlg/Alignment.h"
00004 
00005 #include "GaudiKernel/INTupleSvc.h"
00006 #include "GaudiKernel/IService.h"
00007 #include "GaudiKernel/IMessageSvc.h"
00008 #include "GaudiKernel/StatusCode.h"
00009 #include "GaudiKernel/ISvcLocator.h"
00010 #include "GaudiKernel/Bootstrap.h"
00011 
00012 #include "EventModel/Event.h"
00013 #include "EventModel/EventHeader.h"
00014 #include "MdcRawEvent/MdcDigi.h"
00015 #include "Identifier/Identifier.h"
00016 #include "Identifier/MdcID.h"
00017 
00018 #include "TF1.h" 
00019 #include "TCanvas.h"
00020 
00021 #include <iostream>
00022 #include <fstream>
00023 #include <iomanip>
00024 #include <string>
00025 #include <cstring>
00026 
00027 using namespace std;
00028 
00029 ResiAlign::ResiAlign(){
00030      //m_ndiv = 6;
00031      m_ndiv = 12;
00032      m_resiCut = 1.0;
00033      for(int i=0; i<NEP; i++) m_npoint[i] = 0;
00034 
00035      for(int lay=0; lay<LAYERNMAX; lay++){
00036           if(lay < 8){
00037                m_docaMin[lay] = 0.6; // mm
00038                m_docaMax[lay] = 4.8; // mm
00039           } else{
00040                m_docaMin[lay] = 0.8; // mm
00041                m_docaMax[lay] = 6.4; // mm
00042           }
00043      }
00044      for(int lay=0; lay<LAYERNMAX; lay++){
00045           if((0==lay) || (7==lay) || (8==lay) || (19==lay) || (20==lay) ||
00046              (35==lay) || (36==lay) || (42==lay) ) m_layBound[lay] = true;
00047           else m_layBound[lay] = false;
00048      }
00049 
00050      m_ncut1 = 0;
00051      m_ncut2 = 0;
00052      m_ncut3 = 0;
00053      m_ncut4 = 0;
00054      m_ncut5 = 0;
00055      m_ncut6 = 0;
00056      m_ncut7 = 0;
00057      m_ncut8 = 0;
00058      m_ncut9 = 0;
00059      m_ncut10 = 0;
00060      m_ncut11 = 0;
00061      m_ncut12 = 0;
00062      m_ncut13 = 0;
00063 }
00064 
00065 ResiAlign::~ResiAlign(){
00066 }
00067 
00068 void ResiAlign::clear(){
00069      delete m_hresAll;
00070      delete m_hresInn;
00071      delete m_hresStp;
00072      delete m_hresOut;
00073      for(int lay=0; lay<LAYERNMAX; lay++) delete m_hresLay[lay];
00074      for(int i=0; i<NEP; i++) delete m_gr[i];
00075 }
00076 
00077 void ResiAlign::initialize(TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
00078                            IMdcCalibFunSvc* mdcFunSvc){
00079      IMessageSvc* msgSvc;
00080      Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00081      MsgStream log(msgSvc, "ResiAlign");
00082      log << MSG::INFO << "ResiAlign::initialize()" << endreq;
00083 
00084      m_hlist = hlist;
00085      m_mdcGeomSvc = mdcGeomSvc;
00086      m_mdcFunSvc = mdcFunSvc;
00087 
00088      double zeast;
00089      for(int lay=0; lay<43; lay++){
00090           zeast = m_mdcGeomSvc->Wire(lay, 0)->Backward().z();
00091           m_zrange[lay][1] = 2.0 * fabs(zeast) / (double)m_ndiv;
00092           m_zrange[lay][0] = -1.0 * m_zrange[lay][1];
00093 
00094           m_radii[lay] = m_mdcGeomSvc->Layer(lay)->Radius();
00095      }
00096 
00097      for(int wir=0; wir<WIRENMAX; wir++){
00098           m_xe[wir] = m_mdcGeomSvc->Wire(wir)->Backward().x();
00099           m_ye[wir] = m_mdcGeomSvc->Wire(wir)->Backward().y();
00100           m_ze[wir] = m_mdcGeomSvc->Wire(wir)->Backward().z();
00101           m_xw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().x();
00102           m_yw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().y();
00103           m_zw[wir] = m_mdcGeomSvc->Wire(wir)->Forward().z();
00104      }
00105 
00106      char hname[200];
00107      int iEP;
00108 
00109      INTupleSvc* ntupleSvc;
00110      Gaudi::svcLocator() -> service("NTupleSvc", ntupleSvc);
00111      for(iEP=0; iEP<=NEP; iEP++){
00112           if(iEP < NEP) sprintf(hname, "FILE137/align%02d", iEP);
00113           else sprintf(hname, "FILE137/alignAll");
00114 
00115           NTuplePtr nt(ntupleSvc, hname);
00116           if( nt ) m_tuple[iEP] = nt;
00117           else{
00118                m_tuple[iEP] = ntupleSvc->book(hname, CLID_ColumnWiseTuple,"align");
00119                if (m_tuple[iEP]) {
00120                     m_tuple[iEP]->addItem ("run", m_iRun[iEP]);
00121                     m_tuple[iEP]->addItem ("evt", m_iEvt[iEP]);
00122                     m_tuple[iEP]->addItem ("resi", m_resi[iEP]);
00123                     m_tuple[iEP]->addItem ("p", m_p[iEP]);
00124                     m_tuple[iEP]->addItem ("pt", m_pt[iEP]);
00125                     m_tuple[iEP]->addItem ("phi", m_phi[iEP]);
00126                     m_tuple[iEP]->addItem ("lay", m_lay[iEP]);
00127                     m_tuple[iEP]->addItem ("lr", m_lr[iEP]);
00128                     m_tuple[iEP]->addItem ("cel", m_cel[iEP]);
00129                }
00130                else {
00131                     log << MSG::FATAL << "Cannot book N-tuple:"
00132                         << long(m_tuple[iEP]) << endmsg;
00133                }
00134           }
00135      }
00136 
00137      m_hnTrk = new TH1F("HNtrack", "", 10, -0.5, 9.5);
00138      m_hlist->Add(m_hnTrk);
00139 
00140      m_hnHit = new TH1F("HNhit", "", 100, -0.5, 99.5);
00141      m_hlist->Add(m_hnHit);
00142 
00143      m_hlayHitmap = new TH1F("Hitmap", "", 43, -0.5, 42.5);
00144      m_hlist->Add(m_hlayHitmap);
00145 
00146      m_hresAll = new TH1F("HResAllInc", "", 200, -1.0, 1.0);
00147      m_hlist->Add(m_hresAll);
00148 
00149      m_hresInn = new TH1F("HResInnInc", "", 200, -1.0, 1.0);
00150      m_hlist->Add(m_hresInn);
00151 
00152      m_hresStp = new TH1F("HResStpInc", "", 200, -1.0, 1.0);
00153      m_hlist->Add(m_hresStp);
00154 
00155      m_hresOut = new TH1F("HResOutInc", "", 200, -1.0, 1.0);
00156      m_hlist->Add(m_hresOut);
00157 
00158      int lay;
00159      for(lay=0; lay<LAYERNMAX; lay++){
00160           sprintf(hname, "Res_Layer%02d", lay);
00161           m_hresLay[lay] = new TH1F(hname, "", 200, -1.0, 1.0);
00162           m_hlist->Add(m_hresLay[lay]);
00163      }
00164 
00165      for(iEP=0; iEP<NEP; iEP++){
00166           m_gr[iEP] = new TGraph();
00167           sprintf(hname, "grResi%02d", iEP);
00168           m_gr[iEP]->SetName(hname);
00169           m_hlist->Add(m_gr[iEP]);
00170      }
00171      m_fevt.open("evt.txt");
00172 }
00173 
00174 bool ResiAlign::fillHist(MdcAliEvent* event){
00175      IMessageSvc* msgSvc;
00176      Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00177      MsgStream log(msgSvc, "ResiAlign");
00178      log << MSG::DEBUG << "ResiAlign::fillHist()" << endreq;
00179 
00180      bool esCutFg = event->getEsCutFlag();
00181      if( ! esCutFg ){
00182           m_ncut1++;
00183           return true;
00184      }
00185 
00186      int i = 0;
00187      int k;
00188 
00189      int trkStat;
00190      double dr;
00191      double phi0;
00192      double kappa;
00193      double dz;
00194      double tanl;
00195      double chisq;
00196      double p;
00197      double pt;
00198 
00199      int nhits;
00200      int lay;
00201      int cel;
00202      int wir;
00203      int lr;
00204      int iEnd;
00205      int iEP;
00206 
00207      double doca;
00208      double resi;
00209      double zhit;
00210      double wphi;
00211      double dphi;
00212      double hitPhi;
00213      double xx;
00214      double yy;
00215      double rr;
00216      int stat;
00217      MdcAliRecTrk* rectrk;
00218      MdcAliRecHit* rechit;
00219      int nhitlay;
00220      bool fgHitLay[LAYERNMAX];
00221 
00222      IDataProviderSvc* eventSvc = NULL;
00223      Gaudi::svcLocator()->service("EventDataSvc", eventSvc);
00224      SmartDataPtr<Event::EventHeader> eventHeader(eventSvc,"/Event/EventHeader");
00225      if (!eventHeader) {
00226           log << MSG::FATAL << "Could not find Event Header" << endreq;
00227           m_ncut3++;
00228           return( StatusCode::FAILURE);
00229      }
00230      int iEvt = eventHeader->eventNumber();
00231      int iRun = eventHeader->runNumber();
00232 
00233      int nTrk = event -> getNTrk();
00234      m_hnTrk->Fill(nTrk);
00235      m_fevt << setw(10) << iRun << setw(10) << iEvt << setw(10) << nTrk << endl;
00236      if((nTrk < m_param.nTrkCut[0]) || (nTrk > m_param.nTrkCut[1])){
00237           m_ncut2++;
00238           return true;
00239      }
00240 
00241      for(i=0; i<nTrk; i++){
00242           rectrk = event->getRecTrk(i);
00243           nhits = rectrk->getNHits();
00244           trkStat = rectrk->getStat();
00245 //        if (0 != trkStat) continue;
00246 
00247           // dr cut
00248           dr = rectrk -> getDr();
00249           if(fabs(dr) > m_param.drCut){ m_ncut4++; continue; }
00250 
00251           phi0 = rectrk -> getPhi0();
00252           kappa = rectrk -> getKappa();
00253 
00254           // dz cut
00255           dz = rectrk -> getDz();
00256           if(fabs(dz) > m_param.dzCut){ m_ncut5++; continue; }
00257 
00258           for(lay=0; lay<LAYERNMAX; lay++){
00259                fgHitLay[lay] = false;
00260           }
00261 
00262           m_hnHit->Fill(nhits);
00263           for(k=0; k<nhits; k++){
00264                rechit = rectrk -> getRecHit(k);
00265                lay = rechit -> getLayid();
00266                fgHitLay[lay] = true;
00267           }
00268 
00269           nhitlay = 0;
00270           for(lay=0; lay<LAYERNMAX; lay++){
00271                if(fgHitLay[lay]) nhitlay++;
00272           }
00273           if(nhitlay < m_param.nHitLayCut){ m_ncut6++; continue; }
00274 
00275           tanl = rectrk -> getTanLamda();
00276           chisq = rectrk -> getChisq();
00277           p = rectrk -> getP();
00278           pt = rectrk -> getPt();
00279 
00280           if((fabs(pt)<m_param.ptCut[0]) || (fabs(pt)>m_param.ptCut[1])){ m_ncut7++; continue;}
00281 
00282           for(k=0; k<nhits; k++){
00283                rechit = rectrk->getRecHit(k);
00284                lay = rechit->getLayid();
00285                cel = rechit->getCellid();
00286                lr = rechit->getLR();
00287                doca = rechit -> getDocaInc();
00288                zhit = rechit->getZhit();
00289 
00290                stat = rechit -> getStat();
00291                if((1 == m_param.hitStatCut) && (1 != stat)){ m_ncut8++; continue; }
00292 
00293                if (1 == m_param.resiType) {
00294                     resi = rechit->getResiExcLR();
00295                } else {
00296                     resi = rechit->getResiIncLR();
00297                }
00298                resi *= -1.0;
00299                if( (1==isnan(resi)) || (fabs(resi) > m_resiCut) ||
00300                    (fabs(doca) > m_docaMax[lay]) || (fabs(doca) < m_docaMin[lay]) ){
00301                     m_ncut9++;
00302                     continue;
00303                }
00304 
00305                if(m_param.fgAdjacLayerCut){
00306                     if(0 == lay){
00307                          if( ! fgHitLay[1] ){ m_ncut10++; continue; }
00308                     } else if(42 == lay){
00309                          if( ! fgHitLay[41] ){ m_ncut11++; continue; }
00310                     } else{
00311                          if( (!fgHitLay[lay-1]) && (!fgHitLay[lay+1]) ){ m_ncut12++; continue; }
00312                          // for boundary layers
00313                          if( m_param.fgBoundLayerCut && m_layBound[lay] && 
00314                              ((!fgHitLay[lay-1]) || (!fgHitLay[lay+1])) ){
00315                               m_ncut13++;
00316                               continue;
00317                          }
00318                     }
00319                }
00320 
00321                m_hlayHitmap->Fill(lay);
00322 
00323                // fill alignment trees
00324                if((zhit < m_zrange[lay][0]) || (zhit > m_zrange[lay][1])){
00325                     wir = m_mdcGeomSvc -> Wire(lay, cel) -> Id();
00326                     xx = (zhit - m_zw[wir]) * (m_xe[wir] - m_xw[wir]) / 
00327                          (m_ze[wir] - m_zw[wir]) + m_xw[wir];
00328                     yy = (zhit - m_zw[wir]) * (m_ye[wir] - m_yw[wir]) / 
00329                          (m_ze[wir] - m_zw[wir]) + m_yw[wir];
00330                     rr = sqrt( (xx * xx) + (yy * yy) );
00331                     dphi = fabs(doca) / m_radii[lay];
00332 
00333                     if( yy >= 0 )  wphi = acos(xx / rr);
00334                     else           wphi = PI2 - acos(xx / rr);
00335                     if(1 == lr)    hitPhi = wphi + dphi; // mention
00336                     else hitPhi = wphi - dphi;
00337                     if(hitPhi < 0) hitPhi += PI2;
00338                     else if(hitPhi > PI2) hitPhi -= PI2;
00339 
00340                     if(zhit < m_zrange[lay][0]) iEnd = 0; // west
00341                     else iEnd = 1; // east
00342                     iEP = Alignment::getEpId(lay, iEnd);
00343 
00344                     m_iRun[iEP] = iRun;
00345                     m_iEvt[iEP] = iEvt;
00346                     m_resi[iEP] = resi;
00347                     m_p[iEP] = p;
00348                     m_pt[iEP] = pt;
00349                     m_phi[iEP] = hitPhi;
00350                     m_lay[iEP] = lay;
00351                     m_lr[iEP] = lr;
00352                     m_cel[iEP] = cel;
00353                     m_tuple[iEP]->write();
00354 
00355                     m_resi[NEP] = resi;
00356                     m_p[NEP] = p;
00357                     m_pt[NEP] = pt;
00358                     m_phi[NEP] = hitPhi;
00359                     m_lay[NEP] = lay;
00360                     m_lr[NEP] = lr;
00361                     m_cel[NEP] = cel;
00362                     m_tuple[NEP]->write();
00363 
00364                     m_hresAll->Fill(resi);
00365                     if(lay < 8) m_hresInn->Fill(resi);
00366                     else if(lay < 20) m_hresStp->Fill(resi);
00367                     else m_hresOut->Fill(resi);
00368                     m_hresLay[lay]->Fill(resi);
00369 
00370                     m_gr[iEP]->SetPoint(m_npoint[iEP], hitPhi, resi);
00371                     m_npoint[iEP]++;
00372                }
00373           }
00374      }
00375 
00376      return true;
00377 }
00378 
00379 void ResiAlign::updateConst(MdcAlignPar* alignPar){
00380      IMessageSvc* msgSvc;
00381      Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
00382      MsgStream log(msgSvc, "ResiAlign");
00383      log << MSG::INFO << "ResiAlign::updateConst()" << endreq;
00384      m_fevt.close();
00385 
00386      int iEP;
00387      double par[3];
00388      double err[3];
00389      double dx;
00390      double dy;
00391      double rz;
00392      double rLayer[] = { 120.225, 205.0, 237.55, 270.175,
00393                          302.625, 334.775, 366.65, 500.0,
00394                          120.225, 205.0, 237.55, 270.175,
00395                          302.625, 334.775, 366.65, 500.0 };
00396 
00397      TCanvas c1("c1", "c1", 10, 10, 700, 500);
00398 
00399      TF1* fResPhi = new TF1("fResPhi", funResi, 0, PI2, 3);
00400      fResPhi->SetParameter(0, 0.0);
00401      fResPhi->SetParameter(1, 0.0);
00402      fResPhi->SetParameter(2, 0.0);
00403 
00404      for(iEP=0; iEP<NEP; iEP++){
00405           if((m_gr[iEP]->GetN()) > 500){
00406                m_gr[iEP]->Fit("fResPhi", "V");
00407                par[0] = fResPhi->GetParameter(0);
00408                par[1] = fResPhi->GetParameter(1);
00409                par[2] = fResPhi->GetParameter(2);
00410 
00411                err[0] = fResPhi->GetParError(0);
00412                err[1] = fResPhi->GetParError(1);
00413                err[2] = fResPhi->GetParError(2);
00414 
00415                dx = -1.0 * par[1];
00416                dy = par[2];
00417                rz = par[0] / rLayer[iEP];
00418 
00419                // assume the shift of the outer section is 0
00420                if (7==iEP || 15==iEP) {
00421                     dx = 0.0;
00422                     dy = 0.0;
00423                     rz = 0.0;
00424                     par[0] = 0.0;
00425                     par[1] = 0.0;
00426                     par[2] = 0.0;
00427                }
00428                alignPar->setDelDx(iEP, dx);
00429                alignPar->setDelDy(iEP, dy);
00430                alignPar->setDelRz(iEP, rz);
00431 
00432                alignPar->setErrDx(iEP, err[1]);
00433                alignPar->setErrDy(iEP, err[2]);
00434                alignPar->setErrRz(iEP, err[0]/rLayer[iEP]);
00435           }
00436      }
00437 
00438      cout << "TrackCut: cut1: " << m_ncut1 << ",  cut2: " << m_ncut2 << ",  cut3: " << m_ncut3
00439           << ",  cut4: " << m_ncut4 << ",  cut5: " << m_ncut5 << ",  cut6: " << m_ncut6
00440           << ",  cut7: " << m_ncut7 << endl;
00441      cout << "HitCut:   cut8: " << m_ncut8 << ",  cut9: " << m_ncut9 << ",  cut10: " << m_ncut10
00442           << ",  cut11: " << m_ncut11 << ",  cut12: " << m_ncut12 << ",  cut13: " << m_ncut13 << endl;
00443 
00444      delete fResPhi;
00445 }
00446 
00447 Double_t ResiAlign::funResi(Double_t* x, Double_t* par){
00448      Double_t val;
00449      val = par[0] + par[1]*sin(x[0]) + par[2]*cos(x[0]);
00450      return val;
00451 }
00452