MuonSimpleAlg.cc

Go to the documentation of this file.

#include "MuonSimpleAlg.h"

#include "Event/Readout.h"
#include "Event/ReadoutPmtCrate.h"
#include "Event/ReadoutPmtChannel.h"
#include "Event/CalibReadout.h"
#include "Event/CalibReadoutPmtCrate.h"
#include "Event/CalibReadoutPmtChannel.h"
#include "DataSvc/ICableSvc.h"
#include "Event/RawEventHeader.h"
#include "DetDesc/IGeometryInfo.h"
#include "DetDesc/DetectorElement.h"
#include "DetHelpers/IPmtGeomInfoSvc.h"
#include "CLHEP/Vector/LorentzVector.h"
#include <math.h>
#include <sstream>

#include "Event/RecHeader.h"
#include "Event/RecTrigger.h"
#include "Event/CalibReadoutHeader.h"
#include "Conventions/Reconstruction.h"
#include "Conventions/Detectors.h"
#include <vector>

using namespace std;
using namespace DayaBay;

MuonSimpleAlg::MuonSimpleAlg(const std::string& name, ISvcLocator* pSvcLocator)
  : DybAlgorithm<DayaBay::UserDataHeader>(name,pSvcLocator)
{
  declareProperty("Site",m_site="DayaBay","Site name (DayaBay/Far)");
  declareProperty("m_timeCutLow", m_timeCutLow=-1700., "Lower limit of the TDC cut to get rid of noise");
  declareProperty("m_timeCutHigh", m_timeCutHigh=-1400., "Higher limit of the TDC cut to get rid of noise");
}

MuonSimpleAlg::~MuonSimpleAlg(){}

StatusCode MuonSimpleAlg::initialize()
{
  StatusCode sc = this->GaudiAlgorithm::initialize();
  if(sc.isFailure()) return sc;

  sc = service("StaticCableSvc",m_cableSvc);
  if( sc.isFailure() ) {
    error() << "Can't get StaticCableSvc" << endreq;
    return sc;
  }

  sc = service("StaticCalibDataSvc",m_calibDataSvc);
  if( sc.isFailure() ) {
    error() << "Can't get StaticCalibDataSvc" << endreq;
    return sc;
  }

  sc = service("CoordSysSvc",mCss);
  if(sc.isFailure()) {
    error()<<"Can't get CoordSysSvc"<<endreq;
    return sc;
  }

  // Make use of svc<>
  IMessageSvc* msg = svc<IMessageSvc>("MessageSvc");
  MsgStream log(msg, name());
  log << MSG::DEBUG << "intializing ......" << endreq;
  
  // Get PmtGeoSvc
  IPmtGeomInfoSvc* pmtGeoSvc;
  sc = service("PmtGeomInfoSvc",pmtGeoSvc);
  if(sc.isFailure()) {
    error()<<"Can't get PmtGeomInfoSvc"<<endreq;
    return sc;
  }
  
  //choose Site
  //Site :: Site_t sitestr = Site :: FromString("DayaBay");
  Site :: Site_t sitestr = Site :: FromString(m_site.c_str());
  //choose DetectorId
  DetectorId::DetectorId_t idetid = DetectorId :: FromString("IWS");
  DetectorId::DetectorId_t odetid = DetectorId :: FromString("OWS");
  //for the inward_facing or outward_facing
  bool inward_facing = true;
  bool outward_facing = false;
  CLHEP::Hep3Vector global;
  Gaudi::XYZPoint pos;

  string searchDet="/dd/Structure/Pool/";
  if(m_site=="Far")     searchDet.append("far-ows");
  if(m_site=="DayaBay") searchDet.append("db-ows");
  m_de = getDet<IDetectorElement>(searchDet);

  int PMTcount = 0;
  //IWS
  for(int iwall_number=1;iwall_number<10;++iwall_number){
    for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,idetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}

      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);

      m_PmtId[pmtid] = PMTcount;
      PMTcount++;
    }
  }
  //Outer Pool OWS facing outward
  for(int iwall_number=1;iwall_number<9;++iwall_number){
    for(int iwall_spot=1;iwall_spot<17;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,outward_facing,sitestr,odetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}

      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);

      m_PmtId[pmtid] = PMTcount;
      PMTcount++;
    }      
  }
  //Outer Pool OWS facing inward
  for(int iwall_number=1;iwall_number<10;++iwall_number){
    for(int iwall_spot=1;iwall_spot<33;++iwall_spot){
      PoolPmtSensor poolPmtId(iwall_number,iwall_spot,inward_facing,sitestr,odetid);
      unsigned int pmtid = poolPmtId.fullPackedData();
      bool p = pmtGeoSvc->get(pmtid); if(!p){break;}

      global = pmtGeoSvc->get(pmtid)->globalPosition();
      Gaudi::XYZPoint global_point(global.x(), global.y(), global.z());
      m_PmtPropMap[pmtid].Position = m_de->geometry()->toLocal(global_point);

      m_PmtId[pmtid] = PMTcount;
      PMTcount++;
    }
  }

  return StatusCode::SUCCESS;
}

StatusCode MuonSimpleAlg::finalize()
{
  return this->GaudiAlgorithm::finalize();
}

StatusCode MuonSimpleAlg::execute()
{
  DayaBay::HeaderObject *muonTag = 0;

  if (exist<DayaBay::HeaderObject>(evtSvc(), "/Event/Tag/Muon/MuonAll")) {
    muonTag = get<DayaBay::HeaderObject>("/Event/Tag/Muon/MuonAll");
    //    calibHdrs = muonTag->findHeaders("/Event/CalibReadout/CalibReadoutHeader");
    std::vector<const DayaBay::IHeader*> calibHdrs = muonTag->findHeaders(CalibReadoutHeader::classID());

    //info() << "Calib header size:" << calibHdrs.size() << endreq;

    Gaudi::XYZPoint p1(0,0,0), p2(0,0,0);
    Gaudi::XYZVector vect;
    int ntrigger = 0;
    vector<const DayaBay::IHeader*>::iterator it;
    for(it=calibHdrs.begin(); it!=calibHdrs.end(); it++){
      const CalibReadoutHeader* calibHeader = dynamic_cast<const CalibReadoutHeader*>(*it);
      
      const CalibReadout* calibReadout = calibHeader->calibReadout();
      if(!calibReadout) {
        info() << "No CalibReadout retrieved this cycle." << endreq;
        continue;
      } 
      if(!calibReadout->detector().isWaterShield()) continue;

      ntrigger++;
      //reset values
      map<unsigned int/*PmtId*/,PmtProp>::iterator pmtit,idend=m_PmtPropMap.end();
      for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
        PmtProp& PoolPmt = pmtit->second;
        PoolPmt.Charge = 0;
        PoolPmt.FirstHit = 999999;
      }

      const CalibReadoutPmtCrate* pCalibPmtCrate = 0;    
      pCalibPmtCrate = dynamic_cast< const DayaBay::CalibReadoutPmtCrate* >(calibReadout);
      if(!pCalibPmtCrate) {
        error()<<"Failed to get CalibReadoutPmtCrate"<<endreq;
        return StatusCode::FAILURE;
      }
      
      Detector det = pCalibPmtCrate->detector();
      //ServiceMode svcMode( pCalibReadoutHdr->context(),0 );
    
      CalibReadoutPmtCrate::PmtChannelReadouts channels = pCalibPmtCrate->channelReadout();
      CalibReadoutPmtCrate::PmtChannelReadouts::iterator ci, ci_end = channels.end();
      for(ci = channels.begin(); ci!=ci_end; ci++) {
        unsigned int multi = ci->size();
        
        double charge = 0;
        for(unsigned int i=0; i<multi; i++){
          if(ci->time(i)>m_timeCutLow && ci->time(i)<m_timeCutHigh){
            if(charge<ci->charge(i)) charge = ci->charge(i);
            //if(hittime>calibChannel->time(i)) hittime = calibChannel->time(i);
          }
        } 

        //int sensor = ci->pmtSensorId().sensorId();
        //if(sensor==0)  continue;//skip fake PMT hits
        //PoolPmtSensor PmtId = PoolPmtSensor(sensor);
        //unsigned int pmtid = PmtId.fullPackedData();
        unsigned int pmtid = ci->pmtSensorId().fullPackedData();
        
        m_PmtPropMap[pmtid].Charge = charge;
        m_PmtPropMap[pmtid].FirstHit = ci->earliestTime();
      }
      
      //find maximum PMT
      double max = 0;
      double maxid = 0;
      for(pmtit=m_PmtPropMap.begin(); pmtit!=idend; pmtit++) {
        PmtProp& PoolPmt = pmtit->second;
        if(max<PoolPmt.Charge){
          max = PoolPmt.Charge;
          maxid = pmtit->first;
        }
      }
      
      //if(ntrigger==1) p1 = m_PmtPropMap[maxid].Position;
      //if(ntrigger==2) p2 = m_PmtPropMap[maxid].Position;
      if(calibReadout->detector().detectorId()==5) p1 = m_PmtPropMap[maxid].Position;
      if(calibReadout->detector().detectorId()==6) p2 = m_PmtPropMap[maxid].Position;
    }
    
    int PositionStatus = 0;
    float PositionX = 0;
    float PositionY = 0;
    float PositionZ = 0;
    int DirectionStatus = 0;
    float DirectionX = 0;
    float DirectionY = 0;
    float DirectionZ = 0;
    
    //when both IWS and OWS triggers exist
    if(p1.Mag2()>0 && p2.Mag2()>0){
      //info() << "simple muon reconstruction using IWS and OWS" << endreq;

      vect = p1-p2;
      PositionStatus = 1;
      PositionX = p1.X();
      PositionY = p1.Y();
      PositionZ = p1.Z();
      DirectionStatus = 1;    
      DirectionX = vect.X();
      DirectionY = vect.Y();
      DirectionZ = vect.Z();
      
    }
    
    // Create output user data header
    UserDataHeader* muonSimple = MakeHeaderObject();
    muonSimple->set("PositionStatus",PositionStatus);
    muonSimple->set("PositionX",PositionX);
    muonSimple->set("PositionY",PositionY);
    muonSimple->set("PositionZ",PositionZ);
    muonSimple->set("DirectionStatus",DirectionStatus);
    muonSimple->set("DirectionX",DirectionX);
    muonSimple->set("DirectionY",DirectionY);
    muonSimple->set("DirectionZ",DirectionZ);
    
  }
  return StatusCode::SUCCESS;
}