Inheritance diagram for CalValsTool:
Definition at line 61 of file CalValsTool.cxx.
|
||||||||||||||||
|
Definition at line 270 of file CalValsTool.cxx.
00273 : ValBase( type, name, parent ) 00274 { 00275 // Declare additional interface 00276 declareInterface<IValsTool>(this); 00277 } |
|
|
Definition at line 69 of file CalValsTool.cxx.
00069 { }
|
|
||||||||||||
|
Definition at line 1345 of file CalValsTool.cxx. References ValBase::globalToLocal(), m_towerPitch, m_xNum, and m_yNum. Referenced by calculate().
01346 {
01347 double edge = 0.;
01348 double x = pos.x();
01349 double y = pos.y();
01350 double x_twr = globalToLocal(x, m_towerPitch, m_xNum);
01351 double y_twr = globalToLocal(y, m_towerPitch, m_yNum);
01352
01353 if( fabs(x_twr) > fabs(y_twr) ) {
01354 edge = m_towerPitch/2. - fabs(x_twr);
01355 view = 0;
01356 }
01357 else {
01358 edge = m_towerPitch/2. - fabs(y_twr);
01359 view = 1;
01360 }
01361 return edge;
01362 }
|
|
|
calculate all values; implemented by each XxxValsTool
Reimplemented from ValBase. Definition at line 667 of file CalValsTool.cxx. References activeDist(), CAL_BkHalf_Ratio, CAL_cfp_alpha, CAL_cfp_calEffRLn, CAL_cfp_calfit_alpha, CAL_cfp_calfit_calEffRLn, CAL_cfp_calfit_energy, CAL_cfp_calfit_fiterrflg, CAL_cfp_calfit_tmax, CAL_cfp_calfit_totChiSq, CAL_cfp_energy, CAL_cfp_fiterrflg, CAL_cfp_tkrRLn, CAL_cfp_tmax, CAL_cfp_totChiSq, CAL_Cntr_RLn, CAL_CsI_RLn, CAL_DeadCnt_Rat, CAL_DeadTot_Rat, CAL_deltaT, CAL_eAveBack, CAL_Edge_Corr, CAL_EdgeEnergy, CAL_eLayer, CAL_EnergyCorr, CAL_EnergyRaw, CAL_Gap_Fraction, CAL_LAT_RLn, CAL_LATEdge, CAL_layer0Ratio, CAL_Leak_Corr, CAL_LkHd_energy, CAL_LkHd_energyErr, CAL_lll_energy, CAL_lll_energyErr, CAL_Long_Rms, CAL_LRms_Asym, CAL_Lyr0_Ratio, CAL_Lyr7_Ratio, CAL_Max_Num_Xtals_In_Layer, CAL_MIP_Diff, CAL_MIP_Ratio, CAL_nLayersRmsBack, CAL_nsaturated, CAL_Num_Xtals, CAL_Num_Xtals_Trunc, CAL_posdir_chisq, CAL_posdir_nlayers, CAL_rmsLayerE, CAL_rmsLayerEBack, CAL_t_Pred, CAL_tkl_energy, CAL_tkl_energyErr, CAL_Top_Gap_Dist, CAL_Tot_RLn, CAL_Total_Corr, CAL_Track_Angle, CAL_Track_DOCA, CAL_track_E_rms, CAL_track_E_rms_trunc, CAL_track_rms, CAL_track_rms_trunc, CAL_Track_Sep, CAL_Trans_Rms, CAL_TS_CAL2_T_100, CAL_TS_CAL2_T_68, CAL_TS_CAL2_T_90, CAL_TS_CAL2_T_95, CAL_TS_CAL2_T_99, CAL_TS_CAL2_TL_100, CAL_TS_CAL2_TL_68, CAL_TS_CAL2_TL_90, CAL_TS_CAL2_TL_95, CAL_TS_CAL2_TL_99, CAL_TS_CAL_T_100, CAL_TS_CAL_T_68, CAL_TS_CAL_T_90, CAL_TS_CAL_T_95, CAL_TS_CAL_T_99, CAL_TS_CAL_TL_100, CAL_TS_CAL_TL_68, CAL_TS_CAL_TL_90, CAL_TS_CAL_TL_95, CAL_TS_CAL_TL_99, CAL_TS_TKR_T_100, CAL_TS_TKR_T_68, CAL_TS_TKR_T_90, CAL_TS_TKR_T_95, CAL_TS_TKR_T_99, CAL_TS_TKR_TL_100, CAL_TS_TKR_TL_68, CAL_TS_TKR_TL_90, CAL_TS_TKR_TL_95, CAL_TS_TKR_TL_99, CAL_TwrEdgeCntr, CAL_TwrEdgeTop, CAL_x0, CAL_xdir, CAL_xdir2, CAL_xEcntr, CAL_xEcntr2, CAL_xPosRmsLastLayer, CAL_Xtal_maxEne, CAL_Xtal_Ratio, CAL_y0, CAL_ydir, CAL_ydir2, CAL_yEcntr, CAL_yEcntr2, CAL_yPosRmsLastLayer, CAL_zdir, CAL_zdir2, CAL_zEcntr, CAL_zEcntr2, Doca::docaOfPoint(), m_calXWidth, m_calYWidth, m_calZBot, m_calZTop, m_detSvc, m_G4PropTool, m_nCsI, m_nLayers, ValBase::m_pEventSvc, m_tkrGeom, m_towerPitch, m_xNum, m_yNum, ValBase::mapIter, ValBase::printHeader(), ValBase::setAnaTupBit(), TSaxisP, TSaxisV, TSfillTS(), TSfillTSdist(), TSgetinterpolationTS(), TSnlog, and TSTS.
00668 {
00669 StatusCode sc = StatusCode::SUCCESS;
00670
00671 // Recover Track associated info.
00672 SmartDataPtr<Event::TkrTrackCol>
00673 pTracks(m_pEventSvc,EventModel::TkrRecon::TkrTrackCol);
00674 SmartDataPtr<Event::TkrVertexCol>
00675 pVerts(m_pEventSvc,EventModel::TkrRecon::TkrVertexCol);
00676
00677 // Recover pointers to CalClusters and Xtals
00678 SmartDataPtr<Event::CalClusterCol>
00679 pCals(m_pEventSvc,EventModel::CalRecon::CalClusterCol);
00680 SmartDataPtr<Event::CalXtalRecCol>
00681 pxtalrecs(m_pEventSvc,EventModel::CalRecon::CalXtalRecCol);
00682
00683 // Recover pointer to CalEventEnergy info
00684 #ifdef PRE_CALMOD
00685 Event::CalEventEnergy* calEventEnergy =
00686 SmartDataPtr<Event::CalEventEnergy>(m_pEventSvc, EventModel::CalRecon::CalEventEnergy);
00687 #else
00688 Event::CalEventEnergyCol * calEventEnergyCol =
00689 SmartDataPtr<Event::CalEventEnergyCol>(m_pEventSvc,EventModel::CalRecon::CalEventEnergyCol) ;
00690 Event::CalEventEnergy * calEventEnergy = 0 ;
00691 if ((calEventEnergyCol!=0)&&(!calEventEnergyCol->empty()))
00692 calEventEnergy = calEventEnergyCol->front() ;
00693 #endif
00694
00695
00696 // If calEventEnergy then fill TkrEventParams
00697 // Note: TkrEventParams initializes to zero in the event of no CalEventEnergy
00698 double m_radLen_Stuff = 0.;
00699 double m_radLen_CntrStuff = 0.;
00700 if (calEventEnergy != 0)
00701 {
00702 // Extraction of results from CalValCorrTool in CalRecon...
00703 for(Event::CalCorToolResultCol::iterator corIter = calEventEnergy->begin(); corIter != calEventEnergy->end(); corIter++)
00704 {
00705 Event::CalCorToolResult corResult = **corIter;
00706
00707 if (corResult.getCorrectionName() == "CalValsCorrTool")
00708 {
00709 CAL_EnergyCorr = corResult["CorrectedEnergy"];
00710 CAL_Leak_Corr = corResult["LeakCorrection"];
00711 CAL_Edge_Corr = corResult["EdgeCorrection"];
00712 CAL_Gap_Fraction = corResult["GapFraction"];
00713 CAL_Total_Corr = corResult["TotalCorrection"];
00714
00715 CAL_CsI_RLn = corResult["CsIRLn"];
00716 CAL_LAT_RLn = corResult["LATRLn"];
00717 m_radLen_Stuff = corResult["StuffRLn"];
00718 m_radLen_CntrStuff = corResult["CntrRLnStuff"];
00719 CAL_Cntr_RLn = corResult["CntrRLn"];
00720 CAL_t_Pred = corResult["PredCntr"];
00721 CAL_deltaT = corResult["DeltaCntr"];
00722 CAL_Tot_RLn = corResult["CALRLn"];
00723
00724 // CAL_x0 = corResult["CalTopX0"]; See below under Imaging Calorimeter Top
00725 // CAL_y0 = corResult["CalTopY0"];
00726
00727 //CAL_RmsE = corResult["RmsE"];
00728 //CAL_numLayersRms = corResult["NumLayersRms"];
00729 }
00730 else if (corResult.getCorrectionName() == "CalFullProfileTool")
00731 {
00732 CAL_cfp_energy = corResult.getParams().getEnergy();
00733 CAL_cfp_totChiSq = corResult["totchisq"];
00734 CAL_cfp_calEffRLn = corResult["cal_eff_RLn"];
00735 CAL_cfp_tkrRLn = corResult["tkr_RLn"];
00736 CAL_cfp_alpha = corResult["alpha"];
00737 CAL_cfp_tmax = corResult["tmax"];
00738 CAL_cfp_fiterrflg = corResult["fitflag"];
00739 CAL_cfp_calfit_energy = corResult["calfit_fit_energy"];
00740 CAL_cfp_calfit_totChiSq = corResult["calfit_totchisq"];
00741 CAL_cfp_calfit_calEffRLn = corResult["calfit_cal_eff_RLn"];
00742 CAL_cfp_calfit_alpha = corResult["calfit_alpha"];
00743 CAL_cfp_calfit_tmax = corResult["calfit_tmax"];
00744 CAL_cfp_calfit_fiterrflg = corResult["calfit_fitflag"];
00745 }
00746 else if (corResult.getCorrectionName() == "CalLastLayerLikelihoodTool")
00747 {
00748 CAL_lll_energy = corResult.getParams().getEnergy();
00749 CAL_lll_energyErr = corResult.getParams().getEnergyErr();
00750 }
00751 else if (corResult.getCorrectionName() == "CalTkrLikelihoodTool")
00752 {
00753 CAL_tkl_energy = corResult.getParams().getEnergy();
00754 CAL_tkl_energyErr = corResult.getParams().getEnergyErr();
00755 }
00756 else if (corResult.getCorrectionName() == "CalLikelihoodManagerTool")
00757 {
00758 CAL_LkHd_energy = corResult.getParams().getEnergy();
00759 CAL_LkHd_energyErr = corResult.getParams().getEnergyErr();
00760 }
00761 }
00762 }
00763
00764 //Make sure we have valid cluster data and some energy
00765 if (!pCals) return sc;
00766 if (pCals->empty()) return sc;
00767 Event::CalCluster* calCluster = pCals->front();
00768 CAL_EnergyRaw = calCluster->getCalParams().getEnergy();
00769 if(CAL_EnergyRaw<1.0) return sc;
00770
00771 for(int i = 0; i<m_nLayers; i++) CAL_eLayer[i] = (*calCluster)[i].getEnergy();
00772
00773 CAL_Trans_Rms = sqrt(calCluster->getRmsTrans()/CAL_EnergyRaw);
00774
00775 float logRLn = 0;
00776 if ((CAL_LAT_RLn - CAL_Cntr_RLn) > 0.0) {
00777 logRLn = log(CAL_LAT_RLn - CAL_Cntr_RLn);
00778 }
00779 if (logRLn > 0.0) {
00780 CAL_Long_Rms = sqrt(calCluster->getRmsLong()/CAL_EnergyRaw) / logRLn;
00781 }
00782 CAL_LRms_Asym = calCluster->getRmsLongAsym();
00783
00784 if(CAL_EnergyRaw>0.0) {
00785 CAL_Lyr0_Ratio = CAL_eLayer[0]/CAL_EnergyRaw;
00786 if(m_nLayers==8) {
00787 CAL_Lyr7_Ratio = CAL_eLayer[7]/CAL_EnergyRaw;
00788 CAL_BkHalf_Ratio = (CAL_eLayer[4]+CAL_eLayer[5]+
00789 CAL_eLayer[6]+CAL_eLayer[7])/CAL_EnergyRaw;
00790 }
00791 }
00792
00793 int no_xtals=0;
00794 CAL_Xtal_maxEne = 0.;
00795
00796 // Local array for #xtals in layer with the most xtals (cut on e>=5MeV);
00797 std::vector<int> xtalCount(m_nLayers,0);
00798
00799 Event::CalXtalRecCol::const_iterator jlog;
00800 if (pxtalrecs) {
00801 // Find Xtal with max. energy
00802 CAL_Num_Xtals = (float)pxtalrecs->size();
00803 for( jlog=pxtalrecs->begin(); jlog != pxtalrecs->end(); ++jlog) {
00804 const Event::CalXtalRecData& recLog = **jlog;
00805 double eneLog = recLog.getEnergy();
00806 if(eneLog > CAL_Xtal_maxEne) CAL_Xtal_maxEne = eneLog;
00807 idents::CalXtalId xtalId = recLog.getPackedId();
00808 int layer = xtalId.getLayer();
00809 if(eneLog>5.0) xtalCount[layer]++;
00810 }
00811
00812 std::vector<int>::const_iterator itC =
00813 std::max_element(xtalCount.begin(),xtalCount.end());
00814 CAL_Max_Num_Xtals_In_Layer = *itC;
00815
00816 // Number of Xtals
00817 no_xtals=pxtalrecs->size();
00818 }
00819 int no_xtals_trunc=calCluster->getNumTruncXtals();
00820 CAL_Xtal_Ratio= (no_xtals>0) ? float(no_xtals_trunc)/no_xtals : 0;
00821 CAL_Num_Xtals_Trunc = float(no_xtals_trunc);
00822
00823 // No use in continuing if too little energy in CAL
00824 if(CAL_EnergyRaw < 5.) return sc;
00825
00826 Point cal_pos = calCluster->getPosition();
00827 Vector cal_dir = calCluster->getDirection();
00828
00829 CAL_xEcntr = cal_pos.x();
00830 CAL_yEcntr = cal_pos.y();
00831 CAL_zEcntr = cal_pos.z();
00832 CAL_xdir = cal_dir.x();
00833 CAL_ydir = cal_dir.y();
00834 CAL_zdir = cal_dir.z();
00835
00836 // Get pos and dir determined using only the transverse position information
00837 CAL_xEcntr2 = calCluster->getCalParams().getxPosxPos();
00838 CAL_yEcntr2 = calCluster->getCalParams().getyPosyPos();
00839 CAL_zEcntr2 = calCluster->getCalParams().getzPoszPos();
00840 CAL_posdir_chisq = calCluster->getCalParams().getxPosyPos();
00841 CAL_posdir_nlayers = calCluster->getCalParams().getxPoszPos();
00842 CAL_nsaturated = calCluster->getCalParams().getyPoszPos();
00843 CAL_xdir2 = calCluster->getCalParams().getxDirxDir();
00844 CAL_ydir2 = calCluster->getCalParams().getyDiryDir();
00845 CAL_zdir2 = calCluster->getCalParams().getzDirzDir();
00846
00847 // Get the lower and upper limits for the CAL in the installed towers
00848 double deltaX = 0.5*(m_xNum*m_towerPitch - m_calXWidth);
00849 double deltaY = 0.5*(m_yNum*m_towerPitch - m_calYWidth);
00850 double calXLo = m_tkrGeom->getLATLimit(0,LOW) + deltaX;
00851 double calXHi = m_tkrGeom->getLATLimit(0,HIGH) - deltaX;
00852 double calYLo = m_tkrGeom->getLATLimit(1,LOW) + deltaY;
00853 double calYHi = m_tkrGeom->getLATLimit(1,HIGH) - deltaY;
00854
00855 // Event axis locations relative to towers and LAT
00856 Point pos0(CAL_x0, CAL_y0, m_calZTop);
00857 int view;
00858 CAL_TwrEdgeTop = activeDist(pos0, view);
00859 CAL_TwrEdgeCntr = activeDist(cal_pos, view);
00860 // Find the distance closest to an edge
00861 double dX = std::max(calXLo-CAL_x0, CAL_x0-calXHi);
00862 double dY = std::max(calYLo-CAL_y0, CAL_y0-calYHi);
00863 CAL_LATEdge = -std::max(dX, dY);
00864
00865 // collect the CAL edge energy
00866 // the edge is larger of the width and length of the layer
00867 // we can do better if this is at all useful.
00868 if (pxtalrecs) {
00869
00870 // do the Cal[X/Y]PosRmsLL here
00871 double eRmsLL = 0;
00872 double xLL = 0;
00873 double xSqLL = 0;
00874 double yLL = 0;
00875 double ySqLL = 0;
00876 int nRmsLL = 0;
00877 bool doLL = (m_nLayers>1 && m_nCsI>1); // true for Glast, false for EGRET
00878
00879 for( jlog=pxtalrecs->begin(); jlog != pxtalrecs->end(); ++jlog) {
00880 const Event::CalXtalRecData& recLog = **jlog;
00881 Point pos = recLog.getPosition();
00882 double eneLog = recLog.getEnergy();
00883 double xPos = pos.x(); double yPos = pos.y();
00884 double minX, minY;
00885 minX = std::min(fabs(xPos-calXLo),fabs(xPos-calXHi));
00886 minY = std::min(fabs(yPos-calYLo),fabs(yPos-calYHi));
00887
00888 // for last-layer rms
00889 idents::CalXtalId id = recLog.getPackedId();
00890 int layer = id.getLayer();
00891 if (eneLog>0 && layer==m_nLayers-1 && doLL) {
00892 nRmsLL++;
00893 eRmsLL += eneLog;
00894 xLL += eneLog*xPos;
00895 yLL += eneLog*yPos;
00896 xSqLL += eneLog*xPos*xPos;
00897 ySqLL += eneLog*yPos*yPos;
00898 }
00899
00900 // for later... no time to check this out now!
00901 /*
00902 idents::CalXtalId id = recLog.getPackedId();
00903 if (id.isX()) {
00904 double minX =
00905 std::min(fabs(xPos-(calXLo+m_deltaSide)),fabs(xPos-(calXHi-m_deltaSide)));
00906 double minY = std::min(fabs(yPos-calYLo),fabs(yPos-calYHi));
00907 } else {
00908 double minX = std::min(fabs(xPos-calXLo),fabs(xPos-calXHi));
00909 double minY =
00910 std::min(fabs(yPos-(calYLo+m_deltaSide)),fabs(yPos-(calYHi-m_deltaSide)));
00911 }
00912 */
00913 if ((minX<_deltaEdge) || (minY<_deltaEdge)) {
00914 double eneLog = recLog.getEnergy();
00915 CAL_EdgeEnergy += eneLog;
00916 }
00917 }
00918
00919 // finish last-layer rms
00920 if(nRmsLL>1) {
00921 double xAveLL = xLL/eRmsLL;
00922 double xVarLL = std::max(0.0, xSqLL/eRmsLL - xAveLL*xAveLL);
00923 double yAveLL = yLL/eRmsLL;
00924 double yVarLL = std::max(0.0, ySqLL/eRmsLL - yAveLL*yAveLL);
00925 CAL_xPosRmsLastLayer = (float) sqrt( xVarLL*nRmsLL/(nRmsLL-1));
00926 CAL_yPosRmsLastLayer = (float) sqrt( yVarLL*nRmsLL/(nRmsLL-1));
00927 }
00928 }
00929
00930 // Compare Track (or vertex) to Cal soltion
00931 Point x0 = cal_pos;
00932 Vector t0 = -cal_dir;
00933 if(calCluster->getRmsLong() < .1 ) { // Trap no-calc. condition
00934 x0 = Point(0., 0., 0.);
00935 t0 = Vector(0., 0., -1.);
00936 }
00937
00938 // If there are tracks, use the best one
00939 int num_tracks = 0;
00940 Point x1, x2;
00941 Vector t1, t2;
00942
00943 if(pTracks) num_tracks = pTracks->size();
00944
00945 Event::TkrTrackColConPtr pTrack1;
00946 Event::TkrTrack* track_1;
00947
00948 if(num_tracks > 0) {
00949 // Get the first track
00950 pTrack1 = pTracks->begin();
00951 track_1 = *pTrack1;
00952
00953 // Get the start and direction
00954 x1 = track_1->getInitialPosition();
00955 t1 = track_1->getInitialDirection();
00956 x0 = x1;
00957 t0 = t1;
00958
00959 // Find the distance for energy centroid to track axis
00960 /* old calculation
00961 Vector x_diff = x0 - cal_pos;
00962 double x_diff_sq = x_diff*x_diff;
00963 double x_diff_t0 = x_diff*t0;
00964 CAL_Track_DOCA = sqrt(std::max(0.0, x_diff_sq - x_diff_t0*x_diff_t0));
00965 */
00966 Doca track1(x1, t1);
00967 CAL_Track_DOCA = (float)track1.docaOfPoint(cal_pos);
00968
00969 // Image the top of the Calorimeter - use last hit FILTERED
00970 const Event::TkrTrackParams& tkr1_params = track_1->back()->getTrackParams(Event::TkrTrackHit::FILTERED);
00971 double xSlp = tkr1_params.getxSlope();
00972 double ySlp = tkr1_params.getySlope();
00973 Vector LastDir = Vector(-xSlp, -ySlp, -1).unit();
00974 Point LastHit = track_1->back()->getPoint(Event::TkrTrackHit::FILTERED);
00975
00976 //double calZTop = -46.;
00977 double deltaZ = m_calZTop - LastHit.z();
00978 double topArcLength = deltaZ/LastDir.z();
00979 Point calTopLoc = LastHit + topArcLength*LastDir;
00980 CAL_x0 = calTopLoc.x();
00981 CAL_y0 = calTopLoc.y();
00982
00983 // Now create an active distance type variable using the tower pitch
00984 double integer_part;
00985 double deltaX_crack = modf(fabs(CAL_x0)/m_towerPitch, &integer_part);
00986 if(deltaX_crack > .5) deltaX_crack = 1. - deltaX_crack;
00987 double deltaY_crack = modf(fabs(CAL_y0)/m_towerPitch, &integer_part);
00988 if(deltaY_crack > .5) deltaY_crack = 1. - deltaY_crack;
00989 CAL_Top_Gap_Dist = m_towerPitch*std::min(deltaX_crack, deltaY_crack);
00990
00991 if(num_tracks > 1) {
00992 // Get the second track
00993 pTrack1++;
00994 Event::TkrTrack* track_1 = *pTrack1;
00995
00996 // Get the start and direction
00997 x2 = track_1->getInitialPosition();
00998 t2 = track_1->getInitialDirection();
00999
01000 Point x2Top = x2 + (pos0.z()-x2.z())/t2.z() * t2;
01001 CAL_Track_Sep = (x2Top - pos0).mag();
01002 }
01003
01004
01005 // If vertexed - use first vertex
01006 if(pVerts && pVerts->size()>0) {
01007 Event::TkrVertex* vertex = *(pVerts->begin());
01008 x0 = vertex->getPosition();
01009 t0 = vertex->getDirection();
01010 }
01011
01012 // try Bill's new vars...
01013 if (pxtalrecs) {
01014 //make a map of xtal energy by doca
01015 // we need this so that we can drop the largest entries for the truncated calc.
01016 std::multimap<double, double> docaMap;
01017 std::multimap<double, double>::iterator mapIter;
01018
01019 for( jlog=pxtalrecs->begin(); jlog != pxtalrecs->end(); ++jlog) {
01020 const Event::CalXtalRecData& recLog = **jlog;
01021 Point pos = recLog.getPosition();
01022 double doca = track1.docaOfPoint(pos);
01023 double energy = recLog.getEnergy();
01024
01025 std::pair<double, double> docaPair(doca, energy);
01026 docaMap.insert(docaPair);
01027 }
01028
01029 unsigned int mapSize = docaMap.size();
01030 int truncSize = std::min((int)((mapSize*truncFraction)+0.5),(int) mapSize);
01031
01032 // make up the vars by iterating over the map
01033 // for the truncated var, stop after truncFraction of the xtals
01034 int mapCount = 0;
01035 // 4 measures of the Track-Cal RMS
01036 float totalE1 = 0.0;
01037 float totalE2 = 0.0;
01038 for(mapIter=docaMap.begin(); mapIter!=docaMap.end();++mapIter) {
01039 ++mapCount;
01040 double doca = mapIter->first;
01041 double ene = mapIter->second;
01042 double docaSq = doca*doca;
01043 CAL_track_rms += (float) docaSq;
01044 CAL_track_E_rms += (float) docaSq*ene;
01045 totalE1 += (float) ene;
01046 if (mapCount<=truncSize) {
01047 CAL_track_rms_trunc += (float) docaSq;
01048 CAL_track_E_rms_trunc += (float) docaSq*ene;
01049 totalE2 += (float) ene;
01050 }
01051 }
01052 CAL_track_rms = sqrt(CAL_track_rms/mapSize);
01053 if (totalE1>0.0) {
01054 CAL_track_E_rms =
01055 sqrt(std::max(0.0f, CAL_track_E_rms)/totalE1);
01056 }
01057 CAL_track_rms_trunc = sqrt(CAL_track_rms_trunc/truncSize);
01058 if (totalE2>0.0) {
01059 CAL_track_E_rms_trunc =
01060 sqrt(std::max(0.0f, CAL_track_E_rms_trunc)/totalE2);
01061 }
01062 }
01063 }
01064
01065 // Find angle between Track and Cal. Moment axis
01066 // Note: the direction in Cal is opposite to tracking!
01067 if(num_tracks>0 && fabs(cal_dir.x()) < 1.) {
01068 double cosCalt0 = std::min(1., -t0*cal_dir);
01069 cosCalt0 = std::max(cosCalt0, -1.); // just in case...
01070 CAL_Track_Angle = acos(cosCalt0);
01071 } else {
01072 CAL_Track_Angle = -.1f;
01073 }
01074
01075 // Energy compared to muon
01076 CAL_MIP_Diff = CAL_EnergyRaw- 12.07*CAL_CsI_RLn;
01077 const double minRadLen = 0.1;
01078 CAL_MIP_Ratio = CAL_EnergyRaw/(12.07*std::max(CAL_CsI_RLn*1., minRadLen));
01079
01080 // Ratios of rad. len. in material other then CsI
01081 CAL_DeadTot_Rat = m_radLen_Stuff/std::max(minRadLen, CAL_Tot_RLn*1.);
01082 CAL_DeadCnt_Rat = m_radLen_CntrStuff/std::max(minRadLen, CAL_Cntr_RLn*1.);
01083
01084 // Here we do the CAL_RmsE calculation
01085
01086 Point xEnd;
01087 Vector tEnd;
01088 double arclen;
01089
01090 try {
01091 // get the last point on the best track
01092 if(num_tracks>0) {
01093 Event::TkrTrackHitVecConItr hitIter = (*track_1).end();
01094 hitIter--;
01095 const Event::TkrTrackHit* hit = *hitIter;
01096 xEnd = hit->getPoint(Event::TkrTrackHit::SMOOTHED);
01097 tEnd = hit->getDirection(Event::TkrTrackHit::SMOOTHED);
01098 double eCosTheta = fabs(tEnd.z());
01099
01100 // find the parameters to start the swim
01101 double deltaZ = xEnd.z() - m_calZBot;
01102 arclen = fabs(deltaZ/eCosTheta);
01103
01104 // do the swim
01105 m_G4PropTool->setStepStart(xEnd, tEnd);
01106 m_G4PropTool->step(arclen);
01107
01108 // collect the radlens by layer
01109 int numSteps = m_G4PropTool->getNumberSteps();
01110 std::vector<double> rlCsI(m_nLayers, 0.0);
01111 std::vector<bool> useLayer(m_nLayers, true);
01112 idents::VolumeIdentifier volId;
01113 idents::VolumeIdentifier prefix = m_detSvc->getIDPrefix();
01114 int istep = 0;
01115 for(; istep < numSteps; ++istep) {
01116 volId = m_G4PropTool->getStepVolumeId(istep);
01117 volId.prepend(prefix);
01118 bool inXtal = ( volId.size()>7 && volId[0]==0
01119 && volId[3]==0 && volId[7]==0 ? true : false );
01120 if(inXtal) {
01121 int layer = volId[4];
01122 double radLen_step = m_G4PropTool->getStepRadLength(istep);
01123 rlCsI[layer] += radLen_step;
01124 }
01125 }
01126
01127 double eTot = 0;
01128 double eNorm0 = 0;
01129 double e2 = 0;
01130 int CAL_nLayersRms = 0;
01131 int layer;
01132
01133 for (layer=0; layer<m_nLayers; ++layer) {
01134 if (rlCsI[layer]<0.5) useLayer[layer] = false;
01135 if (CAL_eLayer[layer]<0.05*CAL_EnergyRaw || CAL_EnergyRaw<=0) {
01136 useLayer[layer] = false;
01137 }
01138 }
01139
01140 for (layer=0; layer<m_nLayers; ++layer) {
01141 if(!useLayer[layer]) continue;
01142 CAL_nLayersRms++;
01143 double eNorm = CAL_eLayer[layer]/rlCsI[layer];
01144 if(layer==0) {
01145 eNorm0 = eNorm;
01146 } else {
01147 CAL_nLayersRmsBack++;
01148 }
01149 eTot += eNorm;
01150 e2 += eNorm*eNorm;
01151 //std::cout << "layer " << layer << ", r.l. " << rlCsI[layer]
01152 // << ", eLayer " << CAL_eLayer[layer]
01153 // << ", rlNorm " << rlCsI[layer]*eCosTheta) << ", eL_norm "
01154 // << eNorm << std::endl;
01155 }
01156
01157 double eAve = 0;
01158 double eAveBack = 0;
01159 if(CAL_nLayersRms>1) {
01160 eAve = eTot/CAL_nLayersRms;
01161 CAL_rmsLayerE =
01162 (float)sqrt((e2 - CAL_nLayersRms*eAve*eAve)/(CAL_nLayersRms-1))/eAve;
01163 //std::cout << "eRms " << CAL_rmsLayerE/eAve
01164 // << ", " << CAL_nLayersRms << " layers" << std::endl;
01165 }
01166 if(CAL_nLayersRmsBack>1) {
01167 eAveBack = (eTot-eNorm0)/(CAL_nLayersRmsBack);
01168 CAL_rmsLayerEBack =
01169 (float) sqrt((e2 - eNorm0*eNorm0 - (CAL_nLayersRmsBack)*eAveBack*eAveBack)
01170 /(CAL_nLayersRmsBack-1))/eAve;
01171 //std::cout << "eRmsTrunc " << CAL_rmsLayerEBack/eAveBack
01172 // << ", " << nLayersRmsBack << " layers"<<std::endl;
01173 CAL_eAveBack = eAveBack;
01174 CAL_layer0Ratio = eNorm0/eAveBack;
01175 }
01176 }
01177 } catch( std::exception& /*e*/ ) {
01178 MsgStream log(msgSvc(), name());
01179 printHeader(log);
01180 setAnaTupBit();
01181 log << "See previous exception message." << endreq;
01182 log << " Skipping the Cal_rmsE calculation and resetting" << endreq;
01183
01184 CAL_layer0Ratio = _badFloat;
01185 CAL_eAveBack = _badFloat;
01186 CAL_nLayersRmsBack = _badInt;
01187 CAL_rmsLayerEBack = _badFloat;
01188 } catch(...) {
01189 MsgStream log(msgSvc(), name());
01190 printHeader(log);
01191 setAnaTupBit();
01192 log << "Unknown exception: see previous exception message, if any." << endreq;
01193 log << " Skipping the Cal_rmsE calculation" << endreq;
01194 log << "Initial track parameters: pos: " << xEnd << endreq
01195 << "dir: " << tEnd << " arclen: " << arclen << endreq;
01196
01197 CAL_layer0Ratio = _badFloat;
01198 CAL_eAveBack = _badFloat;
01199 CAL_nLayersRmsBack = _badInt;
01200 CAL_rmsLayerEBack = _badFloat;
01201 }
01202
01203 //
01204 // Estimations of the transverse size (Philippe Bruel)
01205 //
01206
01207 //
01208 // Perform the estimation with main axis = cal axis
01209 //
01210 TSaxisP = calCluster->getPosition();
01211 Vector TSaxisVin = calCluster->getDirection();
01212 TSaxisV = TSaxisVin.unit();
01213 //
01214 // Filling TSdist...
01215 //
01216 TSfillTSdist(pxtalrecs);
01217
01218 //int i;
01219
01220 if(TSnlog>0)
01221 {
01222 //
01223 // fill CAL_TS_CAL_T_
01224 //
01225 TSfillTS(0);
01226 CAL_TS_CAL_T_68 = (float)TSgetinterpolationTS(0.68);
01227 CAL_TS_CAL_T_90 = (float)TSgetinterpolationTS(0.90);
01228 CAL_TS_CAL_T_95 = (float)TSgetinterpolationTS(0.95);
01229 CAL_TS_CAL_T_99 = (float)TSgetinterpolationTS(0.99);
01230 CAL_TS_CAL_T_100 = (float)TSTS[TSnlog-1];
01231 //
01232 // fill CAL_TS_CAL_TL_
01233 //
01234 TSfillTS(1);
01235 CAL_TS_CAL_TL_68 = (float)TSgetinterpolationTS(0.68);
01236 CAL_TS_CAL_TL_90 = (float)TSgetinterpolationTS(0.90);
01237 CAL_TS_CAL_TL_95 = (float)TSgetinterpolationTS(0.95);
01238 CAL_TS_CAL_TL_99 = (float)TSgetinterpolationTS(0.99);
01239 CAL_TS_CAL_TL_100 = (float)TSTS[TSnlog-1];
01240 }
01241
01242 //
01243 // Perform the estimation with centroid and axis determined only with transverse information
01244 //
01245 if(CAL_posdir_nlayers>=4)
01246 {
01247 TSaxisP = Point(CAL_xEcntr2,CAL_yEcntr2,CAL_zEcntr2);
01248 TSaxisVin = Vector(CAL_xdir2,CAL_ydir2,CAL_zdir2);
01249 TSaxisV = TSaxisVin.unit();
01250 //
01251 // Filling TSdist...
01252 //
01253 TSfillTSdist(pxtalrecs);
01254
01255 if(TSnlog>0)
01256 {
01257 //
01258 // fill CAL_TS_CAL2_T_
01259 //
01260 TSfillTS(0);
01261 CAL_TS_CAL2_T_68 = (float)TSgetinterpolationTS(0.68);
01262 CAL_TS_CAL2_T_90 = (float)TSgetinterpolationTS(0.90);
01263 CAL_TS_CAL2_T_95 = (float)TSgetinterpolationTS(0.95);
01264 CAL_TS_CAL2_T_99 = (float)TSgetinterpolationTS(0.99);
01265 CAL_TS_CAL2_T_100 = (float)TSTS[TSnlog-1];
01266 //
01267 // fill CAL_TS_CAL2_TL_
01268 //
01269 TSfillTS(1);
01270 CAL_TS_CAL2_TL_68 = (float)TSgetinterpolationTS(0.68);
01271 CAL_TS_CAL2_TL_90 = (float)TSgetinterpolationTS(0.90);
01272 CAL_TS_CAL2_TL_95 = (float)TSgetinterpolationTS(0.95);
01273 CAL_TS_CAL2_TL_99 = (float)TSgetinterpolationTS(0.99);
01274 CAL_TS_CAL2_TL_100 = (float)TSTS[TSnlog-1];
01275 }
01276 }
01277
01278 //
01279 // Perform the estimation with main axis = best track
01280 //
01281 int mynumtracks = 0;
01282 if(pTracks) mynumtracks = pTracks->size();
01283 if(mynumtracks>0)
01284 {
01285 // Get the first track
01286 pTrack1 = pTracks->begin();
01287 track_1 = *pTrack1;
01288 // Get the start and direction
01289 TSaxisP = track_1->getInitialPosition();
01290 TSaxisVin = track_1->getInitialDirection();
01291 TSaxisV = TSaxisVin.unit();
01292 //
01293 // Filling TSdist...
01294 //
01295 TSfillTSdist(pxtalrecs);
01296
01297 if(TSnlog>0)
01298 {
01299 //
01300 // fill CAL_TS_TKR_T_
01301 //
01302 TSfillTS(0);
01303 CAL_TS_TKR_T_68 = (float)TSgetinterpolationTS(0.68);
01304 CAL_TS_TKR_T_90 = (float)TSgetinterpolationTS(0.90);
01305 CAL_TS_TKR_T_95 = (float)TSgetinterpolationTS(0.95);
01306 CAL_TS_TKR_T_99 = (float)TSgetinterpolationTS(0.99);
01307 CAL_TS_TKR_T_100 = (float)TSTS[TSnlog-1];
01308 //
01309 // fill CAL_TS_TKR_TL_
01310 //
01311 TSfillTS(1);
01312 CAL_TS_TKR_TL_68 = (float)TSgetinterpolationTS(0.68);
01313 CAL_TS_TKR_TL_90 = (float)TSgetinterpolationTS(0.90);
01314 CAL_TS_TKR_TL_95 = (float)TSgetinterpolationTS(0.95);
01315 CAL_TS_TKR_TL_99 = (float)TSgetinterpolationTS(0.99);
01316 CAL_TS_TKR_TL_100 = (float)TSTS[TSnlog-1];
01317 }
01318 }
01319
01320 // throw an exception
01321 //int ii = 1;
01322 //int j = 0;
01323 //int k = ii/j;
01324 //k++;
01325
01326 return sc;
01327 }
|
|
|
gets the CAL info from detModel
Definition at line 1329 of file CalValsTool.cxx. References m_calXWidth, m_calYWidth, m_calZBot, m_calZTop, and m_tkrGeom. Referenced by initialize().
01330 {
01331 m_calZTop = m_tkrGeom->calZTop();
01332 m_calZBot = m_tkrGeom->calZBot();
01333 m_calXWidth = m_tkrGeom->calXWidth();
01334 m_calYWidth = m_tkrGeom->calYWidth();
01335 /*
01336 double calModuleWid, calXtalLen;
01337 m_detSvc->getNumericConstByName("CalModuleWidth", &calModuleWid);
01338 m_detSvc->getNumericConstByName("CsILength", &calXtalLen);
01339 double m_deltaSide = calModuleWid - calXtalLen;
01340 */
01341
01342 return StatusCode::SUCCESS;
01343 }
|
|
|
Reimplemented from ValBase. Definition at line 471 of file CalValsTool.cxx. References ValBase::addItem(), CAL_BkHalf_Ratio, CAL_cfp_alpha, CAL_cfp_calEffRLn, CAL_cfp_calfit_alpha, CAL_cfp_calfit_calEffRLn, CAL_cfp_calfit_energy, CAL_cfp_calfit_fiterrflg, CAL_cfp_calfit_tmax, CAL_cfp_calfit_totChiSq, CAL_cfp_energy, CAL_cfp_fiterrflg, CAL_cfp_tkrRLn, CAL_cfp_tmax, CAL_cfp_totChiSq, CAL_Cntr_RLn, CAL_CsI_RLn, CAL_DeadCnt_Rat, CAL_DeadTot_Rat, CAL_deltaT, CAL_eAveBack, CAL_Edge_Corr, CAL_EdgeEnergy, CAL_eLayer, CAL_EnergyCorr, CAL_EnergyRaw, CAL_Gap_Fraction, CAL_LAT_RLn, CAL_LATEdge, CAL_layer0Ratio, CAL_Leak_Corr, CAL_LkHd_energy, CAL_LkHd_energyErr, CAL_lll_energy, CAL_lll_energyErr, CAL_Long_Rms, CAL_LRms_Asym, CAL_Lyr0_Ratio, CAL_Lyr7_Ratio, CAL_Max_Num_Xtals_In_Layer, CAL_MIP_Diff, CAL_MIP_Ratio, CAL_nLayersRmsBack, CAL_nsaturated, CAL_Num_Xtals, CAL_Num_Xtals_Trunc, CAL_posdir_chisq, CAL_posdir_nlayers, CAL_rmsLayerE, CAL_rmsLayerEBack, CAL_t_Pred, CAL_tkl_energy, CAL_tkl_energyErr, CAL_Top_Gap_Dist, CAL_Tot_RLn, CAL_Total_Corr, CAL_Track_Angle, CAL_Track_DOCA, CAL_track_E_rms, CAL_track_E_rms_trunc, CAL_track_rms, CAL_track_rms_trunc, CAL_Track_Sep, CAL_Trans_Rms, CAL_TS_CAL2_T_100, CAL_TS_CAL2_T_68, CAL_TS_CAL2_T_90, CAL_TS_CAL2_T_95, CAL_TS_CAL2_T_99, CAL_TS_CAL_T_100, CAL_TS_CAL_T_68, CAL_TS_CAL_T_90, CAL_TS_CAL_T_95, CAL_TS_CAL_T_99, CAL_TS_CAL_TL_100, CAL_TS_CAL_TL_68, CAL_TS_CAL_TL_90, CAL_TS_CAL_TL_95, CAL_TS_CAL_TL_99, CAL_TS_TKR_T_100, CAL_TS_TKR_T_68, CAL_TS_TKR_T_90, CAL_TS_TKR_T_95, CAL_TS_TKR_T_99, CAL_TS_TKR_TL_100, CAL_TS_TKR_TL_68, CAL_TS_TKR_TL_90, CAL_TS_TKR_TL_95, CAL_TS_TKR_TL_99, CAL_TwrEdgeCntr, CAL_TwrEdgeTop, CAL_x0, CAL_xdir, CAL_xdir2, CAL_xEcntr, CAL_xEcntr2, CAL_xPosRmsLastLayer, CAL_Xtal_maxEne, CAL_Xtal_Ratio, CAL_y0, CAL_ydir, CAL_ydir2, CAL_yEcntr, CAL_yEcntr2, CAL_yPosRmsLastLayer, CAL_zdir, CAL_zdir2, CAL_zEcntr, CAL_zEcntr2, getCalInfo(), ValBase::initialize(), m_detSvc, m_G4PropTool, m_nCsI, m_nLayers, m_tkrGeom, m_towerPitch, m_xNum, m_yNum, and zeroVals().
00472 {
00473 StatusCode sc = StatusCode::SUCCESS;
00474
00475 MsgStream log(msgSvc(), name());
00476
00477 if( ValBase::initialize().isFailure()) return StatusCode::FAILURE;
00478
00479 // get the services
00480
00481 if( serviceLocator() ) {
00482
00483 // find GlastDevSvc service
00484 if (service("GlastDetSvc", m_detSvc, true).isFailure()){
00485 log << MSG::ERROR << "Couldn't find the GlastDetSvc!" << endreq;
00486 return StatusCode::FAILURE;
00487 }
00488
00489 m_detSvc->getNumericConstByName("CALnLayer", &m_nLayers);
00490 m_detSvc->getNumericConstByName("nCsIPerLayer", &m_nCsI);
00491
00492 // find TkrGeometrySvc service
00493 if (service("TkrGeometrySvc", m_tkrGeom, true).isFailure()){
00494 log << MSG::ERROR << "Couldn't find the TkrGeometrySvc!" << endreq;
00495 return StatusCode::FAILURE;
00496 }
00497
00498 IToolSvc* toolSvc = 0;
00499 if(service("ToolSvc", toolSvc, true).isFailure()) {
00500 log << MSG::ERROR << "Couldn't find the ToolSvc!" << endreq;
00501 return StatusCode::FAILURE;
00502 }
00503
00504 if(!toolSvc->retrieveTool("G4PropagationTool", m_G4PropTool)) {
00505 log << MSG::ERROR << "Couldn't find the G4PropagationTool!" << endreq;
00506 return StatusCode::FAILURE;
00507 }
00508
00509
00510 m_towerPitch = m_tkrGeom->towerPitch();
00511 m_xNum = m_tkrGeom->numXTowers();
00512 m_yNum = m_tkrGeom->numYTowers();
00513
00514 if (getCalInfo().isFailure()) {
00515 log << MSG::ERROR << "Couldn't initialize the CAL constants" << endreq;
00516 return StatusCode::FAILURE;
00517 }
00518 } else {
00519 return StatusCode::FAILURE;
00520 }
00521
00522 // load up the map
00523
00524 addItem("CalEnergyRaw", &CAL_EnergyRaw);
00525 addItem("CalEnergyCorr", &CAL_EnergyCorr);
00526
00527 addItem("CalLeakCorr", &CAL_Leak_Corr);
00528 addItem("CalEdgeCorr", &CAL_Edge_Corr);
00529 addItem("CalTotalCorr", &CAL_Total_Corr);
00530
00531 addItem("CalCsIRLn", &CAL_CsI_RLn);
00532 addItem("CalTotRLn", &CAL_Tot_RLn);
00533 addItem("CalCntRLn", &CAL_Cntr_RLn);
00534 addItem("CalLATRLn", &CAL_LAT_RLn);
00535 addItem("CalDeadTotRat", &CAL_DeadTot_Rat);
00536 addItem("CalDeadCntRat", &CAL_DeadCnt_Rat);
00537
00538 addItem("CalTPred", &CAL_t_Pred);
00539 addItem("CalDeltaT", &CAL_deltaT);
00540
00541 addItem("CalGapFraction",&CAL_Gap_Fraction);
00542 addItem("CalTwrEdgeCntr",&CAL_TwrEdgeCntr);
00543 addItem("CalTwrEdge", &CAL_TwrEdgeTop);
00544 addItem("CalLATEdge", &CAL_LATEdge);
00545 addItem("CalEdgeEnergy", &CAL_EdgeEnergy);
00546 addItem("CalTrackDoca", &CAL_Track_DOCA);
00547 addItem("CalTrackAngle", &CAL_Track_Angle);
00548 addItem("CalTrackSep", &CAL_Track_Sep);
00549
00550 addItem("CalELayer0", &CAL_eLayer[0]);
00551 addItem("CalELayer1", &CAL_eLayer[1]);
00552 addItem("CalELayer2", &CAL_eLayer[2]);
00553 addItem("CalELayer3", &CAL_eLayer[3]);
00554 addItem("CalELayer4", &CAL_eLayer[4]);
00555 addItem("CalELayer5", &CAL_eLayer[5]);
00556 addItem("CalELayer6", &CAL_eLayer[6]);
00557 addItem("CalELayer7", &CAL_eLayer[7]);
00558 addItem("CalLyr0Ratio", &CAL_Lyr0_Ratio);
00559 addItem("CalLyr7Ratio", &CAL_Lyr7_Ratio);
00560 addItem("CalBkHalfRatio",&CAL_BkHalf_Ratio);
00561
00562 addItem("CalXtalsTrunc", &CAL_Num_Xtals_Trunc);
00563 addItem("CalNumXtals", &CAL_Num_Xtals);
00564 addItem("CalXtalRatio", &CAL_Xtal_Ratio);
00565 addItem("CalXtalMaxEne", &CAL_Xtal_maxEne);
00566 addItem("CalMaxNumXtalsInLayer",
00567 &CAL_Max_Num_Xtals_In_Layer);
00568
00569 addItem("CalTransRms", &CAL_Trans_Rms);
00570 addItem("CalLongRms", &CAL_Long_Rms);
00571 addItem("CalLRmsAsym", &CAL_LRms_Asym);
00572
00573 addItem("CalMIPDiff", &CAL_MIP_Diff);
00574 addItem("CalMIPRatio", &CAL_MIP_Ratio);
00575
00576 addItem("CalXEcntr", &CAL_xEcntr);
00577 addItem("CalYEcntr", &CAL_yEcntr);
00578 addItem("CalZEcntr", &CAL_zEcntr);
00579 addItem("CalXDir", &CAL_xdir);
00580 addItem("CalYDir", &CAL_ydir);
00581 addItem("CalZDir", &CAL_zdir);
00582 addItem("CalX0", &CAL_x0);
00583 addItem("CalY0", &CAL_y0);
00584 addItem("CalTopGapDist", &CAL_Top_Gap_Dist);
00585
00586 addItem("CalXEcntr2", &CAL_xEcntr2);
00587 addItem("CalYEcntr2", &CAL_yEcntr2);
00588 addItem("CalZEcntr2", &CAL_zEcntr2);
00589 addItem("CalXDir2", &CAL_xdir2);
00590 addItem("CalYDir2", &CAL_ydir2);
00591 addItem("CalZDir2", &CAL_zdir2);
00592 addItem("CalPosDirChisq", &CAL_posdir_chisq);
00593 addItem("CalPosDirNLayers", &CAL_posdir_nlayers);
00594 addItem("CalNSaturated", &CAL_nsaturated);
00595
00596 addItem("CalTrkXtalRms", &CAL_track_rms);
00597 addItem( |