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 "AtlfastAlgs/GlobalEventData.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/D0TRACRU",     m_ntra, m_d0Truth);
        status = nt1->addItem ("PTRACKS/Z0TRACRU",     m_ntra, m_z0Truth);
        status = nt1->addItem ("PTRACKS/PHITRACRU",    m_ntra, m_phiTruth);
        status = nt1->addItem ("PTRACKS/THETATRACRU",  m_ntra, m_thetaTruth);
        status = nt1->addItem ("PTRACKS/PINVTRACRU",   m_ntra, m_pInvTruth);
        status = nt1->addItem ("PTRACKS/RADTRACRU",    m_ntra, m_radiusTruth);
        status = nt1->addItem ("PTRACKS/CURVTRACRU",   m_ntra, m_curvTruth);
        
        status = nt1->addItem ("PTRACKS/D0TRA",        m_ntra, m_d0Track);
        status = nt1->addItem ("PTRACKS/Z0TRA",        m_ntra, m_z0Track);
        status = nt1->addItem ("PTRACKS/PHITRA",       m_ntra, m_phiTrack);
        status = nt1->addItem ("PTRACKS/THETATRA",     m_ntra, m_thetaTrack);
        status = nt1->addItem ("PTRACKS/PINVTRA",      m_ntra, m_pInvTrack);
        status = nt1->addItem ("PTRACKS/RADTRA",       m_ntra, m_radiusTrack);
        status = nt1->addItem ("PTRACKS/CURVTRA",      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;

    GlobalEventData* ged = GlobalEventData::Instance();
    //int randSeed = ged->randSeed() ;
    // load the location of the MC in StoreGate
    m_mcLocation       = ged -> mcLocation();

    m_tesIO = new TesIO(m_mcLocation, ged->justHardScatter());

    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_thetaTruth[m_ntra]  = truthParam.theta();
      m_pInvTruth[m_ntra]   = truthParam.invPCharge();
      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_thetaTrack[m_ntra]  = smearParam.theta();
      m_pInvTrack[m_ntra]   = smearParam.invPCharge();
      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

---