ResiAlign.cxx

Go to the documentation of this file.

#include "MdcAlignAlg/ResiAlign.h"
#include "MdcAlignAlg/Alignment.h"
#include "MdcAlignAlg/MdcAlignAlg.h"
 
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/IMessageSvc.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/IService.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/StatusCode.h"
 
#include "EventModel/Event.h"
#include "EventModel/EventHeader.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "MdcRawEvent/MdcDigi.h"
 
#include "TCanvas.h"
#include "TF1.h"
 
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <string>
 
using namespace std;
 
ResiAlign::ResiAlign() {
  // m_ndiv = 6;
  m_ndiv    = 12;
  m_resiCut = 1.0;
  for ( int i = 0; i < NEP; i++ ) m_npoint[i] = 0;
 
  for ( int lay = 0; lay < LAYERNMAX; lay++ )
  {
    if ( lay < 8 )
    {
      m_docaMin[lay] = 0.6; // mm
      m_docaMax[lay] = 4.8; // mm
    }
    else
    {
      m_docaMin[lay] = 0.8; // mm
      m_docaMax[lay] = 6.4; // mm
    }
  }
  for ( int lay = 0; lay < LAYERNMAX; lay++ )
  {
    if ( ( 0 == lay ) || ( 7 == lay ) || ( 8 == lay ) || ( 19 == lay ) || ( 20 == lay ) ||
         ( 35 == lay ) || ( 36 == lay ) || ( 42 == lay ) )
      m_layBound[lay] = true;
    else m_layBound[lay] = false;
  }
 
  m_ncut1  = 0;
  m_ncut2  = 0;
  m_ncut3  = 0;
  m_ncut4  = 0;
  m_ncut5  = 0;
  m_ncut6  = 0;
  m_ncut7  = 0;
  m_ncut8  = 0;
  m_ncut9  = 0;
  m_ncut10 = 0;
  m_ncut11 = 0;
  m_ncut12 = 0;
  m_ncut13 = 0;
}
 
ResiAlign::~ResiAlign() {}
 
void ResiAlign::clear() {
  delete m_hresAll;
  delete m_hresInn;
  delete m_hresStp;
  delete m_hresOut;
  for ( int lay = 0; lay < LAYERNMAX; lay++ )
  {
    delete m_hresLay[lay];
    for ( int i = 0; i < 4; i++ ) delete m_hresLay_LR[lay][i];
  }
  for ( int i = 0; i < NEP; i++ ) delete m_gr[i];
}
 
void ResiAlign::initialize( TObjArray* hlist, IMdcGeomSvc* mdcGeomSvc,
                            IMdcCalibFunSvc* mdcFunSvc, IMdcUtilitySvc* mdcUtilitySvc ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "ResiAlign" );
  log << MSG::INFO << "ResiAlign::initialize()" << endmsg;
 
  m_hlist         = hlist;
  m_mdcGeomSvc    = mdcGeomSvc;
  m_mdcFunSvc     = mdcFunSvc;
  m_mdcUtilitySvc = mdcUtilitySvc;
 
  double zeast;
  for ( int lay = 0; lay < 43; lay++ )
  {
    zeast            = m_mdcGeomSvc->Wire( lay, 0 )->Backward().z();
    m_zrange[lay][1] = 2.0 * fabs( zeast ) / (double)m_ndiv;
    m_zrange[lay][0] = -1.0 * m_zrange[lay][1];
 
    m_radii[lay] = m_mdcGeomSvc->Layer( lay )->Radius();
  }
 
  for ( int wir = 0; wir < WIRENMAX; wir++ )
  {
    m_xe[wir] = m_mdcGeomSvc->Wire( wir )->Backward().x();
    m_ye[wir] = m_mdcGeomSvc->Wire( wir )->Backward().y();
    m_ze[wir] = m_mdcGeomSvc->Wire( wir )->Backward().z();
    m_xw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().x();
    m_yw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().y();
    m_zw[wir] = m_mdcGeomSvc->Wire( wir )->Forward().z();
  }
 
  char hname[200];
  int  iEP;
 
  INTupleSvc* ntupleSvc;
  Gaudi::svcLocator()->service( "NTupleSvc", ntupleSvc );
  for ( iEP = 0; iEP <= NEP; iEP++ )
  {
    if ( iEP < NEP ) sprintf( hname, "FILE137/align%02d", iEP );
    else sprintf( hname, "FILE137/alignAll" );
 
    NTuplePtr nt( ntupleSvc, hname );
    if ( nt ) m_tuple[iEP] = nt;
    else
    {
      m_tuple[iEP] = ntupleSvc->book( hname, CLID_ColumnWiseTuple, "align" );
      if ( m_tuple[iEP] )
      {
        m_tuple[iEP]->addItem( "run", m_iRun[iEP] );
        m_tuple[iEP]->addItem( "evt", m_iEvt[iEP] );
        m_tuple[iEP]->addItem( "resi", m_resi[iEP] );
        m_tuple[iEP]->addItem( "p", m_p[iEP] );
        m_tuple[iEP]->addItem( "pt", m_pt[iEP] );
        m_tuple[iEP]->addItem( "phi", m_phi[iEP] );
        m_tuple[iEP]->addItem( "lay", m_lay[iEP] );
        m_tuple[iEP]->addItem( "lr", m_lr[iEP] );
        m_tuple[iEP]->addItem( "cel", m_cel[iEP] );
      }
      else { log << MSG::FATAL << "Cannot book N-tuple:" << long( m_tuple[iEP] ) << endmsg; }
    }
  }
 
  m_hnTrk = new TH1F( "HNtrack", "", 10, -0.5, 9.5 );
  m_hlist->Add( m_hnTrk );
 
  m_hnHit = new TH1F( "HNhit", "", 100, -0.5, 99.5 );
  m_hlist->Add( m_hnHit );
 
  m_hlayHitmap = new TH1F( "Hitmap", "", 43, -0.5, 42.5 );
  m_hlist->Add( m_hlayHitmap );
 
  m_hresAll = new TH1F( "HResAllInc", "", 200, -1.0, 1.0 );
  m_hlist->Add( m_hresAll );
 
  m_hresInn = new TH1F( "HResInnInc", "", 200, -1.0, 1.0 );
  m_hlist->Add( m_hresInn );
 
  m_hresStp = new TH1F( "HResStpInc", "", 200, -1.0, 1.0 );
  m_hlist->Add( m_hresStp );
 
  m_hresOut = new TH1F( "HResOutInc", "", 200, -1.0, 1.0 );
  m_hlist->Add( m_hresOut );
 
  int lay;
  for ( lay = 0; lay < LAYERNMAX; lay++ )
  {
    sprintf( hname, "Res_Layer%02d", lay );
    m_hresLay[lay] = new TH1F( hname, "", 200, -1.0, 1.0 );
    m_hlist->Add( m_hresLay[lay] );
 
    for ( int i = 0; i < 4; i++ )
    {
      if ( 0 == i ) sprintf( hname, "Resi_Lay%02d_Up_L", lay );
      else if ( 1 == i ) sprintf( hname, "Resi_Lay%02d_Up_R", lay );
      else if ( 2 == i ) sprintf( hname, "Resi_Lay%02d_Dw_L", lay );
      else sprintf( hname, "Resi_Lay%02d_Dw_R", lay );
      m_hresLay_LR[lay][i] = new TH1F( hname, "", 200, -1.0, 1.0 );
      m_hlist->Add( m_hresLay_LR[lay][i] );
    }
  }
 
  for ( iEP = 0; iEP < NEP; iEP++ )
  {
    m_gr[iEP] = new TGraph();
    sprintf( hname, "grResi%02d", iEP );
    m_gr[iEP]->SetName( hname );
    m_hlist->Add( m_gr[iEP] );
  }
  m_fevt.open( "evt.txt" );
}
 
bool ResiAlign::fillHist( MdcAliEvent* event ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "ResiAlign" );
  log << MSG::DEBUG << "ResiAlign::fillHist()" << endmsg;
 
  bool esCutFg = event->getEsCutFlag();
  if ( !esCutFg )
  {
    m_ncut1++;
    return true;
  }
 
  int i = 0;
  int k;
 
  int    trkStat;
  double dr;
  double phi0;
  double kappa;
  double dz;
  double tanl;
  double chisq;
  double p;
  double pt;
 
  int nhits;
  int lay;
  int cel;
  int wir;
  int lr;
  int iEnd;
  int iEP;
 
  double        doca;
  double        resi;
  double        zhit;
  double        wphi;
  double        dphi;
  double        hitPhi;
  double        xx;
  double        yy;
  double        rr;
  int           stat;
  MdcAliRecTrk* rectrk;
  MdcAliRecHit* rechit;
  int           nhitlay;
  bool          fgHitLay[LAYERNMAX];
 
  IDataProviderSvc* eventSvc = NULL;
  Gaudi::svcLocator()->service( "EventDataSvc", eventSvc );
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc, "/Event/EventHeader" );
  if ( !eventHeader )
  {
    log << MSG::FATAL << "Could not find Event Header" << endmsg;
    m_ncut3++;
    return ( StatusCode::FAILURE );
  }
  int iEvt = eventHeader->eventNumber();
  int iRun = eventHeader->runNumber();
 
  int nTrk = event->getNTrk();
  m_hnTrk->Fill( nTrk );
  m_fevt << setw( 10 ) << iRun << setw( 10 ) << iEvt << setw( 10 ) << nTrk << endl;
  if ( ( nTrk < m_param.nTrkCut[0] ) || ( nTrk > m_param.nTrkCut[1] ) )
  {
    m_ncut2++;
    return true;
  }
 
  for ( i = 0; i < nTrk; i++ )
  {
    rectrk  = event->getRecTrk( i );
    nhits   = rectrk->getNHits();
    trkStat = rectrk->getStat();
    //    if (0 != trkStat) continue;
 
    // dr cut
    dr = rectrk->getDr();
    if ( fabs( dr ) > m_param.drCut )
    {
      m_ncut4++;
      continue;
    }
 
    phi0  = rectrk->getPhi0();
    kappa = rectrk->getKappa();
 
    // dz cut
    dz = rectrk->getDz();
    if ( fabs( dz ) > m_param.dzCut )
    {
      m_ncut5++;
      continue;
    }
 
    for ( lay = 0; lay < LAYERNMAX; lay++ ) { fgHitLay[lay] = false; }
 
    m_hnHit->Fill( nhits );
    for ( k = 0; k < nhits; k++ )
    {
      rechit        = rectrk->getRecHit( k );
      lay           = rechit->getLayid();
      fgHitLay[lay] = true;
    }
 
    nhitlay = 0;
    for ( lay = 0; lay < LAYERNMAX; lay++ )
    {
      if ( fgHitLay[lay] ) nhitlay++;
    }
    //    cout << "hitlay " << nhitlay << endl;
    if ( nhitlay < m_param.nHitLayCut )
    {
      m_ncut6++;
      continue;
    }
 
    tanl  = rectrk->getTanLamda();
    chisq = rectrk->getChisq();
    p     = rectrk->getP();
    pt    = rectrk->getPt();
 
    if ( ( fabs( pt ) < m_param.ptCut[0] ) || ( fabs( pt ) > m_param.ptCut[1] ) )
    {
      m_ncut7++;
      continue;
    }
 
    HepVector helix = rectrk->getHelix();
    for ( k = 0; k < nhits; k++ )
    {
      rechit = rectrk->getRecHit( k );
      lay    = rechit->getLayid();
      cel    = rechit->getCellid();
      lr     = rechit->getLR();
      doca   = rechit->getDocaInc();
      zhit   = rechit->getZhit();
 
      stat = rechit->getStat();
      if ( ( 1 == m_param.hitStatCut ) && ( 1 != stat ) )
      {
        m_ncut8++;
        continue;
      }
 
      if ( 1 == m_param.resiType ) { resi = rechit->getResiExcLR(); }
      else { resi = rechit->getResiIncLR(); }
      resi *= -1.0;
      // maqm if( (1==isnan(resi)) || (fabs(resi) > m_resiCut) ||
      if ( ( 1 == std::isnan( resi ) ) || ( fabs( resi ) > m_resiCut ) ||
           ( fabs( doca ) > m_docaMax[lay] ) || ( fabs( doca ) < m_docaMin[lay] ) )
      {
        m_ncut9++;
        continue;
      }
 
      if ( m_param.fgAdjacLayerCut )
      {
        if ( 0 == lay )
        {
          if ( !fgHitLay[1] )
          {
            m_ncut10++;
            continue;
          }
        }
        else if ( 42 == lay )
        {
          if ( !fgHitLay[41] )
          {
            m_ncut11++;
            continue;
          }
        }
        else
        {
          if ( ( !fgHitLay[lay - 1] ) && ( !fgHitLay[lay + 1] ) )
          {
            m_ncut12++;
            continue;
          }
          // for boundary layers
          if ( m_param.fgBoundLayerCut && m_layBound[lay] &&
               ( ( !fgHitLay[lay - 1] ) || ( !fgHitLay[lay + 1] ) ) )
          {
            m_ncut13++;
            continue;
          }
        }
      }
      m_hlayHitmap->Fill( lay );
 
      // fill alignment trees
      if ( ( zhit < m_zrange[lay][0] ) || ( zhit > m_zrange[lay][1] ) )
      {
        wir = m_mdcGeomSvc->Wire( lay, cel )->Id();
        xx  = ( zhit - m_zw[wir] ) * ( m_xe[wir] - m_xw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
             m_xw[wir];
        yy = ( zhit - m_zw[wir] ) * ( m_ye[wir] - m_yw[wir] ) / ( m_ze[wir] - m_zw[wir] ) +
             m_yw[wir];
        rr   = sqrt( ( xx * xx ) + ( yy * yy ) );
        dphi = fabs( doca ) / m_radii[lay];
 
        if ( yy >= 0 ) wphi = acos( xx / rr );
        else wphi = PI2 - acos( xx / rr );
        if ( 1 == lr ) hitPhi = wphi + dphi;
        else hitPhi = wphi - dphi;
        if ( hitPhi < 0 ) hitPhi += PI2;
        else if ( hitPhi > PI2 ) hitPhi -= PI2;
 
        if ( zhit < m_zrange[lay][0] ) iEnd = 0; // west
        else iEnd = 1;                           // east
        iEP = Alignment::getEpId( lay, iEnd );
 
        m_iRun[iEP] = iRun;
        m_iEvt[iEP] = iEvt;
        m_resi[iEP] = resi;
        m_p[iEP]    = p;
        m_pt[iEP]   = pt;
        m_phi[iEP]  = hitPhi;
        m_lay[iEP]  = lay;
        m_lr[iEP]   = lr;
        m_cel[iEP]  = cel;
        m_tuple[iEP]->write();
 
        m_resi[NEP] = resi;
        m_p[NEP]    = p;
        m_pt[NEP]   = pt;
        m_phi[NEP]  = hitPhi;
        m_lay[NEP]  = lay;
        m_lr[NEP]   = lr;
        m_cel[NEP]  = cel;
        m_tuple[NEP]->write();
 
        m_hresAll->Fill( resi );
        if ( lay < 8 ) m_hresInn->Fill( resi );
        else if ( lay < 20 ) m_hresStp->Fill( resi );
        else m_hresOut->Fill( resi );
        m_hresLay[lay]->Fill( resi );
        if ( 0 == lr )
        {
          if ( hitPhi < 3.14 ) m_hresLay_LR[lay][0]->Fill( resi ); // up L
          else m_hresLay_LR[lay][2]->Fill( resi );                 // down L
        }
        else
        {
          if ( hitPhi < 3.14 ) m_hresLay_LR[lay][1]->Fill( resi ); // up R
          else m_hresLay_LR[lay][3]->Fill( resi );                 // down R
        }
 
        m_gr[iEP]->SetPoint( m_npoint[iEP], hitPhi, resi );
        m_npoint[iEP]++;
      }
    }
  }
 
  return true;
}
 
void ResiAlign::updateConst( MdcAlignPar* alignPar ) {
  IMessageSvc* msgSvc;
  Gaudi::svcLocator()->service( "MessageSvc", msgSvc );
  MsgStream log( msgSvc, "ResiAlign" );
  log << MSG::INFO << "ResiAlign::updateConst()" << endmsg;
  m_fevt.close();
 
  int    iEP;
  double par[3];
  double err[3];
  double dx;
  double dy;
  double rz;
  double rLayer[] = { 120.225, 205.0, 237.55, 270.175, 302.625, 334.775, 366.65, 500.0,
                      120.225, 205.0, 237.55, 270.175, 302.625, 334.775, 366.65, 500.0 };
 
  TCanvas c1( "c1", "c1", 10, 10, 700, 500 );
 
  TF1* fResPhi = new TF1( "fResPhi", funResi, 0, PI2, 3 );
  fResPhi->SetParameter( 0, 0.0 );
  fResPhi->SetParameter( 1, 0.0 );
  fResPhi->SetParameter( 2, 0.0 );
 
  for ( iEP = 0; iEP < NEP; iEP++ )
  {
    if ( ( m_gr[iEP]->GetN() ) > 500 )
    {
      m_gr[iEP]->Fit( "fResPhi", "V" );
      par[0] = fResPhi->GetParameter( 0 );
      par[1] = fResPhi->GetParameter( 1 );
      par[2] = fResPhi->GetParameter( 2 );
 
      err[0] = fResPhi->GetParError( 0 );
      err[1] = fResPhi->GetParError( 1 );
      err[2] = fResPhi->GetParError( 2 );
 
      dx = -1.0 * par[1];
      dy = par[2];
      rz = par[0] / rLayer[iEP];
 
      // assume the shift of the outer section is 0
      if ( 7 == iEP || 15 == iEP )
      {
        dx     = 0.0;
        dy     = 0.0;
        rz     = 0.0;
        par[0] = 0.0;
        par[1] = 0.0;
        par[2] = 0.0;
      }
      alignPar->setDelDx( iEP, dx );
      alignPar->setDelDy( iEP, dy );
      alignPar->setDelRz( iEP, rz );
 
      alignPar->setErrDx( iEP, err[1] );
      alignPar->setErrDy( iEP, err[2] );
      alignPar->setErrRz( iEP, err[0] / rLayer[iEP] );
    }
  }
 
  cout << "TrackCut: cut1: " << m_ncut1 << ",  cut2: " << m_ncut2 << ",  cut3: " << m_ncut3
       << ",  cut4: " << m_ncut4 << ",  cut5: " << m_ncut5 << ",  cut6: " << m_ncut6
       << ",  cut7: " << m_ncut7 << endl;
  cout << "HitCut:   cut8: " << m_ncut8 << ",  cut9: " << m_ncut9 << ",  cut10: " << m_ncut10
       << ",  cut11: " << m_ncut11 << ",  cut12: " << m_ncut12 << ",  cut13: " << m_ncut13
       << endl;
 
  delete fResPhi;
}
 
Double_t ResiAlign::funResi( Double_t* x, Double_t* par ) {
  Double_t val;
  val = par[0] + par[1] * sin( x[0] ) + par[2] * cos( x[0] );
  return val;
}