DedxCalibEvent.cxx

Go to the documentation of this file.

#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/StatusCode.h"
 
#include "EvTimeEvent/RecEsTime.h"
#include "EventModel/EventHeader.h"
#include "Identifier/Identifier.h"
#include "Identifier/MdcID.h"
#include "MdcRecEvent/RecMdcDedx.h"
#include "MdcRecEvent/RecMdcDedxHit.h"
#include "MdcRecEvent/RecMdcHit.h"
#include "MdcRecEvent/RecMdcKalHelixSeg.h"
#include "MdcRecEvent/RecMdcKalTrack.h"
#include "MdcRecEvent/RecMdcTrack.h"
 
#include "EventModel/Event.h"
#include "EventModel/EventModel.h"
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
 
#include <TMath.h>
 
#include "DedxCalibAlg/DedxCalibEvent.h"
DECLARE_COMPONENT( DedxCalibEvent )
 
typedef std::vector<int> Vint;
 
using namespace std;
using CLHEP::HepVector;
static int evt_count( 0 ), evt_threshold( 0 );
 
DedxCalibEvent::DedxCalibEvent( const std::string& name, ISvcLocator* pSvcLocator )
    : DedxCalib( name, pSvcLocator ) {
  // Declare the properties
  declareProperty( "CutWire", cut_wire = -1 );
  declareProperty( "Count", m_count = 1000000 );
  declareProperty( "Gap", m_gap = 1 );
}
 
void DedxCalibEvent::initializing() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibEvent::initializing()" << endmsg;
 
  StatusCode status;
  NTuplePtr  nt1( ntupleSvc(), "FILE100/n103" );
  if ( nt1 ) m_nt1 = nt1;
  else
  {
    m_nt1 = ntupleSvc()->book( "FILE100/n103", CLID_ColumnWiseTuple, "dEdx per track" );
    if ( m_nt1 )
    {
      status = m_nt1->addItem( "ptrk", m_ptrk );
      status = m_nt1->addItem( "ptrk_t", m_ptrk_t );
      status = m_nt1->addItem( "sintheta", m_sintheta );
      status = m_nt1->addItem( "costheta", m_costheta );
      status = m_nt1->addItem( "charge", m_charge );
      status = m_nt1->addItem( "runNO", m_runNO );
      status = m_nt1->addItem( "runFlag", m_runFlag );
      status = m_nt1->addItem( "evtNO", m_evtNO );
      status = m_nt1->addItem( "t0", m_t0 );
      status = m_nt1->addItem( "trackId", m_trackId );
      status = m_nt1->addItem( "poca_x", m_poca_x );
      status = m_nt1->addItem( "poca_y", m_poca_y );
      status = m_nt1->addItem( "poca_z", m_poca_z );
      status = m_nt1->addItem( "recalg", m_recalg );
      status = m_nt1->addItem( "nhit", m_nhit );
      status = m_nt1->addItem( "nhits", m_nhits );
      status = m_nt1->addItem( "usedhit", m_usedhit );
 
      status = m_nt1->addItem( "ndedxhit", m_nphlisthit, 0, 100 );
      status = m_nt1->addIndexedItem( "dEdx_hit", m_nphlisthit, m_dEdx_hit );
      status = m_nt1->addIndexedItem( "pathlength_hit", m_nphlisthit, m_pathlength_hit );
      status = m_nt1->addIndexedItem( "wid_hit", m_nphlisthit, m_wid_hit );
      status = m_nt1->addIndexedItem( "layid_hit", m_nphlisthit, m_layid_hit );
      status = m_nt1->addIndexedItem( "dd_in_hit", m_nphlisthit, m_dd_in_hit );
      status = m_nt1->addIndexedItem( "eangle_hit", m_nphlisthit, m_eangle_hit );
      status = m_nt1->addIndexedItem( "zhit_hit", m_nphlisthit, m_zhit_hit );
 
      // status = m_nt1->addItem("dEdx_meas_hit", m_dEdx_meas_hit);
      status = m_nt1->addItem( "dEdx_meas", m_dEdx_meas );
      // status = m_nt1->addItem("dEdx_meas_esat", m_dEdx_meas_esat);
      // status = m_nt1->addItem("dEdx_meas_norun", m_dEdx_meas_norun);
 
      status = m_nt1->addItem( "type", m_parttype );
      status = m_nt1->addItem( "chidedx_e", m_chidedxe );
      status = m_nt1->addItem( "chidedx_mu", m_chidedxmu );
      status = m_nt1->addItem( "chidedx_pi", m_chidedxpi );
      status = m_nt1->addItem( "chidedx_k", m_chidedxk );
      status = m_nt1->addItem( "chidedx_p", m_chidedxp );
      status = m_nt1->addItem( "partid", 5, m_probpid );
      status = m_nt1->addItem( "expectid", 5, m_expectid );
      status = m_nt1->addItem( "sigmaid", 5, m_sigmaid );
    }
  }
 
  NTuplePtr nt2( ntupleSvc(), "FILE100/n102" );
  if ( nt2 ) m_nt2 = nt2;
  else
  {
    m_nt2 = ntupleSvc()->book( "FILE100/n102", CLID_RowWiseTuple, "dE/dx per hit" );
    if ( m_nt2 )
    {
      status = m_nt2->addItem( "charge", m_charge1 );
      status = m_nt2->addItem( "adc_raw", m_phraw );
      status = m_nt2->addItem( "exraw", m_exraw );
      status = m_nt2->addItem( "runNO", m_runNO1 );
      status = m_nt2->addItem( "evtNO", m_evtNO1 );
      status = m_nt2->addItem( "runFlag", m_runFlag1 );
      status = m_nt2->addItem( "wire", m_wire );
      status = m_nt2->addItem( "doca_in", m_doca_in );
      status = m_nt2->addItem( "doca_ex", m_doca_ex );
      status = m_nt2->addItem( "driftdist", m_driftdist );
      status = m_nt2->addItem( "eangle", m_eangle );
      status = m_nt2->addItem( "zhit", m_zhit );
      status = m_nt2->addItem( "costheta1", m_costheta1 );
      status = m_nt2->addItem( "path_rphi", m_pathL );
      status = m_nt2->addItem( "layer", m_layer );
      status = m_nt2->addItem( "ptrk1", m_ptrk1 );
      status = m_nt2->addItem( "ptrk_hit", m_ptrk_hit );
      status = m_nt2->addItem( "t01", m_t01 );
      status = m_nt2->addItem( "tdc_raw", m_tdc_raw );
      status = m_nt2->addItem( "driftT", m_driftT );
      status = m_nt2->addItem( "localwid", m_localwid );
      status = m_nt2->addItem( "trackId1", m_trackId1 );
    }
  }
}
 
void DedxCalibEvent::genNtuple() {
  MsgStream log( msgSvc(), name() );
  log << MSG::INFO << "DedxCalibEvent::genNtuple()" << endmsg;
 
  SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
  if ( !eventHeader )
  {
    log << MSG::INFO << "Could not find Event Header" << endmsg;
    return;
  }
  int eventNO = eventHeader->eventNumber();
  int runNO   = eventHeader->runNumber();
  // only save part of the events, in each cluster only m_count events saved, then jump with a
  // gap
  // std::cout << "eventNO " << eventNO << "  evt_count " << evt_count << std::endl;
  if ( eventNO < evt_threshold ) return;
  evt_count++;
  if ( evt_count == m_count )
  {
    evt_threshold = eventNO + m_gap;
    evt_count     = 0;
  }
  // std::cout << "evt_threshold " << evt_threshold << "  evt_count " << evt_count <<
  // std::endl;
 
  log << MSG::INFO << "now begin the event:  runNO= " << runNO << " eventNO= " << eventNO
      << endmsg;
 
  int runFlag = 0; // data type flag
  if ( runNO < RUN0 ) runFlag = 0;
  if ( runNO >= RUN1 && runNO < RUN2 ) runFlag = 1;
  if ( runNO >= RUN2 && runNO < RUN3 ) runFlag = 2;
  if ( runNO >= RUN3 && runNO < RUN4 ) runFlag = 3;
  if ( runNO >= RUN4 && runNO < RUN5 ) runFlag = 4; // jpsi
  if ( runNO >= RUN5 && runNO < RUN6 ) runFlag = 5; // psipp
  if ( runNO >= RUN6 && runNO < RUN7 ) runFlag = 6; // psi4040, psip, jpsi...
  if ( runNO >= RUN7 ) runFlag = 7;                 // psip
 
  double                     tes        = -999.0;
  int                        esTimeflag = -1;
  SmartDataPtr<RecEsTimeCol> aevtimeCol( eventSvc(), "/Event/Recon/RecEsTimeCol" );
  if ( ( !aevtimeCol ) || ( aevtimeCol->size() == 0 ) )
  {
    tes        = -9999.0;
    esTimeflag = -1;
  }
  else
  {
    RecEsTimeCol::iterator iter_evt = aevtimeCol->begin();
    for ( ; iter_evt != aevtimeCol->end(); iter_evt++ )
    {
      tes        = ( *iter_evt )->getTest();
      esTimeflag = ( *iter_evt )->getStat();
    }
  }
  if ( runFlag == 2 )
  {
    if ( tes > 1000 ) return;
  }
  else if ( runFlag == 3 )
  {
    if ( tes > 700 ) return;
  }
  else
  {
    if ( tes > 1400 ) return;
  }
 
  SmartDataPtr<EvtRecEvent> evtRecEvent( eventSvc(), "/Event/EvtRec/EvtRecEvent" );
  if ( !evtRecEvent )
  {
    log << MSG::ERROR << "EvtRecEvent not found" << endmsg;
    return;
  }
 
  SmartDataPtr<EvtRecTrackCol> evtRecTrkCol( eventSvc(), "/Event/EvtRec/EvtRecTrackCol" );
  if ( !evtRecTrkCol )
  {
    log << MSG::ERROR << "EvtRecTrackCol" << endmsg;
    return;
  }
 
  Vint iGood;
  iGood.clear();
  int    nCharge = 0;
  double db = 0, dz = 0, pt0 = 0, charge0 = 0;
  for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    RecMdcDedx*         it    = ( *itTrk )->mdcDedx();
    if ( it == 0 ) continue;
    HepVector a;
    //  if((*itTrk)->isMdcKalTrackValid()) cout << "mdckaltrack is valid" << endl;
    //  else cout << "mdckaltrack is not valid" << endl;
 
    if ( ( *itTrk )->isMdcKalTrackValid() )
    {
      RecMdcKalTrack* trk = it->getMdcKalTrack();
      if ( ParticleFlag > -1 && ParticleFlag < 5 )
      {
        DstMdcKalTrack* dstTrk = trk;
        a                      = dstTrk->getZHelix( ParticleFlag );
      }
      else a = trk->getZHelix();
    }
    else if ( ( *itTrk )->isMdcTrackValid() )
    {
      RecMdcTrack* trk = it->getMdcTrack();
      a                = trk->helix();
    }
    else continue;
    db      = a[0];
    dz      = a[3];
    pt0     = fabs( 1.0 / a[2] );
    charge0 = ( a[2] > 0 ) ? 1 : -1;
 
    // cout<<"db: "<<db<<" dz: "<<dz<<" pt0: "<<pt0<<" charge0: "<<charge0<<endl;
    if ( fabs( dz ) >= VZ0CUT ) continue;
    if ( db >= VR0CUT ) continue;
    if ( pt0 >= PT0HighCut ) continue;
    if ( pt0 <= PT0LowCut ) continue;
    iGood.push_back( it->trackId() );
    nCharge += charge0;
  }
 
  double poca_x = 0, poca_y = 0, poca_z = 0;
  float  sintheta = 0, costheta = 0, ptrk = 0, ptrk_t = 0, charge = 0, trackId = 0;
  int    Nhit = 0, usedhit = 0, Nhits = 0, Nphlisthit = 0;
  double dEdx_meas_hit = 0, dEdx_meas = 0, dEdx_meas_esat = 0, dEdx_meas_norun = 0;
  double dEdx_hit[100] = { 0 }, pathlength_hit[100] = { 0 }, wid_hit[100] = { 0 },
         layid_hit[100] = { 0 }, dd_in_hit[100] = { 0 }, eangle_hit[100] = { 0 },
         zhit_hit[100]  = { 0 };
  int        trk_recalg = -1;
  Identifier mdcid;
 
  for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
  {
    EvtRecTrackIterator itTrk = evtRecTrkCol->begin() + i;
    RecMdcDedx*         it    = ( *itTrk )->mdcDedx();
    if ( it == 0 ) continue;
    bool flag = false;
    for ( unsigned int i = 0; i < iGood.size(); i++ )
    {
      if ( it->trackId() == iGood[i] ) flag = true;
    }
    if ( flag == false ) continue;
 
    HepVector    a;
    HepSymMatrix tkErrM;
    if ( ( *itTrk )->isMdcKalTrackValid() )
    {
      poca_x = it->getMdcKalTrack()->x(); // get poca, default pid is pion; change pid using
                                          // setPidType();
      poca_y = it->getMdcKalTrack()->y();
      poca_z = it->getMdcKalTrack()->z();
 
      RecMdcKalTrack* trk = it->getMdcKalTrack();
      if ( ParticleFlag > -1 && ParticleFlag < 5 )
      {
        DstMdcKalTrack* dstTrk = trk;
        a                      = dstTrk->getFHelix( ParticleFlag );
        tkErrM                 = dstTrk->getFError( ParticleFlag );
      }
      else
      {
        a      = trk->getFHelix();
        tkErrM = trk->getFError();
      }
    }
    else if ( ( *itTrk )->isMdcTrackValid() )
    {
      poca_x = it->getMdcTrack()->x();
      poca_y = it->getMdcTrack()->y();
      poca_z = it->getMdcTrack()->z();
 
      RecMdcTrack* trk = it->getMdcTrack();
      a                = trk->helix();
      tkErrM           = trk->err();
    }
    else continue;
 
    sintheta = sin( M_PI_2 - atan( a[4] ) );
    costheta = cos( M_PI_2 - atan( a[4] ) );
    ptrk_t   = 1.0 / fabs( a[2] );
    ptrk     = ptrk_t * sqrt( 1 + a[4] * a[4] );
    charge   = ( a[2] > 0 ) ? 1 : -1;
 
    Nhit = it->numTotalHits();               // total hits on track used as sample;
    Nhits = ( it->getVecDedxHits() ).size(); // dedx hits on this track, they are put in phlist
                                             // if layid>3
    usedhit = it->numGoodHits(); // hits after truncting phlist and used in cal dE/dx value;
    trk_recalg = it->status();
    trackId    = it->trackId();
 
    if ( m_eventType == "isBhabha" )
    {
      if ( runFlag == 3 && ( ptrk > 1.88 || ptrk < 1.80 ) ) continue;
      if ( runFlag == 4 && ( ptrk > 1.72 || ptrk < 1.45 ) ) continue;
      if ( runFlag == 5 && ( ptrk > 2.00 || ptrk < 1.70 ) ) continue;
      if ( runFlag == 6 && ( ptrk > 1.90 || ptrk < 1.00 ) ) continue;
      if ( runFlag == 7 && ( ptrk > 1.90 || ptrk < 0.90 ) ) continue;
      if ( abs( costheta ) > 0.9 ) continue;
 
      if ( Nhit < 20 ) continue;
      if ( usedhit < 6 ) continue;
    }
 
    int    layid = 0, localwid = 0, w0id = 0, wid = 0, lr = 0;
    double p_hit = 0, adc_raw = 0, tdc_raw = 0, zhit = 0, driftd = 0, driftD = 0, driftT = 0,
           dd_in = 0, dd_ex = 0, eangle = 0;
    double ph = 0, pathlength = 0, rphi_path = 0;
    long   k = 0;
 
    DedxHitRefVec           gothits   = it->getVecDedxHits();
    DedxHitRefVec::iterator it_gothit = gothits.begin();
    for ( ; it_gothit != gothits.end(); it_gothit++ )
    {
      if ( ( *it_gothit )->isMdcHitValid() )
      {
        RecMdcHit* itor = ( *it_gothit )->getMdcHit();
        mdcid           = itor->getMdcId();
        layid           = MdcID::layer( mdcid );
        localwid        = MdcID::wire( mdcid );
        w0id            = geosvc->Layer( layid )->Wirst();
        wid             = w0id + localwid;
        adc_raw         = itor->getAdc();
        tdc_raw         = itor->getTdc();
        zhit            = itor->getZhit();
 
        lr = itor->getFlagLR();
        if ( lr == 2 ) cout << "lr = " << lr << endl;
        if ( lr == 0 || lr == 2 ) driftD = itor->getDriftDistLeft();
        else driftD = itor->getDriftDistRight();
        driftd = abs( driftD );
        dd_in  = itor->getDoca();
        dd_ex  = itor->getDoca();
        if ( lr == 0 || lr == 2 )
        {
          dd_in = -abs( dd_in );
          dd_ex = -abs( dd_ex );
        }
        if ( lr == 1 )
        {
          dd_in = abs( dd_in );
          dd_ex = abs( dd_ex );
        }
        driftT = itor->getDriftT();
        eangle = itor->getEntra();
        eangle = eangle / M_PI;
      }
      else if ( ( *it_gothit )->isMdcKalHelixSegValid() )
      {
        RecMdcKalHelixSeg* itor     = ( *it_gothit )->getMdcKalHelixSeg();
        HepVector          a_hit_in = itor->getHelixIncl();
        p_hit = 1.0 / fabs( a_hit_in( 3 ) ) * sqrt( 1 + a_hit_in( 5 ) * a_hit_in( 5 ) );
 
        mdcid    = itor->getMdcId();
        layid    = MdcID::layer( mdcid );
        localwid = MdcID::wire( mdcid );
        w0id     = geosvc->Layer( layid )->Wirst();
        wid      = w0id + localwid;
        adc_raw  = itor->getAdc();
        tdc_raw  = itor->getTdc();
        zhit     = itor->getZhit();
 
        lr = itor->getFlagLR();
        if ( lr == 2 ) cout << "lr = " << lr << endl;
        driftD = itor->getDD();
        driftd = abs( driftD );
        driftT = itor->getDT();
        dd_in  = itor->getDocaIncl(); // getDocaIncl() include fit unused hit
        dd_ex  = itor->getDocaExcl(); // getDocaExcl() exclude fit unused hit
        if ( lr == -1 || lr == 2 )
        {
          dd_in = dd_in;
          dd_ex = dd_ex;
        }
        else if ( lr == 1 )
        {
          dd_in = -dd_in;
          dd_ex = -dd_ex;
        }
        eangle = itor->getEntra();
        eangle = eangle / M_PI;
      }
      else continue;
 
      pathlength = ( *it_gothit )->getPathLength();
      rphi_path  = pathlength * sintheta;
      ph         = ( *it_gothit )->getDedx();
      if ( layid > 3 )
      {
        dEdx_hit[k]       = adc_raw;
        pathlength_hit[k] = pathlength;
        wid_hit[k]        = wid;
        layid_hit[k]      = layid;
        dd_in_hit[k]      = dd_in;
        eangle_hit[k]     = eangle;
        zhit_hit[k]       = zhit;
 
        k++;
      }
 
      // cout<<"begin to Fill Ntuple n102!!!!!!!!!"<<endl;
      m_charge1   = charge;
      m_t01       = tes;
      m_driftT    = driftT;
      m_tdc_raw   = tdc_raw;
      m_phraw     = adc_raw;
      m_exraw     = ph;
      m_localwid  = localwid;
      m_wire      = wid;
      m_runNO1    = runNO;
      m_evtNO1    = eventNO;
      m_runFlag1  = runFlag;
      m_doca_in   = dd_in;
      m_doca_ex   = dd_ex;
      m_driftdist = driftD;
      m_eangle    = eangle;
      m_zhit      = zhit;
      m_costheta1 = costheta;
      m_pathL     = pathlength;
      m_layer     = layid;
      m_ptrk1     = ptrk;
      m_ptrk_hit  = p_hit;
      m_trackId1  = trackId;
      StatusCode status;
      if ( cut_wire > 0 )
      {
        if ( cut_wire == m_wire ) status = m_nt2->write();
      }
      else status = m_nt2->write();
      if ( status.isFailure() ) log << MSG::ERROR << "Cannot fill Ntuple n102" << endmsg;
    }
 
    Nphlisthit      = k; // dedx hits on this track, exclude the first 3 layers
    dEdx_meas       = it->probPH();
    dEdx_meas_esat  = it->getDedxEsat();
    dEdx_meas_norun = it->getDedxNoRun();
    dEdx_meas_hit   = it->getDedxHit();
 
    // cout<<"begin to Fill Ntuple n103!!!!!!!"<<endl;
    m_poca_x   = poca_x;
    m_poca_y   = poca_y;
    m_poca_z   = poca_z;
    m_ptrk_t   = ptrk_t;
    m_ptrk     = ptrk;
    m_sintheta = sintheta;
    m_costheta = costheta;
    m_charge   = charge;
    m_runNO    = runNO;
    m_runFlag  = runFlag;
    m_evtNO    = eventNO;
    m_t0       = tes;
    m_trackId  = trackId;
    m_recalg   = trk_recalg;
 
    m_nhit       = Nhit;
    m_nhits      = Nhits;
    m_nphlisthit = Nphlisthit;
    m_usedhit    = usedhit;
    for ( int j = 0; j < Nphlisthit; j++ )
    {
      m_dEdx_hit[j]       = dEdx_hit[j];
      m_pathlength_hit[j] = pathlength_hit[j];
      m_wid_hit[j]        = wid_hit[j];
      m_layid_hit[j]      = layid_hit[j];
      m_dd_in_hit[j]      = dd_in_hit[j];
      m_eangle_hit[j]     = eangle_hit[j];
      m_zhit_hit[j]       = zhit_hit[j];
    }
 
    // m_dEdx_meas_hit = dEdx_meas_hit;
    m_dEdx_meas = dEdx_meas;
    // m_dEdx_meas_esat = dEdx_meas_esat;
    // m_dEdx_meas_norun = dEdx_meas_norun;
 
    m_parttype  = it->particleId();
    m_chidedxe  = it->chiE();
    m_chidedxmu = it->chiMu();
    m_chidedxpi = it->chiPi();
    m_chidedxk  = it->chiK();
    m_chidedxp  = it->chiP();
    for ( int i = 0; i < 5; i++ )
    {
      m_probpid[i]  = it->getPidProb( i );
      m_expectid[i] = it->getDedxExpect( i );
      m_sigmaid[i]  = it->getSigmaDedx( i );
    }
    StatusCode status;
    if ( cut_wire < 0 ) status = m_nt1->write();
    if ( status.isFailure() ) { log << MSG::ERROR << "Cannot fill Ntuple n103" << endmsg; }
  }
  // cout<<"track iteration ended!!!!!!!!!!"<<endl;
}