RecNtupleAlg.cpp Source File

00001 // Implements ntuple writing algorithm
00002 
00003 #include "TkrRecon/RecNtupleAlg.h"
00004 
00005 #include "geometry/Ray.h"
00006 #include "geometry/Plane.h"
00007 #include "geometry/Intersection.h"
00008 #include "geometry/Point.h"
00009 #include "geometry/Vector.h"
00010 
00011 #include <algorithm>
00012 inline static double sqr(double x) {return x*x;}
00013 using namespace std;
00014 
00015 static const AlgFactory<RecNtupleAlg>  Factory;
00016 const IAlgFactory& RecNtupleAlgFactory = Factory;
00017 
00018 
00019 //###############################
00020 RecNtupleAlg::RecNtupleAlg(const std::string& name, ISvcLocator* pSvcLocator) :
00021               Algorithm(name, pSvcLocator)
00022 //###############################
00023 {
00024     declareProperty("tupleName", m_tupleName="");
00025 }
00026 
00027 //###############################
00028 StatusCode RecNtupleAlg::initialize() 
00029 //###############################
00030 {       
00031         StatusCode sc = StatusCode::SUCCESS;
00032 
00033     MsgStream log(msgSvc(), name());
00034     log << MSG::INFO << "initialize" << endreq;
00035 
00036     // Use the Job options service to set the Algorithm's parameters
00038     setProperties();
00039 
00040     // get a pointer to our ntupleWriterSvc
00041     sc = serviceLocator()->getService("ntupleWriterSvc", 
00042         IID_INTupleWriterSvc, reinterpret_cast<IInterface*&>( ntupleWriteSvc));
00043 
00044     if( sc.isFailure() ) 
00045     {
00046         log << MSG::ERROR << "RecNtupleAlg failed to get the ntupleWriterSvc" << endreq;
00047         return sc;
00048     }
00049 
00050     // get a pointer to our tracker geometry service
00051     sc = serviceLocator()->getService("TkrGeometrySvc", 
00052         IID_ITkrGeometrySvc, reinterpret_cast<IInterface*&>( pGeometry));
00053 
00054     if( sc.isFailure() ) 
00055     {
00056         log << MSG::ERROR << "RecNtupleAlg failed to get the TkrGeometrySvc" << endreq;
00057         return sc;
00058     }
00059 
00060         return sc;
00061 }
00062 //--------------------------------------------------------
00063 //     Execution Function fills the ntuple
00064 //--------------------------------------------------------
00065 //###############################
00066 StatusCode RecNtupleAlg::execute() 
00067 //###############################
00068 {       
00070     StatusCode sc = StatusCode::SUCCESS;
00071 
00072     //Create the tuple class
00073     RecTupleValues nTuple;
00074     
00075     //Retrieve pointers to the data stored in the TDS
00076     ICsIClusterList* pCalClusters = SmartDataPtr<ICsIClusterList>(eventSvc(),"/Event/CalRecon/CsIClusterList");
00077         SiClusters*      pSiClusters  = SmartDataPtr<SiClusters>(eventSvc(),"/Event/TkrRecon/SiClusters");
00078     SiRecObjs*       pSiRecObjs   = SmartDataPtr<SiRecObjs>(eventSvc(),"/Event/TkrRecon/SiRecObjs");
00079 
00080     if ((sc = nTuple.calcTupleValues(pCalClusters, pSiClusters, pSiRecObjs, pGeometry)).isFailure()) return sc;
00081 
00082     return nTuple.fillTupleValues(ntupleWriteSvc, m_tupleName.c_str());
00083 }
00084 //###############################
00085 StatusCode RecNtupleAlg::finalize() 
00086 //###############################
00087 {       
00088         StatusCode sc = StatusCode::SUCCESS;
00089 
00090     return sc;
00091 }
00092 
00093 //
00094 //Implementation of the RecTupleValues class begins here
00095 //
00096 
00097 //################################
00098 RecTupleValues::RecTupleValues()
00099 //################################
00100 {
00101     //Calculated in calcSkirtVars
00102     Rec_Tkr_SkirtX = -9999.;
00103     Rec_Tkr_SkirtY = -9999.;
00104 
00105     //Calculated in TowerBoundaries
00106     Rec_Conv_Twr_Dist = 0;
00107     Rec_Fit_xV = 0;
00108     Rec_Fit_yV = 0;
00109     //Calculated in ActiveDistance
00110     Rec_Active_Dist = 0;
00111     //Calculated in ExtraHits
00112     Rec_fst_Hit_Count = 0;
00113     Rec_Surplus_Hit_Ratio = 0;
00114     Rec_Outside_Hit_Ratio = 0;
00115     Rec_showerHits1 = 0;
00116     Rec_showerHits2 = 0;
00117     Rec_Sum_Hits = 0;
00118     //Calclulated in EnergyCorrection
00119     Rec_CsI_Energy = 0.3;   //Assumed energy by tracking fit
00120     Rec_CsI_Corr_Energy = 0;
00121 };
00122 
00123 //################################
00124 StatusCode RecTupleValues::calcTupleValues(ICsIClusterList* pCalClusters, SiClusters* pSiClusters, SiRecObjs* pRecObjs, ITkrGeometrySvc* pGeom)
00125 //################################
00126 {
00127         StatusCode sc = StatusCode::SUCCESS;
00128 
00129     //Make sure we have valid reconstructed data
00130     if (pRecObjs)
00131     {
00132         unsigned int  ngammas = pRecObjs->numGammas();
00133 
00134         if (ngammas == 0)   return sc;
00135 
00136         // Extract the total energy from the calorimeter
00137         if (pCalClusters)
00138         {
00139             ICsICluster* pCalClus = pCalClusters->Cluster(0);
00140             Rec_CsI_Energy       = pCalClus->energySum() / 1000; //GeV for now
00141             if (Rec_CsI_Energy < 0.001) Rec_CsI_Energy = 0.3;
00142         }
00143 
00144         // Right now we are assuming that the first gamma is the "right" gamma
00145         GFgamma* pGamma  = pRecObjs->Gamma(0);
00146 
00147         // If the gamma has no tracks then no point in going on 
00148         if (pGamma->empty()) return sc;
00149 
00150         // Call the routines to calculate the values
00151         calcSkirtVars(pGamma);
00152         calcTowerBoundaries(pGamma, pGeom);
00153         calcActiveDistance(pGamma, pGeom);
00154         calcExtraHits(pSiClusters, pGamma, pGeom);
00155 
00156         // This depends on the above having been called first !!
00157         calcEnergyCorrection(pGamma);
00158     }
00159     
00160     return sc;
00161 }
00162 
00163 //################################
00164 StatusCode RecTupleValues::fillTupleValues(INTupleWriterSvc* pSvc, const char* pName)
00165 //################################
00166 {
00167         StatusCode sc = StatusCode::SUCCESS;
00168 
00169     if (pSvc)
00170     {
00171         if ((sc = pSvc->addItem(pName, "REC_Tkr_SkirtX",        Rec_Tkr_SkirtX       )).isFailure()) return sc;
00172         if ((sc = pSvc->addItem(pName, "REC_Tkr_SkirtY",        Rec_Tkr_SkirtY       )).isFailure()) return sc;
00173         if ((sc = pSvc->addItem(pName, "REC_Conv_Twr_Dist",     Rec_Conv_Twr_Dist    )).isFailure()) return sc;
00174         if ((sc = pSvc->addItem(pName, "REC_Fit_xV",            Rec_Fit_xV           )).isFailure()) return sc;
00175         if ((sc = pSvc->addItem(pName, "REC_Fit_yV",            Rec_Fit_yV           )).isFailure()) return sc;
00176         if ((sc = pSvc->addItem(pName, "REC_Active_Dist",       Rec_Active_Dist      )).isFailure()) return sc;
00177         if ((sc = pSvc->addItem(pName, "REC_fst_Hit_Count",     Rec_fst_Hit_Count    )).isFailure()) return sc;
00178         if ((sc = pSvc->addItem(pName, "REC_Surplus_Hit_Ratio", Rec_Surplus_Hit_Ratio)).isFailure()) return sc;
00179         if ((sc = pSvc->addItem(pName, "REC_Outside_Hit_Ratio", Rec_Outside_Hit_Ratio)).isFailure()) return sc;
00180         if ((sc = pSvc->addItem(pName, "REC_showerHits1",       Rec_showerHits1      )).isFailure()) return sc;
00181         if ((sc = pSvc->addItem(pName, "REC_showerHits2",       Rec_showerHits2      )).isFailure()) return sc;
00182         if ((sc = pSvc->addItem(pName, "REC_Sum_Hits",          Rec_Sum_Hits         )).isFailure()) return sc;
00183         if ((sc = pSvc->addItem(pName, "REC_CsI_Corr_Energy",   Rec_CsI_Corr_Energy  )).isFailure()) return sc;
00184     }
00185     
00186     return sc;
00187 }
00188 
00189 
00190 //########################################################
00191 void RecTupleValues::calcSkirtVars(GFgamma* pGamma)
00192 //########################################################
00193 { 
00194     GFtrack* pBestX  = pGamma->getXpair()->getBest();
00195     GFtrack* pBestY  = pGamma->getYpair()->getBest();
00196 
00197     double lastHitX  = 0.;
00198     double lastHitXZ = 0.;
00199     double lastHitY  = 0.;
00200     double lastHitYZ = 0.;
00201 
00202     double lastSlpX  = 1.;
00203     double lastSlpY  = 1.;
00204 
00205     //Extract the position and slope of the last track hit in the X projection
00206     if (!pBestX->empty())
00207     {
00208         int    nHits   = pBestX->numDataPoints();
00209         KalHit lastHit = pBestX->kplanelist[nHits-1].getHit(KalHit::FIT);
00210 
00211         lastHitX  = lastHit.getPar().getPosition();
00212         lastSlpX  = lastHit.getPar().getSlope();
00213         lastHitXZ = pBestX->kplanelist[nHits-1].getZPlane();
00214     }
00215 
00216     //Extract the position and slope of the last track hit in the Y projection
00217     if (!pBestY->empty())
00218     {
00219         int    nHits   = pBestY->numDataPoints();
00220         KalHit lastHit = pBestY->kplanelist[nHits-1].getHit(KalHit::FIT);
00221 
00222         lastHitY  = lastHit.getPar().getPosition();
00223         lastSlpY  = lastHit.getPar().getSlope();
00224         lastHitYZ = pBestY->kplanelist[nHits-1].getZPlane();
00225     }
00226 
00227     //Determine a plane for the gap in the skirt
00228     Point        skirt(0., 0., -14.8751);   
00229     Vector       normalVec(0., 0., 1.);
00230     Plane        p(skirt, normalVec);
00231 
00232     //Find the intersection of the X track with this plane
00233     Vector       dir = Vector(lastSlpX, 0., 1.).unit();
00234     Point        lastHit(lastHitX, 0., lastHitXZ);
00235     Ray          r(lastHit, dir);
00236     Intersection intersectX(r, p);
00237 
00238     //The call to get the intersection is:
00239     // Intersection::distance(double maxStep, int sign)
00240     //where maxStep is the maximum distance you want the intersecting alg to search
00241     //before giving up on finding an intersection
00242     //If sign + the leaving solutions are returned, if sign - entering sol.
00243     double       dist = intersectX.distance(100., 1);
00244 
00245     if (dist < FLT_MAX)
00246     {
00247         Point pos = r.position(dist);
00248         Rec_Tkr_SkirtX = pos.x();
00249     }
00250 
00251     //Repeat for the Y track 
00252     dir       = Vector(0., lastSlpY, 1.).unit();
00253     lastHit   = Point(0., lastHitY, lastHitYZ);
00254     r         = Ray(lastHit, dir);
00255     Intersection intersectY(r, p);
00256 
00257     //The call to get the intersection is:
00258     // Intersection::distance(double maxStep, int sign)
00259     //where maxStep is the maximum distance you want the intersecting alg to search
00260     //before giving up on finding an intersection
00261     //If sign + the leaving solutions are returned, if sign - entering sol.
00262     dist = intersectY.distance(100., 1);
00263 
00264 
00265     if (dist < FLT_MAX)
00266     {
00267         Point pos = r.position(dist);
00268         Rec_Tkr_SkirtY = pos.y();
00269     }
00270 
00271     return;    
00272 }
00273     
00274 //########################################################
00275 void RecTupleValues::calcTowerBoundaries(GFgamma* pGamma, ITkrGeometrySvc* pGeom)
00276 //########################################################
00277 {
00278     Point vertex = pGamma->getFirstHit();
00279     double x_ist = vertex.x();
00280     double y_ist = vertex.y();
00281 
00282     double tower_width = pGeom->towerPitch();
00283     double num_2       = pGeom->numXTowers()/2.;
00284     
00285     Rec_Fit_xV = fmod((x_ist + tower_width*num_2), tower_width);
00286     Rec_Fit_xV = Rec_Fit_xV - tower_width*.5;
00287     Rec_Fit_yV = fmod((y_ist + tower_width*num_2), tower_width);
00288     Rec_Fit_yV = Rec_Fit_yV - tower_width*.5;
00289     
00290     Rec_Fit_xV = fabs(Rec_Fit_xV);
00291     Rec_Fit_yV = fabs(Rec_Fit_yV);
00292 
00293     Rec_Conv_Twr_Dist   = (Rec_Fit_xV > Rec_Fit_yV) ? Rec_Fit_xV : Rec_Fit_yV;
00294     
00295     return;
00296 }
00297 
00298 
00299 //########################################################
00300 void RecTupleValues::calcActiveDistance(GFgamma* pGamma, ITkrGeometrySvc* pGeom)
00301 //########################################################
00302 {  
00303     Rec_Active_Dist  = -20.;    
00304 
00305 /* //THB: This depends on simulation environment, functionality needs to be restored
00306     Point firstHit = pGamma->getFirstHit();
00307     Ray prjRay(firstHit, -t0);
00308     Point x_prj;
00309     
00310     GlastMed *glast = GlastMed::instance();
00311     
00312     int istLayer = pGamma->firstLayer();
00313     
00314     float firstProp = 0.;
00315     if (istLayer > 0) { //go to next gap up + espsilon
00316         firstProp = pGeom->trayHeight()/fabs(t0.z());
00317         x_prj = prjRay.position(firstProp);
00318     }
00319     else {
00320         x_prj = firstHit;
00321         firstProp = 0.;
00322     }
00323     
00324     const Medium *med0 = glast->inside(x_prj);
00325     prjRay = Ray(x_prj, -t0);
00326     const Medium *detMed = 0;
00327     if(med0){
00328         firstProp = med0->distanceToLeave(prjRay, detMed, 100.);
00329         activeDist = -2. - firstProp;
00330     }
00331     
00332     Detector *det =  0;
00333     if(detMed) {
00334         med0 = detMed;
00335         x_prj = prjRay.position(firstProp);
00336         prjRay = Ray(x_prj, -t0);
00337         firstProp = med0->distanceToLeave(prjRay, detMed, 100.);
00338         activeDist = -4. - firstProp;
00339         if(detMed) det = detMed->detector();
00340         if(det) {
00341             if(!strcmp(det->nameOf(), "SiDetector")) {
00342                 x_prj = prjRay.position(firstProp);
00343                 activeDist = static_cast<const SiDetectorMed*>(det)->data().insideActiveArea(x_prj);
00344             }
00345         }
00346     }
00347 */    
00348     return;
00349 }
00350 
00351 
00352 //########################################################
00353 void RecTupleValues::calcExtraHits(SiClusters* pSiClusters, GFgamma* pGamma, ITkrGeometrySvc* pGeom)
00354 //########################################################
00355 {  
00356     // Initialization    
00357     double norma = 1.;
00358 
00359     // Set to the total energy in the calorimeter
00360     double CsICorrEnergy = Rec_CsI_Energy;
00361 
00362     // Retrieve gamma vertex and direction
00363     Point  x0         = pGamma->vertex();
00364     Vector t0         = pGamma->direction();
00365 
00366     Point  x1         = GFdata::doVertex(pGamma->getBest(SiCluster::X)->ray(),
00367                                          pGamma->getBest(SiCluster::Y)->ray());
00368     Vector t1         = GFdata::doDirection(pGamma->getBest(SiCluster::X)->direction(),
00369                                             pGamma->getBest(SiCluster::Y)->direction());
00370 
00371         //      Find number of hits around first hit: 5 sigma out
00372     int    firstLayer = pGamma->firstLayer();
00373     Point  firstHit   = pGamma->vertex();
00374 
00375     int    diffXY     = pGamma->getXpair()->firstLayer()-pGamma->getYpair()->firstLayer();
00376     double dz         = pGamma->getBest(SiCluster::X)->vertex().z()  
00377                       - pGamma->getBest(SiCluster::Y)->vertex().z();
00378     
00379     double hitRadFact = fabs(t0.z()) + sqrt(1. - t0.z()*t0.z())/fabs(t0.z());
00380 
00381     if(hitRadFact > 3.) hitRadFact = 3.;
00382 
00383     int    ipln       = pGeom->numPlanes() - firstLayer; 
00384     double radLenEff  = (pGeom->pbRadLen(max(0,ipln-1)) + 2.*pGeom->siThickness()/9.36 +.003)/fabs(t1.z());
00385     double thetaCone  = .014/(CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
00386     double minSprd    = 5.* thetaCone * pGeom->trayHeight(); // 5 sigma cut
00387     double dltSprd    = minSprd;
00388     double norm       = sqrt(t0.x()*t0.x() + t0.y()*t0.y());
00389     double xFact      = 1., 
00390            yFact      = 1.;
00391 
00392     if(norm > 10000.*FLT_MIN) {
00393         xFact = 1. + (hitRadFact - 1.)*fabs(t0.x())/norm;
00394         yFact = 1. + (hitRadFact - 1.)*fabs(t0.y())/norm;
00395     }
00396 
00397     double dltX       = dltSprd*xFact;
00398     double dltY       = dltSprd*yFact;
00399     
00400     // Compute the First_hit_count variable
00401     Rec_fst_Hit_Count = pSiClusters->numberOfHitsNear(firstLayer, dltX, dltY, firstHit);
00402 
00403     // separate the special case of 3 into 3 or 2.5
00404     if(Rec_fst_Hit_Count == 3 && diffXY == 0) { // Check if extra hits on the far side of the first gap
00405         SiCluster::view far_axis = SiCluster::X;
00406         double dR = dltX;
00407         if (dz>0) {
00408             far_axis = SiCluster::Y;
00409             dR = dltY;
00410         }
00411         int ihit_second = pSiClusters->numberOfHitsNear(far_axis, firstLayer, dltX, firstHit);
00412         if (ihit_second == 2) Rec_fst_Hit_Count -=.5;
00413     }
00414     if(Rec_fst_Hit_Count == 1 && diffXY == 0) { // Check if extra hits on the far side of the first gap
00415         SiCluster::view far_axis = SiCluster::X;
00416         double dR = dltX;
00417         if (dz>0) {
00418             far_axis = SiCluster::Y;
00419             dR = dltY;
00420         }
00421         int ihit_second = pSiClusters->numberOfHitsNear(far_axis, firstLayer, dltX, firstHit);
00422         if (ihit_second == 1) Rec_fst_Hit_Count -=.5;
00423     }
00424     
00425     //  Find the number of hits around the rest of the gamma trajectory
00426     double deltaS    = pGeom->trayHeight()/fabs(t0.z());
00427     double disp      =  deltaS;
00428     int    lastLayer = (int)(4+2.*log(CsICorrEnergy/.01) + firstLayer);
00429     if(lastLayer > pGeom->numPlanes()) lastLayer = pGeom->numPlanes();
00430     
00431     Rec_Sum_Hits = pSiClusters->numberOfHitsNear(firstLayer, .5*dltX, .5*dltY, firstHit);
00432 
00433     Rec_showerHits1 = Rec_Sum_Hits;
00434     Rec_showerHits2 = 0;
00435     
00436     if(firstLayer > 11) {
00437         Rec_showerHits1 = 0; 
00438         Rec_showerHits2 = Rec_Sum_Hits;
00439     }
00440     
00441     if(Rec_Sum_Hits < 2.) Rec_Sum_Hits = 2.;
00442     if(lastLayer - firstLayer < 5 && lastLayer == pGeom->numPlanes()) Rec_Sum_Hits +=.5;
00443     double outHits = 0;
00444     
00445     xFact *= sqrt(t0.z() *t0.z() + t0.x()*t0.x());
00446     yFact *= sqrt(t0.z() *t0.z() + t0.y()*t0.y());
00447     
00448     int i;
00449     for(i=firstLayer+1; i<pGeom->numPlanes(); i++) {
00450         ipln = pGeom->numPlanes() - i;
00451         radLenEff  = (pGeom->pbRadLen(max(0,ipln-1)) + 2.*pGeom->siThickness()/9.36 +.003)/fabs(t1.z());
00452         double thetaMS = .014/(CsICorrEnergy/2.)*sqrt(radLenEff)*(1.+.038*log(radLenEff));
00453         double s_total = (i-firstLayer) * deltaS;
00454         double xSprd = thetaMS * s_total * xFact * 2.021; // = 3.5 sigma/sqrt(3)
00455         double ySprd = thetaMS * s_total * yFact * 2.021;
00456         Point trkHit((Point)(disp*t0 + x0));
00457         if (xSprd > 50) xSprd = 50.;
00458         if (ySprd > 50) ySprd = 50.;
00459         double nearHits = pSiClusters->numberOfHitsNear(i, xSprd, ySprd, trkHit);
00460 
00461         if(i< 12) Rec_showerHits1 +=  nearHits;
00462         else      Rec_showerHits2 +=  nearHits;
00463         
00464         if( i < lastLayer) {
00465             Rec_Sum_Hits +=     nearHits;
00466             outHits +=  pSiClusters->numberOfHitsNear(i, 50.,   50.,   trkHit) - nearHits;
00467         }
00468         disp += deltaS;
00469     }
00470     
00471     norma = lastLayer - firstLayer;
00472 
00473     Rec_Surplus_Hit_Ratio = Rec_Sum_Hits / norma;
00474     Rec_Outside_Hit_Ratio = outHits / norma;
00475 
00476     return;
00477 }
00478 
00479 
00480 //########################################################
00481 void RecTupleValues::calcEnergyCorrection(GFgamma* pGamma)
00482 //########################################################
00483 { 
00484     Rec_CsI_Corr_Energy = Rec_CsI_Energy;
00485 
00486     if (Rec_Conv_Twr_Dist > 0)
00487     {
00488         GFtrack* _Xbest = pGamma->getXpair()->getBest();
00489         GFtrack* _Ybest = pGamma->getYpair()->getBest();
00490         Vector    t0    = pGamma->direction();
00491 
00492 
00493         double hit_energy = 0.0003*Rec_showerHits1      + .0013*Rec_showerHits2;  //Energy in Hits
00494         double num_Layers_Crossed = _Xbest->numDataPoints()+ _Xbest->numGaps();
00495 
00496         if(_Xbest->numDataPoints() < _Ybest->numDataPoints()) 
00497          num_Layers_Crossed = _Ybest->numDataPoints()+ _Ybest->numGaps();
00498 
00499         double x_ing_energy = .0016*num_Layers_Crossed; // Take out observed layer correlation
00500         double edg_corr   =     (Rec_Conv_Twr_Dist < 10.) ? 1 : .68 + (.069  - .0036* Rec_Conv_Twr_Dist )*Rec_Conv_Twr_Dist; // Edge of Tower correction 
00501         double trk_energy = hit_energy + x_ing_energy;  // Total Tracker energy 
00502     
00503         trk_energy /= (fabs(t0.z()) > .2) ? fabs(t0.z()) : .2; // Note: t0.z < 0
00504         Rec_CsI_Corr_Energy = (Rec_CsI_Energy + trk_energy)/edg_corr;
00505     }
00506 
00507     return;    
00508 }
00509 
00510 
00511