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

In This Package:

QMLFTool.cc
Go to the documentation of this file.
00001 #include "QMLFTool.h"
00002 
00003 #include "DetHelpers/IPmtGeomInfoSvc.h"
00004 #include "DetHelpers/ICoordSysSvc.h"
00005 #include "DataSvc/ICalibDataSvc.h"
00006 #include "Conventions/Detectors.h"
00007 #include "Conventions/Reconstruction.h"
00008 #include "Event/RecTrigger.h"
00009 #include "Event/CalibReadoutPmtCrate.h"
00010 
00011 #include "TMinuit.h"
00012 #include "TGraph.h"
00013 #include "TMath.h"
00014 
00015 #include "DetDesc/Material.h"
00016 #include "DetDesc/Surface.h"
00017 #include "DetDesc/IGeometryInfo.h"
00018 #include "DetDesc/TabulatedProperty.h"
00019 #include <fstream>
00020 #include <cstdlib>
00021 
00022 #define NREF_ORDER 3
00023 
00024 //int  tttttttt = 0;
00025 QMLFTool* QMLFTool::current_instance = 0;
00026 
00027 void QMLFTool::NCLL_FCN(int& npar, double* grad, double& fval, double* xval, int iflag)
00028 {
00029     current_instance->ncll_fcn(npar, grad, fval, xval, iflag);
00030 }
00031 
00032 QMLFTool::QMLFTool(const std::string& type,
00033                    const std::string& name,
00034                    const IInterface* parent)
00035   : GaudiTool(type, name, parent),
00036     m_opPara(1.25E+04, 0.98, 0.98, 9000, 3.14*103*103),
00037     m_geomPara(2010., -2010.)
00038 {
00039     declareInterface<IReconTool>(this);
00040 
00041     declareProperty("PmtGeomSvcName", m_pmtGeomSvcName = "PmtGeomInfoSvc", 
00042             "Name of Pmt Geometry Information Service");
00043     declareProperty("PmtCalibDataSvcName", m_pmtCalibDataSvcName = "CalibDataSvc", 
00044             "Name of Pmt Calib Data Information Service");
00045 
00046     // declare property for those optical parameters
00047     declareProperty("LLFMode", m_LLFMode = 1, "Likelihood function mode");
00048     declareProperty("addTimeLLF", m_addTimeLLF = false, "Add time likelihood or not");
00049     declareProperty("pmtChgCut", m_pmtChgCut = 12., "The charge cut for PMTs selection criteria");
00050     declareProperty("opLocation", m_opLocation="DetDesc", "Optical parameters location");
00051     declareProperty("geomLocation", m_geomLocation="DetDesc", "Optical parameters location");
00052     declareProperty("AbsLength", m_opPara.m_attl=-1., "Effective absorption length");
00053     declareProperty("TopReflectivity", m_opPara.m_topRef=-1., "Reflectivity of top reflector");
00054     declareProperty("BotReflectivity", m_opPara.m_botRef=-1., "Reflectivity of bottom reflector");
00055     declareProperty("TopRefZ", m_geomPara.m_topRefZ=0., "Z position of top reflector");
00056     declareProperty("BotRefZ", m_geomPara.m_botRefZ=0., "Z position of bottom reflector");
00057     declareProperty("m_timeCutLow", m_timeCutLow=-1700., "Lower limit of the TDC cut to get rid of noise");
00058     declareProperty("m_timeCutHigh", m_timeCutHigh=-1400., "Higher limit of the TDC cut to get rid of noise");
00059     declareProperty("InitSource", m_src = "/Event/Rec/AdScaled", "The initial value for fitting");
00060 }
00061 
00062 QMLFTool::~QMLFTool()
00063 {
00064 }
00065 
00066 StatusCode QMLFTool::initialize()
00067 {
00068     // Get Pmt Geometry Service
00069     m_pmtGeomSvc = svc<IPmtGeomInfoSvc>(m_pmtGeomSvcName, true);
00070 
00071     // Get Pmt Calib Data Service
00072     m_pmtCalibDataSvc = svc<ICalibDataSvc>(m_pmtCalibDataSvcName, true);
00073 
00074     // Get Parameters
00075     getOpPara();
00076     getGeomPara();
00077     getPmtEff();
00078 
00079     // Show Optical and Geom Parameters
00080     printPara();
00081 
00082     return StatusCode::SUCCESS;
00083 }
00084 
00085 StatusCode QMLFTool::finalize()
00086 {
00087     return StatusCode::SUCCESS;
00088 }
00089 
00090 StatusCode QMLFTool::reconstruct(const DayaBay::CalibReadout& calibReadout,
00091                                  DayaBay::RecTrigger& recTrigger)
00092 {
00093     //tttttttt = 0;
00094     if ( ! calibReadout.detector().isAD() ) {
00095         debug() << "Not an AD readout; ignoring detector " 
00096                 << calibReadout.detector().detName() << endreq;
00097         recTrigger.setPositionStatus( ReconStatus::kNotProcessed);
00098         recTrigger.setEnergyStatus( ReconStatus::kNotProcessed);
00099         return StatusCode::SUCCESS;
00100     }
00101 
00102     m_calibCrate = dynamic_cast<const DayaBay::CalibReadoutPmtCrate*>(&calibReadout);
00103     if ( !m_calibCrate ) {
00104         error() << "Incorrect type of readout crate for detector "
00105                 << calibReadout.detector().detName() << endreq;
00106         recTrigger.setPositionStatus(ReconStatus::kBadReadout);
00107         recTrigger.setEnergyStatus(ReconStatus::kBadReadout);
00108         return StatusCode::FAILURE;
00109     }
00110 
00111     // Initial the results from user specified source
00112     const DayaBay::RecHeader* srcRecHeader = get<DayaBay::RecHeader>(m_src);
00113     //if ( ! srcRecHeader ) {
00114     // // error handling
00115     //}
00116     const DayaBay::RecTrigger& srcTrigger = srcRecHeader->recTrigger();
00117     const CLHEP::HepLorentzVector& srcPos = srcTrigger.position();
00118     double xini = srcPos.x();
00119     double yini = srcPos.y();
00120     double zini = srcPos.z();
00121     double eini = srcTrigger.energy();
00122     //double tini = recTrigger.triggerTime().GetNanoSec();
00123 
00124     // Set this QMLFTool instance as the current_instance
00125     current_instance = this;
00126 
00127     // fill the local buffers
00128     Site::Site_t site = m_calibCrate->detector().site();
00129     DetectorId::DetectorId_t  detid = m_calibCrate->detector().detectorId();
00130     m_svcMode = ServiceMode(m_calibCrate->header()->context(), 0);
00131     if ( site != Site::kFar) {
00132         m_ADNo = (site>>1)*2 + detid - 1;
00133     }
00134 
00135     for ( int local_id = 0; local_id < 192; ++local_id ) {
00136         m_obsq[local_id] = 0.0;
00137 
00138         int ring = local_id/24 + 1;
00139         int col  = local_id%24 + 1;
00140         DayaBay::AdPmtSensor adPmtId(ring, col, site, detid);
00141         const DayaBay::PmtCalibData* pmtData =
00142             m_pmtCalibDataSvc->pmtCalibData(adPmtId, m_svcMode);
00143         // TODO: optimize : if (m_vPmtData[local_id] == pmtData) break;
00144         // m_vPmtData[4][192];
00145         m_vPmtData[local_id] = pmtData;
00146         if ( pmtData != 0 && pmtData->m_status == DayaBay::PmtCalibData::kGood ) {
00147             m_chgLLFUsedPMTs[local_id] = adPmtId.fullPackedData();
00148         }
00149         else {
00150             m_chgLLFUsedPMTs[local_id] = 0;
00151         }
00152     }
00153     DayaBay::CalibReadoutPmtCrate::PmtChannelReadouts::const_iterator chanIter, 
00154         chanEnd = m_calibCrate->channelReadout().end();
00155     for(chanIter=m_calibCrate->channelReadout().begin(); chanIter != chanEnd; ++chanIter) {
00156         const DayaBay::CalibReadoutPmtChannel& channel = *chanIter;
00157         const DayaBay::AdPmtSensor pmtId(channel.pmtSensorId().fullPackedData());
00158         if ( ! pmtId.is8inch() ) continue;
00159         const std::vector<double>& vTime = channel.time();
00160         const std::vector<double>& vChrg = channel.charge();
00161         double maxCharge = 0.0;
00162         int local_id = (pmtId.ring()-1)*24 + pmtId.column() - 1;
00163         std::vector<double>::const_iterator iT = vTime.begin(), iC = vChrg.begin(), iEnd = vChrg.end();
00164         while ( iC != iEnd ) {
00165             if ( *iT > m_timeCutLow && *iT < m_timeCutHigh ) {
00166                 if ( maxCharge < *iC ) maxCharge = *iC;
00167             }
00168             ++iC; ++iT;
00169         }
00170         m_obsq[local_id] = maxCharge;
00171     }
00172     // End of filling local buffers
00173 
00174     //
00175     // Set parameters for Minuit minimization
00176     static const int Npar = 4;
00177     int    ierflg = 0;
00178     double standin = 0.0;
00179     double maxCalls = 2000.;
00180     double tolerance = 1.0e-3;
00181     double arglist[2];
00182 
00183     //Create TMinuit
00184     TMinuit * recMinuit = new TMinuit(Npar);
00185     recMinuit->SetFCN(QMLFTool::NCLL_FCN);
00186     recMinuit->SetPrintLevel(-1);
00187 
00188     // Initial TMinuit Fitting Parameter
00189     recMinuit->mnparm(0, "xpos",   xini,  100.0, 0.0, 0.0, ierflg);
00190     recMinuit->mnparm(1, "ypos",   yini,  100.0, 0.0, 0.0, ierflg);
00191     recMinuit->mnparm(2, "zpos",   zini,  10.0, m_geomPara.m_botRefZ, m_geomPara.m_topRefZ, ierflg);
00192     recMinuit->mnparm(3, "energy", eini,   0.2, 0.001, 1000.0, ierflg);
00193 
00194     // No Warnings
00195     recMinuit->mnexcm("SET NOW", &standin, 0, ierflg);
00196 
00197     // Set error Definition
00198     // 1 for Chi square
00199     // 0.5 for negative log likelihood
00200     if(m_LLFMode==1)
00201         recMinuit->SetErrorDef(1);
00202     else 
00203         recMinuit->SetErrorDef(0.5);
00204 
00205     // Minimization strategy
00206     // 1 standard;  2 try to improve minimum (slower)
00207     arglist[0] = 2;
00208     recMinuit->mnexcm("SET STR", arglist, 1, ierflg); 
00209 
00210     //Do Minuit fit!  
00211     arglist[0] = maxCalls;
00212     arglist[1] = tolerance;
00213     recMinuit->mnexcm("MIGrad", arglist, 2, ierflg);
00214     //recMinuit->mnexcm("MINIMIZE", arglist, 2, ierflg);
00215     //recMinuit->mnexcm("SIMPlex", arglist, 2, ierflg);
00216 
00217     //Scan is usefull to check FCN
00218     // arglist[0] = 0;
00219     // recMinuit->mnexcm("SCAN", arglist, 1, ierflg);
00220 
00221     //Collect the output
00222     //How to draw the n-sigma contour of a Minuit fit?
00223     //http://root.cern.ch/cgi-bin/print_hit_bold.pl/root/html508/examples/fitcont.C.html?minuit#first_hit
00224 
00225     // record fitter's status
00226     double ln0, edm, errdef;
00227     int nvpar, nparx, icstat;
00228     recMinuit->mnstat(ln0, edm, errdef, nvpar, nparx, icstat);
00229 
00230     // Get Fitted Parameter
00231     double x_fit, x_fit_err;
00232     double y_fit, y_fit_err;
00233     double z_fit, z_fit_err;
00234     double e_fit, e_fit_err;
00235 
00236     recMinuit->GetParameter(0, x_fit, x_fit_err);
00237     recMinuit->GetParameter(1, y_fit, y_fit_err);
00238     recMinuit->GetParameter(2, z_fit, z_fit_err);
00239     recMinuit->GetParameter(3, e_fit, e_fit_err);
00240 
00241     // error matrix; mnemat()
00242     double emat[Npar][Npar];
00243     CLHEP::HepMatrix matrix(Npar, Npar);
00244     recMinuit->mnemat(&emat[0][0], Npar);
00245     for(int row=0; row<Npar; row++) {
00246         for(int col=0; col<Npar; col++) {
00247             matrix[row][col] = emat[row][col];
00248         }
00249     }
00250 
00251     CLHEP::HepLorentzVector vertex(x_fit, y_fit, z_fit, 0.);
00252     recTrigger.setPosition(vertex);
00253     recTrigger.setEnergy(e_fit);
00254     recTrigger.setPositionQuality(ln0);
00255     recTrigger.setEnergyQuality(ln0);
00256     recTrigger.setErrorMatrix(matrix);
00257 
00258     delete recMinuit;
00259     //std::cout << "ncll_fcn called: " << tttttttt << std::endl;
00260 
00261     return StatusCode::SUCCESS;
00262 }
00263 
00264 void QMLFTool::ncll_fcn(int& /*npar*/, double* /*grad*/, double& fval, double* xval, int /*iflag*/)
00265 {
00266     //++tttttttt;
00267     // Always calculate the value of the function, FVAL,
00268     // which is usually a chisquare or log likelihood.
00269 
00270     // Global coefficient (QE, Yield...). Nominal 'absolute' QE: 0.19
00271     //static const double c_n = 200.; // light speed in liquid, mm/s
00272     static const double ECAL = m_opPara.m_photoCathodeArea/4.0/TMath::Pi() * m_opPara.m_yield * 0.19;
00273     static const double fap0 =  0.999946; // 0.9984;
00274     static const double fap1 =  0.101046; // 0.51547;
00275     static const double fap2 = -1.04014;  //-2.2329;
00276     static const double fap3 =  1.01381;  // 2.18043;
00277     static const double fap4 = -0.410953; //-0.78739;
00278     //static const double ftp0 = -0.042;
00279     //static const double ftp1 = 0.11365;
00280     //static const double ftp2 = 1.0;
00281 
00282     double chg_likelihood = 0.;
00283 
00284     // Loop over calib readout channels
00285     for ( int local_id=0; local_id<192; local_id++ ) {
00286         unsigned int pmtid = m_chgLLFUsedPMTs[local_id];
00287         if ( pmtid == 0 ) continue;
00288 
00289         // TODO: optimize, buffer the pos and norm locally
00290         IPmtGeomInfo* geoInfo = m_pmtGeomSvc->get(pmtid);
00291         const CLHEP::Hep3Vector& pos  = geoInfo->localPosition();
00292         const CLHEP::Hep3Vector& norm = geoInfo->localDirection();
00293 
00294         const DayaBay::PmtCalibData* pmtData = m_vPmtData[local_id];
00295 
00296         double dx = xval[0] - pos.x();
00297         double dy = xval[1] - pos.y();
00298         double dz = xval[2] - pos.z();
00299         double dist_norm = dx*norm.x() + dy*norm.y();
00300 
00301         double dxy2  = dx*dx + dy*dy;
00302         double dist2 = dxy2 + dz*dz;
00303         if ( dist2 < 0.0001 ) dist2 = 0.0001;
00304         double dist = sqrt(dist2);
00305 
00306         double cosa_tmp = dist_norm / dist;
00307         if ( cosa_tmp > 1.0 ) cosa_tmp = 1.0;
00308 
00309         double Theta = TMath::ACos(cosa_tmp);
00310         double S_angle = fap0 + Theta*(fap1 + Theta*(fap2 + Theta*(fap3 + fap4*Theta)));
00311         if ( S_angle < 1.0e-5 ) S_angle = 1.0e-5;
00312 
00313         //************* IAV lid permeation *************
00314         //double transTop = 1.0, transBot = 1.0;
00315         //double dz2IAVtop = xval[2] - 1550;
00316         //double rzTop     = dz2IAVtop/dz;
00317         //if ( rzTop > 0 ) {
00318         //    double xIAVtop = xval[0] - rzTop*dx;
00319         //    double yIAVtop = xval[1] - rzTop*dy;
00320         //    double r2IAVtop = xIAVtop*xIAVtop + yIAVtop*yIAVtop;
00321         //    if ( r2IAVtop < 1560*1560 ) {
00322         //      double cosaTop = dz > 0 ? (dz/dist) : (-dz/dist);
00323         //      transTop = ftp2 - ftp2 / (1.0 + exp((cosaTop-ftp0)/ftp1));
00324         //    }
00325         //}
00326         //double dz2IAVbot = xval[2] + 1550;
00327         //double rzBot     = dz2IAVbot/dz;
00328         //if ( rzBot > 0 ) {
00329         //    double xIAVbot = xval[0] - rzBot*dx;
00330         //    double yIAVbot = xval[1] - rzBot*dy;
00331         //    double r2IAVbot = xIAVbot*xIAVbot + yIAVbot*yIAVbot;
00332         //    if ( r2IAVbot < 1560*1560 ) {
00333         //      double cosaBot = dz > 0 ? (dz/dist) : (-dz/dist);
00334         //      transBot = ftp2 - ftp2 / (1.0 + exp((cosaBot-ftp0)/ftp1));
00335         //    }
00336         //}
00337         //std::cout << transTop << " : " << transBot << std::endl;
00338         //S_angle *= transTop*transBot;
00339         //------------- IAV lid permeation -------------
00340 
00341         double expected_q = S_angle*exp(-dist/m_opPara.m_attl)/dist2;
00342 
00343         // Reflection contribution: NREF_ORDER orders are count
00344         double dz_up = -dz, dz_down = dz;
00345         double ref_up = 1.0, ref_down = 1.0;
00346         double pos2refX2[2] = { 2*(m_geomPara.m_topRefZ - pos.z()),
00347                                 2*(pos.z() - m_geomPara.m_botRefZ) };
00348         for ( int rOrder = 0; rOrder < NREF_ORDER; ++rOrder ) {
00349             int up_down_flag = (rOrder & 1);  // rOrder % 2
00350 
00351             // The init direction points up
00352             dz_up += pos2refX2[up_down_flag];
00353             dist2 = dxy2 + dz_up*dz_up;
00354             if ( dist2 < 0.0001 ) dist2 = 0.0001;
00355             dist = sqrt(dist2);
00356 
00357             cosa_tmp = dist_norm / dist;
00358             if ( cosa_tmp > 1.0 ) cosa_tmp = 1.0;
00359             Theta = TMath::ACos(cosa_tmp);
00360             S_angle = fap0 + Theta*(fap1 + Theta*(fap2 + Theta*(fap3 + fap4*Theta)));
00361             if ( S_angle < 1.0e-5 ) S_angle = 1.0e-5;
00362 
00363             ref_up = up_down_flag ? ref_up*m_opPara.m_botRef : ref_up*m_opPara.m_topRef;
00364             expected_q += ref_up*S_angle*exp(-dist/m_opPara.m_attl)/dist2;
00365 
00366             // The init direction points down
00367             dz_down += pos2refX2[1-up_down_flag];
00368             dist2 = dxy2 + dz_down*dz_down;
00369             if ( dist2 < 0.0001 ) dist2 = 0.0001;
00370             dist = sqrt(dist2);
00371 
00372             cosa_tmp = dist_norm / dist;
00373             if ( cosa_tmp > 1.0 ) cosa_tmp = 1.0;
00374             Theta = TMath::ACos(cosa_tmp);
00375             S_angle = fap0 + Theta*(fap1 + Theta*(fap2 + Theta*(fap3 + fap4*Theta)));
00376             if ( S_angle < 1.0e-5) S_angle = 1.0e-5;
00377 
00378             ref_down = up_down_flag ? ref_down*m_opPara.m_topRef : ref_down*m_opPara.m_botRef;
00379             expected_q += ref_down*S_angle*exp(-dist/m_opPara.m_attl)/dist2;
00380         }
00381 
00382         expected_q *= xval[3] * ECAL * m_pmtEff[m_ADNo][local_id];  //pmtData->m_efficiency;
00383 
00384         if ( m_LLFMode == 1 ) {
00386             double iCharge = m_obsq[local_id];
00387             if ( iCharge != 0 ) {
00388                 chg_likelihood += 2*((expected_q-iCharge) + log(iCharge/expected_q)*iCharge);
00389             }
00390             else {
00391                 chg_likelihood += 2*expected_q;
00392             }
00393         }
00394         else {
00397         }
00398     } // loop for local_id
00399 
00400     fval = chg_likelihood;
00401 }
00402 
00403 void QMLFTool::getOpPara()
00404 {
00405     if ( m_opLocation == "DetDesc" ) {
00406         // Get Optical Parameters from Det/DetDesc
00407         // A tabular was defined in DetDesc, but only one value was get here
00408         // FIXME
00409         // Get the Detector Data Service to access TDS:
00410         // This should be done in your tools initialize()
00411         // Later, look up some material.  If this material is always used you
00412         // should look it up once in initialize() and save it for later use
00413 
00414         if ( m_opPara.m_attl < 0 ) { // absorption length
00415             m_opPara.m_attl = meanOpticalPara<Material>(
00416                     "/dd/Materials/GdDopedLS",
00417                     "ABSLENGTH" );
00418         }
00419         if ( m_opPara.m_topRef < 0 ) { // top ESR reflectivity
00420             m_opPara.m_topRef = meanOpticalPara<Surface>(
00421                     "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceTop",
00422                     "REFLECTIVITY" );
00423         }
00424         if ( m_opPara.m_botRef < 0 ) { // bottom ESR reflectivity
00425             m_opPara.m_botRef = meanOpticalPara<Surface>(
00426                     "/dd/Geometry/AdDetails/AdSurfacesAll/ESRAirSurfaceBot",
00427                     "REFLECTIVITY" );
00428         }
00429     }
00430     else if ( m_opLocation == "OpCalibSvc" ) {
00431     }
00432     else {
00433         warning() <<"No external optical data source, use default value."<< endreq;
00434     }
00435 }
00436 
00437 void QMLFTool::getGeomPara()
00438 {
00439     if ( m_geomPara.m_topRefZ != 0. && m_geomPara.m_botRefZ != 0. ) return;
00440 
00441     if ( m_geomLocation == "DetDesc" ) {
00442         IDataProviderSvc* dds = 0;
00443         StatusCode sc = service( "DetectorDataSvc", dds, true );
00444         if (sc.isFailure()) {
00445             error() << "No DetectorDataSvc available!" << endreq;
00446         }
00447 
00448         // define a pointer to ICoordSysSvc to
00449         // use to cache the service during the job.
00450         ICoordSysSvc* css;
00451         sc = service( "CoordSysSvc", css, true );
00452 
00453         Gaudi::XYZPoint zero(0,0,0);
00454 
00455         if ( m_geomPara.m_topRefZ == 0. ) {
00456             std::string topEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-topesr";
00457             SmartDataPtr<IDetectorElement> topEsr(dds,topEsrPath);
00458             // Set top ESR position
00459             Gaudi::XYZPoint globalPoint = topEsr->geometry()->toGlobal(zero);
00460             IDetectorElement* theAdDe = css->coordSysDE(globalPoint);
00461             Gaudi::XYZPoint adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
00462             m_geomPara.m_topRefZ = adLocalPoint.z();
00463         }
00464 
00465         if ( m_geomPara.m_botRefZ == 0. ) {
00466             std::string botEsrPath = "/dd/Structure/AdReflectorStructure/db-ad/db-ad1-botesr";
00467             SmartDataPtr<IDetectorElement> botEsr(dds,botEsrPath);
00468             // Set bot ESR position
00469             Gaudi::XYZPoint globalPoint = botEsr->geometry()->toGlobal(zero);
00470             IDetectorElement* theAdDe = css->coordSysDE(globalPoint);
00471             Gaudi::XYZPoint adLocalPoint = theAdDe->geometry()->toLocal(globalPoint);
00472             m_geomPara.m_botRefZ = adLocalPoint.z();
00473         }
00474     }
00475     else if ( m_geomLocation == "OnsiteSuveyData" ) {
00476         // Get geometry info from onsite suvey data
00477     }
00478     else {
00479         warning() <<"Unknown geometry data source, use default value."<< endreq;
00480     }
00481 }
00482 
00483 void QMLFTool::getPmtEff()
00484 {
00485     std::string fname = getenv("ADRECROOT");
00486     fname += "/share/AD-PMT-RQE.txt";
00487     std::ifstream fs(fname.c_str());
00488 
00489     int tmp;
00490     for ( int j = 0; j < 192; ++j ) {
00491         fs >> tmp >> tmp;
00492         for ( int i = 0; i < 8; ++i ) {
00493             fs >> m_pmtEff[i][j];
00494         }
00495     }
00496 }
00497 
00498 void QMLFTool::printPara()
00499 {
00500     info() << "Op and Geom parameters used in reconstruction." << endreq;
00501     info() << "Attenuation length : " << m_opPara.m_attl <<"mm"<< endreq;
00502     info() << "Top reflectivity : " << m_opPara.m_topRef << endreq;
00503     info() << "Bottom reflectivity : " << m_opPara.m_botRef << endreq;
00504     info() << "Top reflector Z pos :" << m_geomPara.m_topRefZ <<"mm"<< endreq;
00505     info() << "Bot reflector Z pos :" << m_geomPara.m_botRefZ <<"mm"<< endreq;
00506 }
00507 
00508 template<typename MType>
00509 double QMLFTool::meanOpticalPara(const std::string& location, const std::string& type)
00510 {
00511     IDataProviderSvc* dds = 0;
00512     StatusCode sc = service( "DetectorDataSvc", dds, true );
00513     if (sc.isFailure()) {
00514         error() << "No DetectorDataSvc available!" << endreq;
00515         return -1.;
00516     }
00517 
00518     MType* element = GaudiCommon<AlgTool>::get<MType>(dds, location);
00519     typename MType::Tables& element_tab = element->tabulatedProperties();
00520 
00521     typename MType::Tables::const_iterator tbIter = element_tab.begin(); 
00522     while ( tbIter != element_tab.end() ) {
00523         if ( (*tbIter)->type() == type ) break;
00524         ++tbIter;
00525     }
00526     const TabulatedProperty::Table& par_tab = (*tbIter)->table(); 
00527 
00528     // look up Tables, get an average value for the parameter
00529     int count = 0;
00530     double parSum = 0.0;
00531     for ( TabulatedProperty::Table::const_iterator parIter = par_tab.begin();
00532             parIter != par_tab.end();
00533             parIter++ )
00534     {
00535         if(parIter->first<3.18e-06 && parIter->first>2.07e-06 ) { // 390nm ~ 600nm
00536             count++;
00537             parSum += parIter->second;
00538         }
00539     }
00540 
00541     return parSum/count; // average value between 390nm~600nm
00542 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:13:25 for AdRec by doxygen 1.7.4