In This Package:

ValiNTuple.cc

Go to the documentation of this file.

#include "ValiNTuple.h"

#include "Event/SimHeader.h"
#include "Event/SimHitHeader.h"
#include "Event/SimHitCollection.h"
#include "Event/SimPmtHit.h"
#include "Event/SimHit.h"

#include "Event/SimParticleHistory.h"

#include "Event/SimTrack.h"
#include "Event/SimVertex.h"

#include "Conventions/Detectors.h"
#include "Conventions/Site.h"
#include "Conventions/DetectorId.h"

#include "GaudiKernel/ITHistSvc.h"

#include "TTree.h"

#include <map>
#include <string>
using namespace std;

typedef map<short int, pair<int,string> > detector_bins_t;
static detector_bins_t detector_bins;

typedef struct{
  double hitTime;
  int sensDetId;
  float weight;
  int ancestorPdg;
  double PosX;
  double PosY;
  double PosZ;  
  double wavelt;
} hitStruct;

ValiNTuple::ValiNTuple(const std::string& name, ISvcLocator* pSvcLocator)
    : GaudiAlgorithm(name,pSvcLocator)
    , m_hsvc(0)
    , m_outTree(0)
{
    declareProperty("Location",m_location=DayaBay::SimHeaderLocation::Default,
                    "Location in the TES to get HepMCEvents");
    declareProperty("FilePath",m_filepath="/file1/sim/",
                    "File path with with to register TTree .");

    if (!detector_bins.size()) {
        DayaBay::Detector dets[] = {
            DayaBay::Detector(Site::kDayaBay,DetectorId::kAD1),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kAD2),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kOWS),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kIWS),
            DayaBay::Detector(Site::kDayaBay,DetectorId::kRPC),
            DayaBay::Detector(Site::kLingAo,DetectorId::kAD1),
            DayaBay::Detector(Site::kLingAo,DetectorId::kAD2),
            DayaBay::Detector(Site::kLingAo,DetectorId::kOWS),
            DayaBay::Detector(Site::kLingAo,DetectorId::kIWS),
            DayaBay::Detector(Site::kLingAo,DetectorId::kRPC),
            DayaBay::Detector(Site::kFar,DetectorId::kAD1),
            DayaBay::Detector(Site::kFar,DetectorId::kAD2),
            DayaBay::Detector(Site::kFar,DetectorId::kAD3),
            DayaBay::Detector(Site::kFar,DetectorId::kAD4),
            DayaBay::Detector(Site::kFar,DetectorId::kOWS),
            DayaBay::Detector(Site::kFar,DetectorId::kIWS),
            DayaBay::Detector(Site::kFar,DetectorId::kRPC),
            DayaBay::Detector(),
        };

        detector_bins[0] = pair<int,string>(0,"unknown");
        for (int ind=0; dets[ind].site(); ++ind) {
            detector_bins[dets[ind].siteDetPackedData()] 
                = pair<int,string>(ind+1,dets[ind].detName());
        }
    }    

    
}
ValiNTuple::~ValiNTuple()
{
}

StatusCode ValiNTuple::initialize()
{
    this->GaudiAlgorithm::initialize();

    // Get the histogram service
    if ( service("THistSvc", m_hsvc).isFailure()) {
        error() << " No THistSvc available." << endreq;
        return StatusCode::FAILURE;
    } 

    m_outTree = new TTree("nOutTree","Output from DetSim");
    string outstring= m_filepath+"nOutTree";
    if (m_hsvc->regTree(outstring,m_outTree).isFailure()) {
        error() << "Could not register " << m_filepath+"nOutTree" << endreq;
        delete m_outTree; m_outTree = 0;
        return StatusCode::FAILURE;
    }

    return StatusCode::SUCCESS;
}

StatusCode ValiNTuple::execute()
{

  DayaBay::SimHeader* shead = get<DayaBay::SimHeader>(m_location);

  // Get hit info 
  const DayaBay::SimHitHeader* shith = shead->hits();
  const DayaBay::SimHitHeader::hc_map& hcmap = shith->hitCollection();
  DayaBay::SimHitHeader::hc_map::const_iterator it, done = hcmap.end();

  
  // Get Primary Track Info.
  const DayaBay::SimParticleHistory* shist = shead->particleHistory();
  const std::list<const DayaBay::SimTrack*>& priTrks = shist->primaryTracks();
  std::list<const DayaBay::SimTrack*>::const_iterator priIt, priEnd=priTrks.end();

  //Get UnObservable Statistics info
  const DayaBay::SimUnobservableStatisticsHeader* unobSt = shead->unobservableStatistics();
  const DayaBay::SimUnobservableStatisticsHeader::stat_map& statmap = unobSt->stats();
  DayaBay::SimUnobservableStatisticsHeader::stat_map::const_iterator st, stdone = statmap.end();
  
  
  //  if(statmap.size()){
  for (st=statmap.begin(); st != stdone; ++st) {
    
    info () << " stat first " << st->first << " stat second sum " << st->second.sum() <<endreq;
  }
  //  }
  
  
  int hitsize(0),hitInCont(0);
  double trkLengInOws(0),trkLengInIws(0);
  
  hitsize=hcmap.size();
  info () << " Hit size " << hitsize <<endreq;

  //  hitStruct* hitCol[hitsize];

  hitStruct singleHit;

  //**** hit information**
  m_outTree->Branch("hitTime",&singleHit.hitTime,"hitTime/D");
  m_outTree->Branch("hitDetID",&singleHit.sensDetId,"hitDetID/I");
  m_outTree->Branch("hitWt",&singleHit.weight,"hitWt/F");
  m_outTree->Branch("hitancPDG",&singleHit.ancestorPdg,"hitancPDG/I");
  m_outTree->Branch("hitpX",&singleHit.PosX,"hitpX/D");
  m_outTree->Branch("hitpY",&singleHit.PosY,"hitpY/D");
  m_outTree->Branch("hitpZ",&singleHit.PosZ,"hitpZ/D");
  m_outTree->Branch("hitWavl",&singleHit.wavelt,"hitWavl/D");

  //***  ***
  m_outTree->Branch("hitsize",&hitsize,"hitsize/I");
  m_outTree->Branch("hitInCont",&hitInCont,"hitInCont/I");
  m_outTree->Branch("trkLengInOws",&trkLengInOws,"trkLengInOws/D");
  m_outTree->Branch("trkLengInIws",&trkLengInIws,"trkLengInIws/D");


  info () << " 01010101010101010 " << endreq;

  //  m_outTree.Branch("hitCol", &hitCol,"hitCol[hitsize]/F");
  
  info () << " 000000000 " << endreq;



  // Iterate to get the hit informatioin.
  for (it=hcmap.begin(); it != done; ++it) {
    DayaBay::Detector det(it->first);
    pair<int,string> bin = detector_bins[det.siteDetPackedData()];
    
    info () << "Got hit collection from " << det.detName() 
            << " in bin " << bin.first << endreq;
    
    size_t siz = it->second->collection().size();
    info () << " #hits in this Det " << siz << endreq;
    
    //    const DayaBay::SimHitCollection::hit_container hitcont=it->second->collection();
    //    DayaBay::SimHitCollection::hit_container::const_iterator ithit;

    const std::vector<DayaBay::SimHit*>&  hitcont=it->second->collection();
    std::vector<DayaBay::SimHit*>::const_iterator ithit;

    hitInCont=siz;
    
    if(siz!=0){
      for(ithit=hitcont.begin(); ithit !=hitcont.end();++ithit) {
        singleHit.hitTime=(*ithit)->hitTime();
        //      info () << " Hit time "<< singleHit.hitTime << endreq;
        singleHit.sensDetId=(*ithit)->sensDetId();
        singleHit.weight=(*ithit)->weight();
        if((*ithit)->ancestor().track()){
          singleHit.ancestorPdg=(*ithit)->ancestor().track()->particle(); //particle PDG//generate error
        }
        singleHit.PosX=(*ithit)->localPos().x();
        singleHit.PosY=(*ithit)->localPos().y();
        singleHit.PosZ=(*ithit)->localPos().z();
        //      info () << " x: "<< singleHit.PosX <<singleHit.PosY << singleHit.PosZ << endreq;
        
        const DayaBay::SimPmtHit* pmthit = static_cast<const DayaBay::SimPmtHit*>(*ithit);
        
        //      info () << " 22222222 " << endreq;
        
        //      singleHit.wavelt=(*pmthit).wavelength();
        //      m_outTree->Fill();
      }
    }
    else m_outTree->Fill();
  }



  
  //  m_outTree->Write();
  return StatusCode::SUCCESS;
}

StatusCode ValiNTuple::finalize()
{

    return this->GaudiAlgorithm::finalize();
}

---