---

TrackNtupleMaker.cxx

Go to the documentation of this file.

//  ====================================================================
//  TrackNtupleMaker.cxx
//  --------------------------------------------------------------------
//
//  This is the implementation of Track Ntuple Maker which 
//  creates the combined Ntuple based upon entities in the TES
//  This code is based upon AtlfastSTLNtupAlg by E. Richter-Was 
//  
//
//  Authors   : H.T. Phillips, P. Clarke, P. Sherwood, 
//                 R. Steward, E. Richter-Was
//
//  ====================================================================
#include "AtlfastAlgs/TrackNtupleMaker.h"
#include "AtlfastEvent/CollectionDefs.h"
#include "AtlfastUtils/HeaderPrinter.h"

// Framework include files
#include "GaudiKernel/PropertyMgr.h"
#include "GaudiKernel/INTupleSvc.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/NTuple.h"
#include "GaudiKernel/ObjectList.h"
#include "GaudiKernel/SmartDataPtr.h" //need for ntuple core code

// CLHEP,HepMC
//#include "GeneratorObjects/McEventCollection.h"
//#include "HepMC/GenEvent.h"




#include <cmath> 

namespace Atlfast {
  TrackNtupleMaker::TrackNtupleMaker(const std::string& name, 
                                     ISvcLocator* pSvcLocator) 
    :Algorithm(name, pSvcLocator){
    
    m_trackLocation            = "/Event/AtlfastTracks";
    declareProperty( "TrackLocation", m_trackLocation );
  }
  
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode TrackNtupleMaker::initialize() {
    MsgStream log(msgSvc(), name());
    log << MSG::DEBUG << "Initialising SNM" << endreq;
    StatusCode status;
    
    m_nEvents = 0;   
    
    //
    log << MSG::DEBUG << "Ntuple call 52" << endreq;
    NTuple::Directory* col = 0;
    log << MSG::DEBUG << "About to open ntuple file" << endreq;
    NTupleFilePtr file2(ntupleSvc(), "/NTUPLES/FILE1");
    if ( file2 )    {
      log << MSG::DEBUG << "Opened ntuple file" << endreq;
      if ( !(col=ntupleSvc()->createDirectory("/NTUPLES/FILE1/Atlfast")) )   {
        log << MSG::ERROR << "Cannot create directory /NTUPLES/FILE1/Atlfast" 
            << endreq;  
        return StatusCode::FAILURE;
      }else{
        log << MSG::DEBUG << "Created ntuple directory" << endreq;
      }
    }
    else    {
      log << MSG::ERROR << "Cannot access file /NTUPLES/FILE1" << endreq;  
      return StatusCode::FAILURE;
    }
    // First: A column wise N tuple
    NTuplePtr nt1(ntupleSvc(), "/NTUPLES/FILE1/Atlfast/52");
    if ( !nt1 )    {    // Check if already booked
      nt1 = ntupleSvc()->book (col, 52, CLID_ColumnWiseTuple, "Tracks52");
      if ( nt1 )    {
        log << MSG::DEBUG << "booked ntuple " << endreq;
        // Add an index column
        //tracks
        status = nt1->addItem ("PTRACKS/NTRA"      , m_ntra, 0, 500 );
        status = nt1->addItem ("PTRACKS/KPTRA"     , m_ntra, m_kpTruth);
        status = nt1->addItem ("PTRACKS/KFTRA"     , m_ntra, m_kfTruth);
        /*        
        status = nt1->addItem ("PTRACKS/KPM1TRA"   , m_ntra, m_kpm1tra);
        status = nt1->addItem ("PTRACKS/KFM1TRA"   , m_ntra, m_kfm1tra);
        status = nt1->addItem ("PTRACKS/KPM2TRA"   , m_ntra, m_kpm2tra);
        status = nt1->addItem ("PTRACKS/KFM2TRA"   , m_ntra, m_kfm2tra);
        status = nt1->addItem ("PTRACKS/KPM3TRA"   , m_ntra, m_kpm3tra);
        status = nt1->addItem ("PTRACKS/KFM3TRA"   , m_ntra, m_kfm3tra);
        status = nt1->addItem ("PTRACKS/KPM4TRA"   , m_ntra, m_kpm4tra);
        status = nt1->addItem ("PTRACKS/KFM4TRA"   , m_ntra, m_kfm4tra);
        status = nt1->addItem ("PTRACKS/KPM5TRA"   , m_ntra, m_kpm5tra);
        status = nt1->addItem ("PTRACKS/KFM5TRA"   , m_ntra, m_kfm5tra);
        status = nt1->addItem ("PTRACKS/KPM6TRA"   , m_ntra, m_kpm6tra);
        status = nt1->addItem ("PTRACKS/KFM6TRA"   , m_ntra, m_kfm6tra);
        */
        status = nt1->addItem ("PTRACKS/D0TRA  "      , m_ntra, m_d0Truth);
        status = nt1->addItem ("PTRACKS/Z0TRA  "      , m_ntra, m_z0Truth);
        status = nt1->addItem ("PTRACKS/PHITRA "      , m_ntra, m_phiTruth);
        status = nt1->addItem ("PTRACKS/COTTRA "      , m_ntra, m_cotTruth);
        status = nt1->addItem ("PTRACKS/PTINVTRA "    , m_ntra, m_ptInvTruth);
        status = nt1->addItem ("PTRACKS/RADTRA "     , m_ntra, m_radiusTruth);
        status = nt1->addItem ("PTRACKS/CURVTRA "    , m_ntra, m_curvTruth);
        
        status = nt1->addItem ("PTRACKS/D0TRACRU  "   , m_ntra, m_d0Track);
        status = nt1->addItem ("PTRACKS/Z0TRACRU  "   , m_ntra, m_z0Track);
        status = nt1->addItem ("PTRACKS/PHITRACRU "   , m_ntra, m_phiTrack);
        status = nt1->addItem ("PTRACKS/COTTRACRU "   , m_ntra, m_cotTrack);
        status = nt1->addItem ("PTRACKS/PTINVTRACRU " , m_ntra, m_ptInvTrack);
        status = nt1->addItem ("PTRACKS/RADTRACRU "     , m_ntra, m_radiusTrack);
        status = nt1->addItem ("PTRACKS/CURVTRACRU "    , m_ntra, m_curvTrack);
        
        status = nt1->addItem ("PTRACKS/CORR11TRA "   , m_ntra, m_corr11tra);
        status = nt1->addItem ("PTRACKS/CORR21TRA "   , m_ntra, m_corr12tra);
        status = nt1->addItem ("PTRACKS/CORR31TRA "   , m_ntra, m_corr13tra);
        status = nt1->addItem ("PTRACKS/CORR41TRA "   , m_ntra, m_corr14tra);
        status = nt1->addItem ("PTRACKS/CORR51TRA "   , m_ntra, m_corr15tra);
        status = nt1->addItem ("PTRACKS/CORR22TRA "   , m_ntra, m_corr22tra);
        status = nt1->addItem ("PTRACKS/CORR23TRA "   , m_ntra, m_corr23tra);
        status = nt1->addItem ("PTRACKS/CORR24TRA "   , m_ntra, m_corr24tra);
        status = nt1->addItem ("PTRACKS/CORR25TRA "   , m_ntra, m_corr25tra);
        status = nt1->addItem ("PTRACKS/CORR33TRA "   , m_ntra, m_corr33tra);
        status = nt1->addItem ("PTRACKS/CORR34TRA "   , m_ntra, m_corr34tra);
        status = nt1->addItem ("PTRACKS/CORR35TRA "   , m_ntra, m_corr35tra);
        status = nt1->addItem ("PTRACKS/CORR44TRA "   , m_ntra, m_corr44tra);
        status = nt1->addItem ("PTRACKS/CORR45TRA "   , m_ntra, m_corr45tra);
        status = nt1->addItem ("PTRACKS/CORR55TRA "   , m_ntra, m_corr55tra);
      }
      else    {   // did not manage to book the N tuple....
        return StatusCode::FAILURE;
      }
    }

    log<<MSG::DEBUG<<"returning with status"<< status<<endreq;
    m_tesIO = new TesIO();

    HeaderPrinter hp("Atlfast TrackNtupleMaker:", log);
    hp.add("Track Location           ", m_trackLocation);
    hp.print();

    return StatusCode::SUCCESS;
  }
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode TrackNtupleMaker::execute() {
    // This just makes the code below a bit easier to read (and type)
    MsgStream log(msgSvc(), name());
    
    // Counter of events processed
    if( m_nEvents % 50 == 0) {
      log << MSG::INFO << "    processing event number " 
          << m_nEvents << endreq;
    }     
    
    m_nEvents++;  
    
    StatusCode status;

    //-------------------
    // Retrieve Tracks 
    //-------------------

    std::vector<Track*> tr;
    if( !m_tesIO->copy<TrackCollection> (tr, m_trackLocation) ){
      log <<MSG::DEBUG << "Could not find Tracks"<<endreq;
    }
    log << MSG::DEBUG << "Processing tracks" <<tr.size()<< endreq;

    std::vector<Track*>::const_iterator itr  = tr.begin();
    
    m_ntra = 0;
    for (;itr != tr.end(); ++itr )   {
      log<<MSG::DEBUG<< " dumping track: "<<m_ntra<<endreq;
      log<<MSG::DEBUG<<(*itr)<<endreq<<endreq;
      //get track info for four vector:
      const HepMC::GenParticle* particle = (*itr)->truth();
      
      m_kpTruth[m_ntra]= 0;
      m_kfTruth[m_ntra]= particle->pdg_id();
      /*
      m_kpm1tra[m_ntra]= 0;
      m_kfm1tra[m_ntra]= 0;
      m_kpm2tra[m_ntra]= 0;
      m_kfm2tra[m_ntra]= 0;
      m_kpm3tra[m_ntra]= 0;
      m_kfm3tra[m_ntra]= 0;
      m_kpm4tra[m_ntra]= 0;
      m_kfm4tra[m_ntra]= 0;
      m_kpm5tra[m_ntra]= 0;
      m_kfm5tra[m_ntra]= 0;
      m_kpm6tra[m_ntra]= 0;
      m_kfm6tra[m_ntra]= 0;
      */
      //
      TrackTrajectory truthTraj = (*itr)->truthTrajectory();
      TrackTrajectory traj = (*itr)->trajectory();
      TrackParameters truthParam = truthTraj.parameters();
      TrackParameters smearParam = traj.parameters();
      
      m_d0Truth[m_ntra]     = truthParam.impactParameter();
      m_z0Truth[m_ntra]     = truthParam.zPerigee();
      m_phiTruth[m_ntra]    = truthParam.phi();
      m_cotTruth[m_ntra]    = truthParam.cotTheta();
      m_ptInvTruth[m_ntra]  = truthParam.invPtCharge();
      m_radiusTruth[m_ntra] = truthTraj.radius();
      m_curvTruth[m_ntra]   = truthTraj.curvature();
      //
      m_d0Track[m_ntra]    = smearParam.impactParameter();
      m_z0Track[m_ntra]    = smearParam.zPerigee();
      m_phiTrack[m_ntra]   = smearParam.phi();
      m_cotTrack[m_ntra]   = smearParam.cotTheta();
      m_ptInvTrack[m_ntra] = smearParam.invPtCharge();
      m_radiusTrack[m_ntra]     = traj.radius();
      m_curvTrack[m_ntra]       = traj.curvature();
      //
      m_corr11tra[m_ntra]= ( (*itr)->smearMatrix() )(1,1);
      m_corr12tra[m_ntra]= ( (*itr)->smearMatrix() )(1,2);
      m_corr13tra[m_ntra]= ( (*itr)->smearMatrix() )(1,3);
      m_corr14tra[m_ntra]= ( (*itr)->smearMatrix() )(1,4);
      m_corr15tra[m_ntra]= ( (*itr)->smearMatrix() )(1,5);
      m_corr22tra[m_ntra]= ( (*itr)->smearMatrix() )(2,2);
      m_corr23tra[m_ntra]= ( (*itr)->smearMatrix() )(2,3);
      m_corr24tra[m_ntra]= ( (*itr)->smearMatrix() )(2,4);
      m_corr25tra[m_ntra]= ( (*itr)->smearMatrix() )(2,5);
      m_corr33tra[m_ntra]= ( (*itr)->smearMatrix() )(3,3);
      m_corr34tra[m_ntra]= ( (*itr)->smearMatrix() )(3,4);
      m_corr35tra[m_ntra]= ( (*itr)->smearMatrix() )(3,5);
      m_corr44tra[m_ntra]= ( (*itr)->smearMatrix() )(4,4);
      m_corr45tra[m_ntra]= ( (*itr)->smearMatrix() )(4,5);
      m_corr55tra[m_ntra]= ( (*itr)->smearMatrix() )(5,5);

      if ( ++m_ntra == m_ntra->range().distance() ) break;
    }

    // Write out ntuple
    
    log << MSG::DEBUG << "Writing ntuple part 52" << endreq;
    status = ntupleSvc()->writeRecord("/NTUPLES/FILE1/Atlfast/52");
    if ( !status.isSuccess() ) log <<MSG::ERROR 
                                   <<"Cannot fill id 52"
                                   <<endreq;
    log << MSG::DEBUG << "Done all the ntuple writing" << endreq;
    return StatusCode::SUCCESS;
  }
  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
  
  StatusCode TrackNtupleMaker::finalize() {
    
    MsgStream log(msgSvc(), name());
    log << MSG::DEBUG 
        << "    Number of events processed      = " 
        << m_nEvents << endreq;
    
    return StatusCode::SUCCESS;
  }
  
  
} // end of namesapce bracket

---