DstReformAlg.cxx

Go to the documentation of this file.

#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "EventModel/EventModel.h"
#include "EventModel/Event.h"
 
#include "EvtRecEvent/EvtRecEvent.h"
#include "EvtRecEvent/EvtRecTrack.h"
#include "EvtRecEvent/EvtRecDTag.h"
 
#include "DstReformAlg/DstReformAlg.h"
 
#include "GaudiKernel/SmartIF.h"
#include "GaudiKernel/IJobOptionsSvc.h"
#include "VertexFit/IVertexDbSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/ISvcLocator.h"
#include "CLHEP/Vector/ThreeVector.h"
#include "CLHEP/Vector/LorentzVector.h"
#include "CLHEP/Vector/TwoVector.h"
using CLHEP::Hep3Vector;
using CLHEP::Hep2Vector;
using CLHEP::HepLorentzVector;
#include "CLHEP/Geometry/Point3D.h"
#ifndef ENABLE_BACKWARDS_COMPATIBILITY
   typedef HepGeom::Point3D<double> HepPoint3D;
#endif
#include "CLHEP/Matrix/SymMatrix.h"
#include "VertexFit/Helix.h"
 
#include "RootCnvSvc/RootCnvSvc.h"
#include "RootEventData/TEvtRecObject.h"
#include "RootEventData/TDstEvent.h"
#include "RootEventData/TMdcTrack.h"
#include "RootEventData/TEvtHeader.h"
#include <TFile.h>
#include <TTree.h>
#include <map>
#include <vector>
#include <iostream>
#include <cmath>
using namespace CLHEP;
 
/////////////////////////////////////////////////////////////////////////////
 
DstReformAlg::DstReformAlg(const std::string& name, ISvcLocator* pSvcLocator) :
  Algorithm(name, pSvcLocator) {
  
  declareProperty("OutputDstFile", m_outputFile);
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode DstReformAlg::initialize(){
  MsgStream log(msgSvc(), name());
 
  log << MSG::INFO << "in initialize()" << endmsg;
  m_evtNum=0; 
  Gaudi::svcLocator()->service("TagFilterSvc", m_tagFilterSvc);
 
  m_dstDataType = m_tagFilterSvc->getDstDataType();
  std::cout<<"init DstReformAlg, m_dstDataType: "<< m_dstDataType <<std::endl;
 
  return StatusCode::SUCCESS;
}
 
StatusCode DstReformAlg::execute(){
  UInt_t tagData0 = m_tagFilterSvc->getTagData0();
  UInt_t tagData1 = m_tagFilterSvc->getTagData1();
  UInt_t tagData2 = m_tagFilterSvc->getTagData2();
  UInt_t tagData3 = m_tagFilterSvc->getTagData3();
  UInt_t tagData4 = m_tagFilterSvc->getTagData4();
 
  //tagData0 and tagData1 is for DST reform, tagData0 is first priority
 
  //for DTagSetAlg and LTagSetAlg, if tagData0=3, then reform these events with tagData1(nGoodCharged).
  int dstDataType = m_tagFilterSvc->getDstDataType();
  if( (dstDataType==4) && (tagData0==3) )
  {
    cmap[tagData0+tagData1].push_back(m_evtNum);
  }
  else
    cmap[tagData0].push_back(m_evtNum);
 
  //add tag information to each event
  if (m_dstDataType==1) {
     evtmap[m_evtNum].push_back(tagData0);
     evtmap[m_evtNum].push_back(tagData1);
     evtmap[m_evtNum].push_back(tagData2);
     evtmap[m_evtNum].push_back(tagData3);
     evtmap[m_evtNum].push_back(tagData4);
  } else if (m_dstDataType==2) {
     UInt_t tagData5 = m_tagFilterSvc->getTagData5();
     UInt_t tagData6 = m_tagFilterSvc->getTagData6();
     UInt_t tagData7 = m_tagFilterSvc->getTagData7();
     UInt_t tagData8 = m_tagFilterSvc->getTagData8();
     evtmap[m_evtNum].push_back(tagData1);
     evtmap[m_evtNum].push_back(tagData2);
     evtmap[m_evtNum].push_back(tagData3);
     evtmap[m_evtNum].push_back(tagData4);
     evtmap[m_evtNum].push_back(tagData5);
     evtmap[m_evtNum].push_back(tagData6);
     evtmap[m_evtNum].push_back(tagData7);
     evtmap[m_evtNum].push_back(tagData8);
  } else if ( (m_dstDataType==3) || (m_dstDataType==4) ) {
     UInt_t tagData5 = m_tagFilterSvc->getTagData5();
     UInt_t tagData6 = m_tagFilterSvc->getTagData6();
     UInt_t tagData7 = m_tagFilterSvc->getTagData7();
     UInt_t tagData8 = m_tagFilterSvc->getTagData8();
     UInt_t tagData9 = m_tagFilterSvc->getTagData9();
     evtmap[m_evtNum].push_back(tagData1);
     evtmap[m_evtNum].push_back(tagData2);
     evtmap[m_evtNum].push_back(tagData3);
     evtmap[m_evtNum].push_back(tagData4);
     evtmap[m_evtNum].push_back(tagData5);
     evtmap[m_evtNum].push_back(tagData6);
     evtmap[m_evtNum].push_back(tagData7);
     evtmap[m_evtNum].push_back(tagData8);
     evtmap[m_evtNum].push_back(tagData9);
 
  }
 
  m_evtNum++;
}
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * 
StatusCode DstReformAlg::finalize() {
  MsgStream log(msgSvc(), name());
 
  std::string fn;
  SmartIF<IJobOptionsSvc> iSvc(Gaudi::svcLocator()->service("JobOptionsSvc"));
  if ( iSvc.isValid() ) {
    const std::vector<const Property*>* ps = iSvc->getProperties("EventCnvSvc");
    std::vector<const Property*>::const_iterator it, end = ps->end();
    for ( it = ps->begin(); it != end; ++it ) {
      if ( (*it)->name() == "digiRootInputFile" ) {
        std::string flist=(*it)->toString();
        std::string::size_type p1 = flist.find('"') + 1;
        std::string::size_type p2 = flist.find('"', p1);
        fn = flist.substr(p1, p2-p1);
        std::cout<<"input dst file name: "<<fn<<std::endl;
      }
    }
  }
  TFile* f1 = TFile::Open(fn.c_str(),"READ");
  TFile* f2 = TFile::Open(m_outputFile.c_str(), "RECREATE");
  TTree* t1 = (TTree*)f1->Get("Event");
 
  log << MSG::INFO << "in finalize()" << endmsg;
  TTree* jt1 = (TTree*)(f1->Get("JobInfoTree"));                            
  TTree* jt2 = jt1->CloneTree();                                            
  TTree* t0 = new TTree("Metadata", "");
  Int_t  mode = -1,  begin = 0, end = 0;
  t0->Branch("mode",  &mode,  "mode/I");
  t0->Branch("begin", &begin, "begin/I");
  t0->Branch("end",   &end,   "end/I");
 
  for ( std::map<int, std::vector<long int> >::iterator it = cmap.begin(); it != cmap.end(); ++it ) {
      mode = it->first;
      begin = end;
      end += (it->second).size();
      std::cout<<"mode: "<<mode<< " begin: "<<begin<< " end: "<<end<<std::endl;
      t0->Fill();
  }
                                                                             
  TTree* t2 = t1->CloneTree(0);
  Int_t nEvents = t1->GetEntries();
  std::cout<<"Total nEvents: "<<nEvents<<std::endl;
  Int_t origEntry=0, currEntry=0;
  UInt_t tagData0, tagData1, tagData2, tagData3, tagData4;
  UInt_t tagData5, tagData6, tagData7, tagData8;
  TTree* t3 = new TTree("Entries", "");
  t3->Branch("origEntry", &origEntry, "origEntry/I");
  t3->Branch("currEntry", &currEntry, "currEntry/I");
  t3->Branch("tagData0", &tagData0, "tagData0/i");
  t3->Branch("tagData1", &tagData1, "tagData1/i");
  t3->Branch("tagData2", &tagData2, "tagData2/i");
  t3->Branch("tagData3", &tagData3, "tagData3/i");
  t3->Branch("tagData4", &tagData4, "tagData4/i");
  if (m_dstDataType==2) {
     t3->Branch("tagData5", &tagData5, "tagData5/i");
     t3->Branch("tagData6", &tagData6, "tagData6/i");
     t3->Branch("tagData7", &tagData7, "tagData7/i");
  }
  if ( (m_dstDataType==3) || (m_dstDataType==4) ) {
     t3->Branch("tagData5", &tagData5, "tagData5/i");
     t3->Branch("tagData6", &tagData6, "tagData6/i");
     t3->Branch("tagData7", &tagData7, "tagData7/i");
     t3->Branch("tagData8", &tagData8, "tagData8/i");
  }
 
  for ( std::map<int, std::vector<long int> >::iterator it = cmap.begin(); it != cmap.end(); ++it ) {
     int nNt = it->second.size();
     std::cout << "Writing " << nNt << "\t events with " << it->first << " TagData0..." << std::endl;
     for ( int i = 0; i < nNt; ++i ) {
        origEntry=(it->second)[i];
        t1->GetEntry(origEntry);
        tagData0=(evtmap[origEntry])[0];
        tagData1=(evtmap[origEntry])[1];
        tagData2=(evtmap[origEntry])[2];
        tagData3=(evtmap[origEntry])[3];
        tagData4=(evtmap[origEntry])[4];
        if (m_dstDataType==2) {
           tagData5=(evtmap[origEntry])[5];
           tagData6=(evtmap[origEntry])[6];
           tagData7=(evtmap[origEntry])[7];
        }
        if ( (m_dstDataType==3) || (m_dstDataType==4) ) {
           tagData5=(evtmap[origEntry])[5];
           tagData6=(evtmap[origEntry])[6];
           tagData7=(evtmap[origEntry])[7];
           tagData8=(evtmap[origEntry])[8];
        }
        t2->Fill();
        t3->Fill();
        currEntry++;
     }
     (it->second).clear();
  }
                                                                             
  f2->Write();
  f2->Close();
  return StatusCode::SUCCESS;
}