AnaModule.cxx

Go to the documentation of this file.

#include <iomanip>
#include <TFile.h>
#include <TTree.h>
 
#include <fun4all/Fun4AllReturnCodes.h>
#include <phool/getClass.h>
#include <interface_main/SQHitVector_v1.h>
#include <interface_main/SQTrackVector_v1.h>
#include <interface_main/SQDimuonVector_v1.h>
 
 
#include "AnaModule.h"
 
AnaModule::AnaModule(const std::string& name): SubsysReco(name), p_geomSvc(GeomSvc::instance())
{}
 
AnaModule::~AnaModule()
{}
 
int AnaModule::Init(PHCompositeNode* topNode)
{
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int AnaModule::InitRun(PHCompositeNode* topNode)
{
  int ret = GetNodes(topNode);
  if(ret != Fun4AllReturnCodes::EVENT_OK) return ret;
 
  eventID = 0;
  MakeTree();
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int AnaModule::process_event(PHCompositeNode* topNode)
{
  int nTracklets = trackletVec->size();
  for(int i = 0; i < nTracklets; ++i)
  {
    Tracklet* tracklet = trackletVec->at(i);
    nHits = tracklet->getNHits();
    chisq = tracklet->getChisq();
 
    //very loose cuts here
    if(nHits < 9) continue;
    if(chisq > 20.) continue;
 
    for(auto it = detectorIDs.begin(); it != detectorIDs.end(); ++it)
    {
      detectorID = *it;
      double z_exp = p_geomSvc->getPlanePosition(detectorID);
      x_exp = tracklet->getExpPositionX(z_exp);
      y_exp = tracklet->getExpPositionY(z_exp);
      if(!p_geomSvc->isInPlane(detectorID, x_exp, y_exp)) continue;
 
      elementID_exp = p_geomSvc->getExpElementID(detectorID, tracklet->getExpPositionW(detectorID));
      if(elementID_exp < 1 || elementID_exp > p_geomSvc->getPlaneNElements(detectorID)) continue;
 
      SQHit* hit = findHit(detectorID, elementID_exp);
      elementID_closest = hit == nullptr ? -1 : hit->get_element_id();
 
      saveTree->Fill();
    }
  }
 
  ++eventID;
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int AnaModule::End(PHCompositeNode* topNode)
{
  saveFile->cd();
  saveTree->Write();
  saveFile->Close();
 
  return Fun4AllReturnCodes::EVENT_OK;
}
 
int AnaModule::GetNodes(PHCompositeNode* topNode)
{
  hitVector   = findNode::getClass<SQHitVector>(topNode, "SQHitVector");
  trackletVec = findNode::getClass<TrackletVector>(topNode, "TrackletVector");
  if(!hitVector || !trackletVec)
  {
    return Fun4AllReturnCodes::ABORTRUN;
  }
 
  return Fun4AllReturnCodes::EVENT_OK;
}
 
void AnaModule::MakeTree()
{
  saveFile = new TFile(saveName, "RECREATE");
 
  saveTree = new TTree("save", "Efficiency tree Created by AnaModule");
  saveTree->Branch("eventID", &eventID, "eventID/I");
  saveTree->Branch("detectorID", &detectorID, "detectorID/I");
  saveTree->Branch("elementID_exp", &elementID_exp, "elementID_exp/I");
  saveTree->Branch("elementID_closest", &elementID_closest, "elementID_closest/I");
  saveTree->Branch("x_exp", &x_exp, "x_exp/D");
  saveTree->Branch("y_exp", &y_exp, "y_exp/D");
  saveTree->Branch("nHits", &nHits, "nHits/I");
  saveTree->Branch("chisq", &chisq, "chisq/D");
}
 
void AnaModule::registerDetector(TString name)
{
  detectorIDs.insert(p_geomSvc->getDetectorID(name.Data()));
}
 
SQHit* AnaModule::findHit(int detID, int eleID)
{
  int delta = 999;
  SQHit* hit = nullptr;
  for(SQHitVector::Iter it = hitVector->begin(); it != hitVector->end(); ++it)
  {
    if((*it)->get_detector_id() != detID) continue;
    int delta_temp = abs(eleID - (*it)->get_element_id());
    if(delta > delta_temp)
    {
      delta = delta_temp;
      hit = (*it);
    }
  }
 
  return hit;
}