/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

Public Member Functions | Static Public Member Functions | Private Attributes
AdScaledEnergyTool Class Reference

#include <AdScaledEnergyTool.h>

Inheritance diagram for AdScaledEnergyTool:
Inheritance graph
[legend]
Collaboration diagram for AdScaledEnergyTool:
Collaboration graph
[legend]

List of all members.

Public Member Functions

 AdScaledEnergyTool (const std::string &type, const std::string &name, const IInterface *parent)
virtual ~AdScaledEnergyTool ()
virtual StatusCode reconstruct (const DayaBay::CalibReadout &readout, DayaBay::RecTrigger &recTrigger)
virtual StatusCode initialize ()
virtual StatusCode finalize ()

Static Public Member Functions

static const InterfaceID & interfaceID ()

Private Attributes

std::string m_cableSvcName
ICableSvcm_cableSvc
std::string m_calibStatsLocation
IPmtCalibSvcm_calibDataSvc
double m_intercept
double m_slope
bool m_IsMCIBD
std::string m_pmtGeomSvcName
IPmtGeomInfoSvcm_pmtGeomSvc
std::string m_calibDataSvcName
std::string m_ChannelQualitySvcName
IChannelQualitySvcm_chanqualSvc
std::string m_recDataSvcName
int m_error
IDetCalibSvcm_recDataSvc

Detailed Description

Definition at line 24 of file AdScaledEnergyTool.h.


Constructor & Destructor Documentation

AdScaledEnergyTool::AdScaledEnergyTool ( const std::string &  type,
const std::string &  name,
const IInterface *  parent 
)

Definition at line 26 of file AdScaledEnergyTool.cc.

  : GaudiTool(type,name,parent)
  , m_cableSvc(0)
  , m_calibDataSvc(0)
  , m_pmtGeomSvc(0)
{
    declareInterface< IReconTool >(this) ;   
    declareProperty("RecDataSvcName", m_recDataSvcName="DybDetCalibSvc",
                    "Name of calibration data service");
    declareProperty("PmtGeomSvcName", m_pmtGeomSvcName="PmtGeomInfoSvc",
                    "Name of PMT geometry data service");
    declareProperty("CableSvcName",m_cableSvcName="CableSvc",
                    "Name of service to map between detector, hardware, and electronic IDs");
    declareProperty("CalibDataSvcName", m_calibDataSvcName="DybPmtCalibSvc",
                    "Name of calibration data service");
    declareProperty("CalibStatsLocation",m_calibStatsLocation="/Event/Data/CalibStats",
                    "Location of the CalibStats");
    declareProperty("IsMCIBD", m_IsMCIBD = false, "is MC IBD?");
    declareProperty("ChannelQualitySvcName", m_ChannelQualitySvcName="DybChannelQualitySvc",
                    "Name of channel quality service");
}
AdScaledEnergyTool::~AdScaledEnergyTool ( ) [virtual]

Definition at line 50 of file AdScaledEnergyTool.cc.

{}

Member Function Documentation

StatusCode AdScaledEnergyTool::reconstruct ( const DayaBay::CalibReadout readout,
DayaBay::RecTrigger recTrigger 
) [virtual]

protect from nan values.

save reconstruct results

Implements IReconTool.

Definition at line 68 of file AdScaledEnergyTool.cc.

{
  double neutrino_e = 0;

  if( m_IsMCIBD ) {
    const DayaBay::GenHeader* gen_header = dynamic_cast<const DayaBay::GenHeader*>( ((readout.header())->findHeaders(51201))[0] );

    const HepMC::GenEvent* genevt = gen_header->event();
    for ( HepMC::GenEvent::particle_const_iterator p = genevt->particles_begin(); p != genevt->particles_end(); ++p ) {
      if( (*p)->pdg_id() == -12 ) {
        neutrino_e = ((*p)->momentum()).e();
      }
    }
  }

  if( !readout.detector().isAD() ){
    debug() << "Not an AD readout; ignoring detector " 
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kNotProcessed );
    return StatusCode::SUCCESS;
  }
  const CalibReadoutPmtCrate* pmtReadout 
    = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
  if(!pmtReadout){
    error() << "Incorrect type of readout crate for detector "
            << readout.detector().detName() << endreq;
    recTrigger.setEnergyStatus( ReconStatus::kBadReadout );
    return StatusCode::FAILURE;
  }
  double scaleConstant=160;

  // Context for this data
  int task = 0;
  ServiceMode svcMode(readout.header()->context(), task);
  const DayaBay::DetEnergyReconData* srcEn = m_recDataSvc->detEnergyReconData(svcMode);
  if(!srcEn){
                m_error++;
                if(m_error<10)warning() << "No energy scale found in DB for context " << readout.header()->context().AsString() << " and task " << task<< ". Will use a default value of " << scaleConstant << endreq;
                if(m_error==10) warning() << "No more errors about not finding energy scale in DB will be reported. A summary count of these errors will appear at the end." << endreq;
    }
                                
  else scaleConstant = srcEn->m_PeEvis;
  // Loop over hits and add up charge
  double qSum = 0;
  double qSumX = 0, qSumY = 0, qSumZ = 0;
  double qSquareSum = 0;
  int nHits = 0;
  CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
    chanEnd = pmtReadout->channelReadout().end();
  for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
      chanIter++){
    const CalibReadoutPmtChannel& channel = *chanIter;
    AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
    if( !pmtId.is8inch() ){ continue; }  // Ignore non 8-inch pmt channels

    double peakAdc=channel.earliestCharge(-1650,-1250);    
    qSum += peakAdc;
    const CLHEP::Hep3Vector& pmtPos 
      = m_pmtGeomSvc->get(pmtId.fullPackedData())->localPosition();
    qSumX += peakAdc*pmtPos[0];
    qSumY += peakAdc*pmtPos[1];
    qSumZ += peakAdc*pmtPos[2];
    qSquareSum += peakAdc*peakAdc;
    if(peakAdc>0)nHits++;
  }

  if(qSum<1 ){ //non-physics event
    CLHEP::HepLorentzVector position(0,0,0,0);
    recTrigger.setPosition( position );
    recTrigger.setRawEvis( 0 );
    recTrigger.setEnergy( 0 );
    recTrigger.setEnrec(0);
    recTrigger.setEprec(0);
    recTrigger.setPositionStatus( ReconStatus::kNoHits );
    recTrigger.setEnergyStatus( ReconStatus::kNoHits );
    return StatusCode::SUCCESS;
  } 
  qSumX/=qSum;
  qSumY/=qSum;
  qSumZ/=qSum;
  double qSumR = sqrt(qSumX*qSumX+qSumY*qSumY);
  double recR = qSumR*1.82 - 1.95e-4*qSumR*qSumR;
  if((readout.header()->context()).GetSimFlag()==2)
        recR = qSumR*1.60 - 1.36e-5*qSumR*qSumR;
  double recX = qSumX/qSumR*recR;
  double recY = qSumY/qSumR*recR;
  double recZ = (qSumZ - 1.579e-7*qSumZ*qSumZ*qSumZ)*(3.128-9.64e-4*qSumR);
  if((readout.header()->context()).GetSimFlag()==2)
        recZ = (qSumZ - 1.37e-7*qSumZ*qSumZ*qSumZ)*(2.86-9.02e-4*qSumR);
  if(recR>2500){
        recX = qSumX/qSumR*2500;
        recY = qSumY/qSumR*2500;
  }
  if(recZ>2500)recZ=2500;
  if(recZ<-2500)recZ=-2500;

  //active PMTs
  DayaBay::UserDataHeader* calibStatsHdr = 0;
  if(exist<UserDataHeader>(m_calibStatsLocation)){
    calibStatsHdr = get<UserDataHeader>(m_calibStatsLocation);
  }else{
    error() << "Cannot get CalibStats information" << endreq;
    error() << "AdScaled now require to run CalibStats before this!!!!" << endreq;
    return StatusCode::FAILURE;
  }
  double RingRatio[8] = {0.108,0.121,0.135,0.143,0.145,0.132,0.116,0.099};
  float nActivePMTs = calibStatsHdr->getFloat("nActivePMTs");
  if (nActivePMTs > 0){
        if(nActivePMTs==192)
                scaleConstant = scaleConstant * nActivePMTs/192.;
        else{
                int badPMTRing[192];
                for(int iit=0;iit<192;iit++) badPMTRing[iit]=0;
                const IChannelQuality * cq = m_chanqualSvc->channelQuality(svcMode);
                IChannelQuality::ChannelSet_t chans = cq->channels();
                int NN=0;
                for (int i = 0; i < chans.size(); i++){
                        DayaBay::FeeChannelId chid = chans[i];
                        if (!cq->good(chid)){
                                AdPmtSensor pmtId = m_cableSvc->adPmtSensor(chid,svcMode);
                                badPMTRing[NN] = pmtId.ring();
                                NN++;
                        }
                }
                for(int i=0;i<192;i++){
                        if(badPMTRing[i]!=0)
                                scaleConstant = scaleConstant * (1- RingRatio[badPMTRing[i]-1]/24);
                }
             }
  }
  
  double energy = 0;
  double rawEvis= 0;

  rawEvis = qSum / scaleConstant; //raw visible energy

  int siteId = readout.detector().site();
  int detId = readout.detector().detectorId();
  if     (siteId==1&&detId==1) detId=0;
  else if(siteId==1&&detId==2) detId=1;
  else if(siteId==2&&detId==1) detId=2;
  else if(siteId==4&&detId==1) detId=3;
  else if(siteId==4&&detId==2) detId=4;
  else if(siteId==4&&detId==3) detId=5;
  else if(siteId==4&&detId==4) detId=6;
  else if(siteId==2&&detId==2) detId=7;
  else detId=0;

  double uniCorr = (1+3.487e-08*recR*recR)*(1.0-1.524e-5*recZ-1.675e-8*recZ*recZ+1.982e-12*recZ*recZ*recZ);
  double linCorr = 8.05/CoToNGd[detId];
  if((readout.header()->context()).GetSimFlag()==2){
        uniCorr = (1+corrMC[0]*recR*recR)*(1+corrMC[1]*recZ+corrMC[2]*recZ*recZ+corrMC[3]*recZ*recZ*recZ);
        linCorr = 8.05/8.20;
  }

  double Enrec = linCorr*rawEvis/uniCorr;  //vertex correction

  double Eprec = 0;//positron non-linearity correction

  CLHEP::HepLorentzVector position(recX,recY,recZ,0);
  recTrigger.setPosition( position );
  recTrigger.setRawEvis( rawEvis );
  recTrigger.setEnergy( energy );
  recTrigger.setEnrec(Enrec);
  recTrigger.setEprec(Eprec);
  recTrigger.setPositionStatus( ReconStatus::kGood );
  // Set quality to standard deviation of hits
  recTrigger.setEnergyQuality(sqrt((qSquareSum - qSum*qSum/nHits)/(nHits-1)));
  recTrigger.setEnergyStatus( ReconStatus::kGood );

  if(m_IsMCIBD) {
    recTrigger.setDirectionQuality(neutrino_e);
  }

  return StatusCode::SUCCESS;
}
StatusCode AdScaledEnergyTool::initialize ( ) [virtual]

Definition at line 52 of file AdScaledEnergyTool.cc.

{
  // Get Cable Service
  m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
  m_calibDataSvc = svc<IPmtCalibSvc>(m_calibDataSvcName,true);
  m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName,true);
  m_recDataSvc = svc<IDetCalibSvc>(m_recDataSvcName,true);
  m_chanqualSvc = svc<IChannelQualitySvc>(m_ChannelQualitySvcName,true);
  return StatusCode::SUCCESS;
}
StatusCode AdScaledEnergyTool::finalize ( ) [virtual]

Definition at line 63 of file AdScaledEnergyTool.cc.

{
  return StatusCode::SUCCESS;
}

Member Data Documentation

std::string AdScaledEnergyTool::m_cableSvcName [private]

Definition at line 42 of file AdScaledEnergyTool.h.

Definition at line 45 of file AdScaledEnergyTool.h.

Definition at line 48 of file AdScaledEnergyTool.h.

Definition at line 50 of file AdScaledEnergyTool.h.

Definition at line 52 of file AdScaledEnergyTool.h.

double AdScaledEnergyTool::m_slope [private]

Definition at line 53 of file AdScaledEnergyTool.h.

Definition at line 54 of file AdScaledEnergyTool.h.

std::string AdScaledEnergyTool::m_pmtGeomSvcName [private]

Definition at line 57 of file AdScaledEnergyTool.h.

Definition at line 59 of file AdScaledEnergyTool.h.

Definition at line 61 of file AdScaledEnergyTool.h.

Definition at line 63 of file AdScaledEnergyTool.h.

Definition at line 65 of file AdScaledEnergyTool.h.

std::string AdScaledEnergyTool::m_recDataSvcName [private]

Definition at line 66 of file AdScaledEnergyTool.h.

Definition at line 67 of file AdScaledEnergyTool.h.

Definition at line 69 of file AdScaledEnergyTool.h.


The documentation for this class was generated from the following files:
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:12:37 for AdScaledEnergy by doxygen 1.7.4