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

In This Package:

QsumEnergyTool.cc
Go to the documentation of this file.
00001 #include "QsumEnergyTool.h"
00002 
00003 #include "Event/CalibReadout.h"
00004 #include "Event/CalibReadoutPmtCrate.h"
00005 #include "Event/CalibReadoutPmtChannel.h"
00006 #include "Event/RecTrigger.h"
00007 #include "Event/UserDataHeader.h"
00008 
00009 #include "Conventions/Reconstruction.h"
00010 #include "Conventions/Detectors.h"
00011 #include <iostream>
00012 #include "DataSvc/ICableSvc.h"
00013 #include "DataSvc/ICalibDataSvc.h"
00014 #include "DybKernel/IDetCalibSvc.h"
00015 
00016 #include "CLHEP/Units/SystemOfUnits.h"
00017 #include "CLHEP/Vector/ThreeVector.h"
00018 
00019 #define NODB_ERRORS 15
00020 
00021 using namespace CLHEP;
00022 using namespace DayaBay;
00023 using namespace std;
00024 
00025 QsumEnergyTool::QsumEnergyTool(const std::string& type,
00026                                const std::string& name, 
00027                                const IInterface* parent)
00028   : GaudiTool(type,name,parent)
00029   , m_cableSvc(0)
00030 {
00031 
00032   declareInterface< IReconTool >(this) ;   
00033     declareProperty("CableSvcName",m_cableSvcName="CableSvc",
00034                     "Name of service to map between detector, hardware, and electronic IDs");
00035     declareProperty("RecDataSvcName", m_recDataSvcName="DybDetCalibSvc",
00036                     "Name of calibration data service");
00037     declareProperty("TaskForEnergyScale",m_task_es=1,
00038                     "Task number for DB energy scale lookup (0 for Co60; 1 for SN)");
00039     declareProperty("CalibStatsLocation",m_calibStatsLocation="/Event/Data/CalibStats",
00040                     "Location of the CalibStats");
00041     declareProperty("Charge",m_chargeLeaveName="NominalCharge",
00042                     "Name of the leave that hold the total charge in CalibStats");
00043 
00044     std::string Path = getenv ("DATASVCROOT");
00045     declareProperty("UniformityCorrFileName", m_uniformityCorrFileName = Path+"/share/uni-corr-p14a-v1.root",
00046                     "File name of the charge distribution template");
00047 
00048     declareProperty("UniformityCorrFileNameMC", m_uniformityCorrFileNameMC = Path+"/share/uni-corr-m14a-v1.root",
00049                     "File name of the charge distribution template for MC");
00050 
00051     // This root file include the new AD-by-AD non-uniformity correction.
00052     // Currently, this new non-uniformity correction works only for data.
00053     
00054 }
00055 
00056 QsumEnergyTool::~QsumEnergyTool(){}
00057 
00058 StatusCode QsumEnergyTool::initialize()
00059 {
00060   // Get Cable Service
00061   m_cableSvc = svc<ICableSvc>(m_cableSvcName,true);
00062   m_recDataSvc = svc<IDetCalibSvc>(m_recDataSvcName,true);
00063 
00064   // Initialize error counter
00065   nodb_count=0;
00066   
00067   info() << "Setting up non-uniformity correction using " << m_uniformityCorrFileName << endreq;
00068   m_corrector = new Corrector(m_uniformityCorrFileName.c_str());
00069   
00070   info() << "Setting up non-uniformity correction for MC using " << m_uniformityCorrFileNameMC << endreq;
00071   m_corrector_mc = new Corrector(m_uniformityCorrFileNameMC.c_str());
00072   
00073   return StatusCode::SUCCESS;
00074 
00075 }
00076 
00077 StatusCode QsumEnergyTool::finalize()
00078 {
00079 
00080   if(nodb_count>0){
00081     warning() << "There were " << nodb_count << " warnings about not being able to find the energy scale constant in the database" << endreq;
00082   }
00083 
00084   return StatusCode::SUCCESS;
00085 }
00086 
00087 StatusCode QsumEnergyTool::reconstruct(const CalibReadout& readout,
00088                                        RecTrigger& recTrigger)
00089 {
00090   if( !readout.detector().isAD() ){
00091     debug() << "Not an AD readout; ignoring detector " 
00092             << readout.detector().detName() << endreq;
00093     recTrigger.setEnergyStatus( ReconStatus::kNotProcessed );
00094     return StatusCode::SUCCESS;
00095   }
00096 
00097   const CalibReadoutPmtCrate* pmtReadout 
00098     = dynamic_cast<const CalibReadoutPmtCrate*>(&readout);
00099   if(!pmtReadout){
00100     error() << "Incorrect type of readout crate for detector "
00101             << readout.detector().detName() << endreq;
00102     recTrigger.setEnergyStatus( ReconStatus::kBadReadout );
00103     return StatusCode::FAILURE;
00104   }
00105 
00106 
00107   // Context for this data
00108   // Note: when looking up pmt status (right below) need task of zero. When looking up energy scale the task may change, according to the value of m_task_es (jpochoa)
00109   ServiceMode svcMode(readout.header()->context(), 0);
00110   ServiceMode svcMode_es(readout.header()->context(), m_task_es);
00111 
00112 
00113   // Loop over hits and add up charge
00114   double qSum = 0;
00115   double qSquareSum = 0;
00116   int nHits = 0;
00117 
00118   
00119   DayaBay::UserDataHeader* calibStatsHdr = 0;
00120   if(exist<UserDataHeader>(m_calibStatsLocation)){
00121     calibStatsHdr = get<UserDataHeader>(m_calibStatsLocation);
00122   }else{
00123     error() << "Cannot get CalibStats information" << endreq;
00124     error() << "QsumEnergyTool now require to run CalibStats before this!!!!" << endreq;
00125     return StatusCode::FAILURE;
00126 
00127     // followings are obsolete code
00128     
00129     CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
00130       chanEnd = pmtReadout->channelReadout().end();
00131     for(chanIter=pmtReadout->channelReadout().begin(); chanIter != chanEnd; 
00132         chanIter++){
00133       const CalibReadoutPmtChannel& channel = *chanIter;
00134       AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
00135       //AdPmtSensor pmtId = m_cableSvc->adPmtSensor(channel.channelId(), svcMode);
00136       if( pmtId.bogus() || pmtId.ring() < 1 ){ continue; }  // Ignore 2" pmts and non-pmt channels  (timing reference, noise... etc)
00137          
00138       //if(pmtId.ring()==1&&pmtId.column()==1) cout<<readout.detector().detName()<<' '<<pmtCalib->m_speHigh<<endl;
00139       //    if( pmtCalib->m_status != PmtCalibData::kGood ) { continue; } // this is obsolete
00140       
00141       double peakAdc = channel.earliestCharge(-1650,-1250);
00142       qSum += peakAdc;
00143       qSquareSum += peakAdc*peakAdc;
00144       if(peakAdc>0){
00145         nHits++;
00146       }
00147     }
00148 
00149   }
00150 
00151   
00152   qSum = calibStatsHdr->getFloat(m_chargeLeaveName);
00153   nHits = (int)calibStatsHdr->getInt("nHit");
00154   float nActivePMTs = calibStatsHdr->getFloat("nActivePMTs");
00155   
00156   
00157   if(nHits<1){
00158     recTrigger.setEnergyStatus( ReconStatus::kNoHits );
00159     return StatusCode::SUCCESS;
00160   }
00161 
00162   if(recTrigger.positionStatus()==ReconStatus::kNotConverged){
00163     return StatusCode::SUCCESS;
00164   }
00165 
00166   
00167   //Get position-dependent correction
00168   double posDepCorr=1.;
00169   if(recTrigger.positionStatus()==ReconStatus::kGood){
00170     const CLHEP::Hep3Vector pos = recTrigger.position();
00171     double position_unit = CLHEP::meter;
00172     double posx = pos.x()/position_unit;
00173     double posy = pos.y()/position_unit;
00174     double posz = pos.z()/position_unit;
00175     if ((readout.header()->context()).GetSimFlag()== SimFlag::kMC){
00176       int site = (readout.header()->context()).GetSite();
00177       int detid = (readout.header()->context()).GetDetId();
00178       posDepCorr = m_corrector_mc->get_corr(site,detid,posx,posy,posz);
00179     }else{
00180       // posDepCorr = uniformCorr(posx,posy,posz);//<-- uniformCorr takes reco pos in meters
00181       int site = (readout.header()->context()).GetSite();
00182       int detid = (readout.header()->context()).GetDetId();
00183       posDepCorr = m_corrector->get_corr(site,detid,posx,posy,posz);
00184 
00185     }
00186   }  
00187 
00188   // Perform scaling 
00189   //-->Get number from database
00190   double scaleConstant=170.; //default value to use for now (PEtoEvis)
00191   // Note: change task number when looking-up energy scale number 
00192   const DayaBay::DetEnergyReconData* srcEn = m_recDataSvc->detEnergyReconData(svcMode_es);
00193   if(!srcEn){
00194     nodb_count+=1;
00195     if(nodb_count < NODB_ERRORS){
00196       warning() << "No energy scale found in DB for context " << readout.header()->context().AsString() << " and task " << m_task_es << ". Will use a default value of " << scaleConstant << endreq;
00197     }
00198     if(nodb_count == NODB_ERRORS){
00199       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;
00200     }
00201   } else {
00202     scaleConstant = srcEn->m_PeEvis;
00203   }
00204 
00205   // Apply scale constant based on the number of live PMTs. Aug. 2012, Yasu.
00206   if (nActivePMTs > 0)
00207     scaleConstant = scaleConstant * nActivePMTs/192.;
00208   
00209   //-->do the math
00210   double energy = 0;
00211   double rawEvis = 0;
00212   if(scaleConstant>0) {
00213     rawEvis = (qSum *1./ scaleConstant);
00214     energy = rawEvis*posDepCorr;
00215   }else {
00216     error() << "scaleConstant found to be negative or zero! Something is seriously wrong, so will abort" << endreq;
00217     return StatusCode::FAILURE;
00218   }
00219   //--> storage
00220   recTrigger.setRawEvis( rawEvis );
00221   recTrigger.setEnergy( energy );  
00222   // (set quality to standard deviation of hits)
00223   recTrigger.setEnergyQuality(sqrt((qSquareSum - qSum*qSum/nHits)/(nHits-1)));
00224   recTrigger.setEnergyStatus( ReconStatus::kGood );
00225 
00226   return StatusCode::SUCCESS;
00227 }
00228 
00229 
00230 // Below are correction used for P13A, M13B and before, and no longer used.
00231 // P14A, M14A and later uses new correction function "Corrector()"
00232 
00233 
00234 
00235 // // Position-dependent correction 
00236 // // Important: reconstructed position must be provided in meters 
00237 // double QsumEnergyTool::uniformCorr(double reco_x, double reco_y, double reco_z){
00238   
00239 
00240 //   double reco_r = sqrt(pow(reco_x,2)+pow(reco_y,2));
00241 
00242 //   // Set maximum r and z position to avoid over-correcting
00243   
00244 //   double max_reco_r = 2.5;
00245 //   double max_reco_z = 2.5;
00246 //   double min_reco_z = -2.5;
00247   
00248 //   if (reco_r > max_reco_r) reco_r = max_reco_r;
00249 //   if (reco_z > max_reco_z) reco_z = max_reco_z;
00250 //   if (reco_z < min_reco_z) reco_z = min_reco_z;
00251   
00252 //   double corr = 8.05 / (7.84628 * (1 + 3.41294e-02*reco_r*reco_r) * (1 - 1.21750e-02*reco_z  - 1.64275e-02*reco_z*reco_z + 7.33006e-04*pow(reco_z,3)));
00253 
00254   
00255 //   return corr;
00256 // }
00257 
00258 // double QsumEnergyTool::uniformCorrMC(double reco_x, double reco_y, double reco_z){
00259   
00260 //   double reco_r = sqrt(pow(reco_x,2)+pow(reco_y,2));
00261 
00262 //   // Set maximum r and z position to avoid over-correcting
00263   
00264 //   double max_reco_r = 2.5;
00265 //   double max_reco_z = 2.5;
00266 //   double min_reco_z = -2.5;
00267   
00268 //   if (reco_r > max_reco_r) reco_r = max_reco_r;
00269 //   if (reco_z > max_reco_z) reco_z = max_reco_z;
00270 //   if (reco_z < min_reco_z) reco_z = min_reco_z;
00271   
00272 //   double corr = 7.937 / (7.66412 * (1 + 0.0399542*pow(reco_r,2) + 0.00116836*pow(reco_r,4)) * (1 - 0.00994061*reco_z -  0.0100378*reco_z*reco_z + 0.000827075*pow(reco_z,3)));
00273 
00274 //   return corr;
00275 // }
00276 
00277 
00278 Corrector::Corrector(const char* fname){
00279   infile = new TFile(fname);
00280   char hname[100];
00281   for (int detidx=0; detidx<8; ++detidx){
00282     sprintf(hname, "resp_r2z_%d", detidx+1);
00283     resp_r2z[detidx] = (TH2D*)infile->Get(hname);
00284 
00285     sprintf(hname, "resp_phi_%d", detidx+1);
00286     resp_phi[detidx] = (TF1*)infile->Get(hname);
00287   }
00288 }
00289 
00290 Corrector::~Corrector(){
00291   if (infile != 0 && infile->IsOpen())
00292     infile->Close();
00293   if (infile != 0)
00294     delete infile;
00295 }
00296 
00297 int Corrector::get_ad(int site, int detid){
00298   switch(site) {
00299   case 1:
00300     return detid;
00301     break;
00302   case 2:
00303     if (detid == 1)
00304       return 3;
00305     else if (detid == 2)
00306       return 8;
00307     break;
00308   case 4:
00309     return detid + 3;
00310     break;
00311   default:
00312     return 0;
00313   }
00314   return 0;
00315 }
00316 
00317 // site: 1,2,4 for EH1,2,3 respectively
00318 // detid: 1,2,3,4
00319 // x, y, z vertex in meter
00320 double Corrector::get_corr(int site, int detid,
00321                            double x, double y, double z) {
00322   int detidx = get_ad(site, detid) - 1;
00323   double r2 = x*x + y*y;
00324   double phi = atan2(y, x);
00325 
00326   if (abs(z) >= 2)
00327     z = z > 0 ? 1.999 : -1.999;
00328   if (r2 >= 4)
00329     r2 = 3.999;
00330 
00331   return 1. / resp_r2z[detidx]->Interpolate(r2, z)
00332     / resp_phi[detidx]->Eval(phi);
00333 }
00334 
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:13:18 for QsumEnergy by doxygen 1.7.4