Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2021-04-29 02:32:49 

0001 // Include our own header first
0002 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h"
0003 
0004 // Geometry services
0005 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0006 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0007 
0008 //#define DEBUG
0009 
0010 // MessageLogger
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 
0013 // Magnetic field
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015 
0016 // The template header files
0017 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco.h"
0018 
0019 // Commented for now (3/10/17) until we figure out how to resuscitate 2D template splitter
0020 /// #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
0021 
0022 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0023 
0024 #include <vector>
0025 #include "boost/multi_array.hpp"
0026 
0027 #include <iostream>
0028 
0029 using namespace SiPixelTemplateReco;
0030 //using namespace SiPixelTemplateSplit;
0031 using namespace std;
0032 
0033 namespace {
0034   constexpr float micronsToCm = 1.0e-4;
0035   constexpr int cluster_matrix_size_x = 13;
0036   constexpr int cluster_matrix_size_y = 21;
0037 }  // namespace
0038 
0039 //-----------------------------------------------------------------------------
0040 //  Constructor.
0041 //
0042 //-----------------------------------------------------------------------------
0043 PixelCPETemplateReco::PixelCPETemplateReco(edm::ParameterSet const& conf,
0044                                            const MagneticField* mag,
0045                                            const TrackerGeometry& geom,
0046                                            const TrackerTopology& ttopo,
0047                                            const SiPixelLorentzAngle* lorentzAngle,
0048                                            const SiPixelTemplateDBObject* templateDBobject)
0049     : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr, 1) {
0050   //cout << endl;
0051   //cout << "Constructing PixelCPETemplateReco::PixelCPETemplateReco(...)................................................." << endl;
0052   //cout << endl;
0053 
0054   // Configurable parameters
0055   //DoCosmics_ = conf.getParameter<bool>("DoCosmics"); // Not used in templates
0056   //LoadTemplatesFromDB_ = conf.getParameter<bool>("LoadTemplatesFromDB"); // Moved to Base
0057 
0058   //cout << " PixelCPETemplateReco : (int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
0059   //cout << "field_magnitude = " << field_magnitude << endl;
0060 
0061   // configuration parameter to decide between DB or text file template access
0062 
0063   if (LoadTemplatesFromDB_) {
0064     //cout << "PixelCPETemplateReco: Loading templates from database (DB) --------- " << endl;
0065 
0066     // Initialize template store to the selected ID [Morris, 6/25/08]
0067     if (!SiPixelTemplate::pushfile(*templateDBobject_, thePixelTemp_))
0068       throw cms::Exception("PixelCPETemplateReco")
0069           << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
0070           << (*templateDBobject_).version() << "\n\n";
0071   } else {
0072     //cout << "PixelCPETemplateReco : Loading templates for barrel and forward from ASCII files ----------" << endl;
0073     barrelTemplateID_ = conf.getParameter<int>("barrelTemplateID");
0074     forwardTemplateID_ = conf.getParameter<int>("forwardTemplateID");
0075     templateDir_ = conf.getParameter<int>("directoryWithTemplates");
0076 
0077     if (!SiPixelTemplate::pushfile(barrelTemplateID_, thePixelTemp_, templateDir_))
0078       throw cms::Exception("PixelCPETemplateReco")
0079           << "\nERROR: Template ID " << barrelTemplateID_
0080           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0081 
0082     if (!SiPixelTemplate::pushfile(forwardTemplateID_, thePixelTemp_, templateDir_))
0083       throw cms::Exception("PixelCPETemplateReco")
0084           << "\nERROR: Template ID " << forwardTemplateID_
0085           << " not loaded correctly from text file. Reconstruction will fail.\n\n";
0086   }
0087 
0088   speed_ = conf.getParameter<int>("speed");
0089   LogDebug("PixelCPETemplateReco::PixelCPETemplateReco:") << "Template speed = " << speed_ << "\n";
0090 
0091   UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
0092 }
0093 
0094 //-----------------------------------------------------------------------------
0095 //  Clean up.
0096 //-----------------------------------------------------------------------------
0097 PixelCPETemplateReco::~PixelCPETemplateReco() {}
0098 
0099 std::unique_ptr<PixelCPEBase::ClusterParam> PixelCPETemplateReco::createClusterParam(const SiPixelCluster& cl) const {
0100   return std::make_unique<ClusterParamTemplate>(cl);
0101 }
0102 
0103 //------------------------------------------------------------------
0104 //  Public methods mandated by the base class.
0105 //------------------------------------------------------------------
0106 
0107 //------------------------------------------------------------------
0108 //  The main call to the template code.
0109 //------------------------------------------------------------------
0110 LocalPoint PixelCPETemplateReco::localPosition(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0111   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0112 
0113   if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0114     throw cms::Exception("PixelCPETemplateReco::localPosition :") << "A non-pixel detector type in here?";
0115   //  barrel(false) or forward(true)
0116   const bool fpix = GeomDetEnumerators::isEndcap(theDetParam.thePart);
0117 
0118   int ID = -9999;
0119   if (LoadTemplatesFromDB_) {
0120     int ID0 = templateDBobject_->getTemplateID(theDetParam.theDet->geographicalId());  // just to comapre
0121     ID = theDetParam.detTemplateId;
0122     if (ID0 != ID)
0123       edm::LogError("PixelCPETemplateReco") << " different id" << ID << " " << ID0 << endl;
0124   } else {  // from asci file
0125     if (!fpix)
0126       ID = barrelTemplateID_;  // barrel
0127     else
0128       ID = forwardTemplateID_;  // forward
0129   }
0130   //cout << "PixelCPETemplateReco : ID = " << ID << endl;
0131 
0132   SiPixelTemplate templ(thePixelTemp_);
0133 
0134   // Preparing to retrieve ADC counts from the SiPixeltheClusterParam.theCluster->  In the cluster,
0135   // we have the following:
0136   //   int minPixelRow(); // Minimum pixel index in the x direction (low edge).
0137   //   int maxPixelRow(); // Maximum pixel index in the x direction (top edge).
0138   //   int minPixelCol(); // Minimum pixel index in the y direction (left edge).
0139   //   int maxPixelCol(); // Maximum pixel index in the y direction (right edge).
0140   // So the pixels from minPixelRow() will go into clust_array_2d[0][*],
0141   // and the pixels from minPixelCol() will go into clust_array_2d[*][0].
0142 
0143   int row_offset = theClusterParam.theCluster->minPixelRow();
0144   int col_offset = theClusterParam.theCluster->minPixelCol();
0145 
0146   // Store the coordinates of the center of the (0,0) pixel of the array that
0147   // gets passed to PixelTempReco1D
0148   // Will add these values to the output of  PixelTempReco1D
0149   float tmp_x = float(row_offset) + 0.5f;
0150   float tmp_y = float(col_offset) + 0.5f;
0151 
0152   // Store these offsets (to be added later) in a LocalPoint after tranforming
0153   // them from measurement units (pixel units) to local coordinates (cm)
0154   //
0155   //
0156 
0157   // In case of template reco failure, these are the lorentz drift corrections
0158   // to be applied
0159   float lorentzshiftX = 0.5f * theDetParam.lorentzShiftInCmX;
0160   float lorentzshiftY = 0.5f * theDetParam.lorentzShiftInCmY;
0161 
0162   // ggiurgiu@jhu.edu 12/09/2010 : update call with trk angles needed for bow/kink corrections
0163   LocalPoint lp;
0164 
0165   if (theClusterParam.with_track_angle)
0166     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred);
0167   else {
0168     edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0169                                           << "Should never be here. PixelCPETemplateReco should always be called with "
0170                                              "track angles. This is a bad error !!! ";
0171 
0172     lp = theDetParam.theTopol->localPosition(MeasurementPoint(tmp_x, tmp_y));
0173   }
0174 
0175   // first compute matrix size
0176   int mrow = 0, mcol = 0;
0177   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0178     auto pix = theClusterParam.theCluster->pixel(i);
0179     int irow = int(pix.x);
0180     int icol = int(pix.y);
0181     mrow = std::max(mrow, irow);
0182     mcol = std::max(mcol, icol);
0183   }
0184   mrow -= row_offset;
0185   mrow += 1;
0186   mrow = std::min(mrow, cluster_matrix_size_x);
0187   mcol -= col_offset;
0188   mcol += 1;
0189   mcol = std::min(mcol, cluster_matrix_size_y);
0190   assert(mrow > 0);
0191   assert(mcol > 0);
0192 
0193   float clustMatrix[mrow][mcol];
0194   memset(clustMatrix, 0, sizeof(float) * mrow * mcol);
0195 
0196   // Copy clust's pixels (calibrated in electrons) into clusMatrix;
0197   for (int i = 0; i != theClusterParam.theCluster->size(); ++i) {
0198     auto pix = theClusterParam.theCluster->pixel(i);
0199     int irow = int(pix.x) - row_offset;
0200     int icol = int(pix.y) - col_offset;
0201 
0202     // Gavril : what do we do here if the row/column is larger than cluster_matrix_size_x/cluster_matrix_size_y  ?
0203     // Ignore them for the moment...
0204     if ((irow < mrow) & (icol < mcol))
0205       clustMatrix[irow][icol] = float(pix.adc);
0206   }
0207 
0208   // Make and fill the bool arrays flagging double pixels
0209   bool xdouble[mrow], ydouble[mcol];
0210   // x directions (shorter), rows
0211   for (int irow = 0; irow < mrow; ++irow)
0212     xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX(irow + row_offset);
0213 
0214   // y directions (longer), columns
0215   for (int icol = 0; icol < mcol; ++icol)
0216     ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY(icol + col_offset);
0217 
0218   SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow, mcol};
0219 
0220   // Output:
0221   float nonsense = -99999.9f;  // nonsense init value
0222   theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ =
0223       theClusterParam.templSigmaY_ = nonsense;
0224   // If the template recontruction fails, we want to return 1.0 for now
0225   theClusterParam.templProbY_ = theClusterParam.templProbX_ = theClusterParam.templProbQ_ = 1.0f;
0226   theClusterParam.templQbin_ = 0;
0227   // We have a boolean denoting whether the reco failed or not
0228   theClusterParam.hasFilledProb_ = false;
0229 
0230   float templYrec1_ = nonsense;
0231   float templXrec1_ = nonsense;
0232   float templYrec2_ = nonsense;
0233   float templXrec2_ = nonsense;
0234 
0235   // ******************************************************************
0236   // Do it! Use cotalpha_ and cotbeta_ calculated in PixelCPEBase
0237 
0238   float locBz = theDetParam.bz;
0239   float locBx = theDetParam.bx;
0240 
0241   theClusterParam.ierr = PixelTempReco1D(ID,
0242                                          theClusterParam.cotalpha,
0243                                          theClusterParam.cotbeta,
0244                                          locBz,
0245                                          locBx,
0246                                          clusterPayload,
0247                                          templ,
0248                                          theClusterParam.templYrec_,
0249                                          theClusterParam.templSigmaY_,
0250                                          theClusterParam.templProbY_,
0251                                          theClusterParam.templXrec_,
0252                                          theClusterParam.templSigmaX_,
0253                                          theClusterParam.templProbX_,
0254                                          theClusterParam.templQbin_,
0255                                          speed_,
0256                                          theClusterParam.templProbQ_);
0257 
0258   // ******************************************************************
0259 
0260   // Check exit status
0261   if UNLIKELY (theClusterParam.ierr != 0) {
0262     LogDebug("PixelCPETemplateReco::localPosition")
0263         << "reconstruction failed with error " << theClusterParam.ierr << "\n";
0264 
0265     // Template reco has failed, compute position estimates based on cluster center of gravity + Lorentz drift
0266     // Future improvement would be to call generic reco instead
0267 
0268     // ggiurgiu@jhu.edu, 21/09/2010 : trk angles needed to correct for bows/kinks
0269     if (theClusterParam.with_track_angle) {
0270       theClusterParam.templXrec_ =
0271           theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0272       theClusterParam.templYrec_ =
0273           theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0274     } else {
0275       edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0276                                             << "Should never be here. PixelCPETemplateReco should always be called "
0277                                                "with track angles. This is a bad error !!! ";
0278 
0279       theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0280       theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0281     }
0282   } else if UNLIKELY (UseClusterSplitter_ && theClusterParam.templQbin_ == 0) {
0283     edm::LogError("PixelCPETemplateReco") << " PixelCPETemplateReco: Qbin = 0 but using cluster splitter, we should "
0284                                              "never be here !!!!!!!!!!!!!!!!!!!!!! \n"
0285                                           << "(int)UseClusterSplitter_ = " << (int)UseClusterSplitter_ << endl;
0286 
0287     //ierr =
0288     //PixelTempSplit( ID, fpix, cotalpha_, cotbeta_,
0289     //      clust_array_2d, ydouble, xdouble,
0290     //      templ,
0291     //      templYrec1_, templYrec2_, templSigmaY_, templProbY_,
0292     //      templXrec1_, templXrec2_, templSigmaX_, templProbX_,
0293     //      templQbin_ );
0294 
0295     // Commented for now (3/10/17) until we figure out how to resuscitate 2D template splitter
0296     ///      std::vector< SiPixelTemplateStore2D > thePixelTemp2D_;
0297     ///SiPixelTemplate2D::pushfile(ID, thePixelTemp2D_);
0298     ///SiPixelTemplate2D templ2D_(thePixelTemp2D_);
0299 
0300     theClusterParam.ierr = -123;
0301     /*
0302        float dchisq;
0303        float templProbQ_;
0304        SiPixelTemplateSplit::PixelTempSplit( ID, theClusterParam.cotalpha, theClusterParam.cotbeta,
0305        clust_array_2d,
0306        ydouble, xdouble,
0307        templ,
0308        templYrec1_, templYrec2_, theClusterParam.templSigmaY_, theClusterParam.templProbY_,
0309        templXrec1_, templXrec2_, theClusterParam.templSigmaX_, theClusterParam.templProbX_,
0310        theClusterParam.templQbin_,
0311        templProbQ_,
0312        true,
0313        dchisq,
0314        templ2D_ );
0315        
0316        */
0317     if (theClusterParam.ierr != 0) {
0318       // Template reco has failed, compute position estimates based on cluster center of gravity + Lorentz drift
0319       // Future improvement would be to call generic reco instead
0320 
0321       // ggiurgiu@jhu.edu, 12/09/2010 : trk angles needed to correct for bows/kinks
0322       if (theClusterParam.with_track_angle) {
0323         theClusterParam.templXrec_ =
0324             theDetParam.theTopol->localX(theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred) + lorentzshiftX;
0325         theClusterParam.templYrec_ =
0326             theDetParam.theTopol->localY(theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred) + lorentzshiftY;
0327       } else {
0328         edm::LogError("PixelCPETemplateReco") << "@SUB = PixelCPETemplateReco::localPosition"
0329                                               << "Should never be here. PixelCPETemplateReco should always be called "
0330                                                  "with track angles. This is a bad error !!! ";
0331         theClusterParam.templXrec_ = theDetParam.theTopol->localX(theClusterParam.theCluster->x()) + lorentzshiftX;
0332         theClusterParam.templYrec_ = theDetParam.theTopol->localY(theClusterParam.theCluster->y()) + lorentzshiftY;
0333       }
0334     } else {
0335       // go from micrometer to centimeter
0336       templXrec1_ *= micronsToCm;
0337       templYrec1_ *= micronsToCm;
0338       templXrec2_ *= micronsToCm;
0339       templYrec2_ *= micronsToCm;
0340 
0341       // go back to the module coordinate system
0342       templXrec1_ += lp.x();
0343       templYrec1_ += lp.y();
0344       templXrec2_ += lp.x();
0345       templYrec2_ += lp.y();
0346 
0347       // calculate distance from each hit to the track and choose the hit closest to the track
0348       float distX1 = std::abs(templXrec1_ - theClusterParam.trk_lp_x);
0349       float distX2 = std::abs(templXrec2_ - theClusterParam.trk_lp_x);
0350       float distY1 = std::abs(templYrec1_ - theClusterParam.trk_lp_y);
0351       float distY2 = std::abs(templYrec2_ - theClusterParam.trk_lp_y);
0352       theClusterParam.templXrec_ = (distX1 < distX2 ? templXrec1_ : templXrec2_);
0353       theClusterParam.templYrec_ = (distY1 < distY2 ? templYrec1_ : templYrec2_);
0354     }
0355   }  // else if ( UseClusterSplitter_ && templQbin_ == 0 )
0356 
0357   else  // apparenly this is the good one!
0358   {
0359     // go from micrometer to centimeter
0360     theClusterParam.templXrec_ *= micronsToCm;
0361     theClusterParam.templYrec_ *= micronsToCm;
0362 
0363     // go back to the module coordinate system
0364     theClusterParam.templXrec_ += lp.x();
0365     theClusterParam.templYrec_ += lp.y();
0366 
0367     // Compute the Alignment Group Corrections [template ID should already be selected from call to reco procedure]
0368     if (doLorentzFromAlignment_) {
0369       // Do only if the lotentzshift has meaningfull numbers
0370       if (theDetParam.lorentzShiftInCmX != 0.0 || theDetParam.lorentzShiftInCmY != 0.0) {
0371         // the LA width/shift returned by templates use (+)
0372         // the LA width/shift produced by PixelCPEBase for positive LA is (-)
0373         // correct this by iserting (-)
0374         //float temp1 = -micronsToCm*templ.lorxwidth();  // old
0375         //float temp2 = -micronsToCm*templ.lorywidth();  // does not incl 1/2
0376         float templateLorbiasCmX = -micronsToCm * templ.lorxbias();  // new
0377         float templateLorbiasCmY = -micronsToCm * templ.lorybias();  //incl. 1/2
0378         // now, correctly, we can use the difference of shifts
0379         //theClusterParam.templXrec_ += 0.5*(theDetParam.lorentzShiftInCmX - templateLorbiasCmX);
0380         //theClusterParam.templYrec_ += 0.5*(theDetParam.lorentzShiftInCmY - templateLorbiasCmY);
0381         theClusterParam.templXrec_ += (0.5 * (theDetParam.lorentzShiftInCmX) - templateLorbiasCmX);
0382         theClusterParam.templYrec_ += (0.5 * (theDetParam.lorentzShiftInCmY) - templateLorbiasCmY);
0383         //cout << "Templates: la lorentz offset = "
0384         //   <<(0.5*(theDetParam.lorentzShiftInCmX)-templateLorbiasCmX)
0385         //   <<" "<<templateLorbiasCmX<<" "<<templateLorbiasCmY
0386         //   <<" "<<temp1<<" "<<temp2
0387         //   <<" "<<theDetParam.lorentzShiftInCmX
0388         //   <<" "<<theDetParam.lorentzShiftInCmY
0389         //   << endl; //dk
0390       }  //else {cout<<" LA is 0, disable offset corrections "<<endl;} //dk
0391     }    //else {cout<<" Do not do LA offset correction "<<endl;} //dk
0392   }
0393 
0394   // Save probabilities and qBin in the quantities given to us by the base class
0395   // (for which there are also inline getters).  &&& templProbX_ etc. should be retired...
0396   theClusterParam.probabilityX_ = theClusterParam.templProbX_;
0397   theClusterParam.probabilityY_ = theClusterParam.templProbY_;
0398   theClusterParam.probabilityQ_ = theClusterParam.templProbQ_;
0399   theClusterParam.qBin_ = theClusterParam.templQbin_;
0400 
0401   if (theClusterParam.ierr == 0)  // always true here
0402     theClusterParam.hasFilledProb_ = true;
0403 
0404   return LocalPoint(theClusterParam.templXrec_, theClusterParam.templYrec_);
0405 }
0406 
0407 //------------------------------------------------------------------
0408 //  localError() relies on localPosition() being called FIRST!!!
0409 //------------------------------------------------------------------
0410 LocalError PixelCPETemplateReco::localError(DetParam const& theDetParam, ClusterParam& theClusterParamBase) const {
0411   ClusterParamTemplate& theClusterParam = static_cast<ClusterParamTemplate&>(theClusterParamBase);
0412 
0413   //cout << endl;
0414   //cout << "Set PixelCPETemplate errors .............................................." << endl;
0415 
0416   //cout << "CPETemplate : " << endl;
0417 
0418   //--- Default is the maximum error used for edge clusters.
0419   //--- (never used, in fact: let comment it out, shut up the complains of the static analyzer, and save a few CPU cycles)
0420   //   const float sig12 = 1./sqrt(12.0);
0421   //   float xerr = theDetParam.thePitchX *sig12;
0422   //   float yerr = theDetParam.thePitchY *sig12;
0423   float xerr, yerr;
0424 
0425   // Check if the errors were already set at the clusters splitting level
0426   if (theClusterParam.theCluster->getSplitClusterErrorX() > 0.0f &&
0427       theClusterParam.theCluster->getSplitClusterErrorX() < clusterSplitMaxError_ &&
0428       theClusterParam.theCluster->getSplitClusterErrorY() > 0.0f &&
0429       theClusterParam.theCluster->getSplitClusterErrorY() < clusterSplitMaxError_) {
0430     xerr = theClusterParam.theCluster->getSplitClusterErrorX() * micronsToCm;
0431     yerr = theClusterParam.theCluster->getSplitClusterErrorY() * micronsToCm;
0432 
0433     //cout << "Errors set at cluster splitting level : " << endl;
0434     //cout << "xerr = " << xerr << endl;
0435     //cout << "yerr = " << yerr << endl;
0436   } else {
0437     // If errors are not split at the cluster splitting level, set the errors here
0438 
0439     //cout  << "Errors are not split at the cluster splitting level, set the errors here : " << endl;
0440 
0441     int maxPixelCol = theClusterParam.theCluster->maxPixelCol();
0442     int maxPixelRow = theClusterParam.theCluster->maxPixelRow();
0443     int minPixelCol = theClusterParam.theCluster->minPixelCol();
0444     int minPixelRow = theClusterParam.theCluster->minPixelRow();
0445 
0446     //--- Are we near either of the edges?
0447     bool edgex = (theDetParam.theRecTopol->isItEdgePixelInX(minPixelRow) ||
0448                   theDetParam.theRecTopol->isItEdgePixelInX(maxPixelRow));
0449     bool edgey = (theDetParam.theRecTopol->isItEdgePixelInY(minPixelCol) ||
0450                   theDetParam.theRecTopol->isItEdgePixelInY(maxPixelCol));
0451 
0452     if (theClusterParam.ierr != 0) {
0453       // If reconstruction fails the hit position is calculated from cluster center of gravity
0454       // corrected in x by average Lorentz drift. Assign huge errors.
0455       //xerr = 10.0 * (float)theClusterParam.theCluster->sizeX() * xerr;
0456       //yerr = 10.0 * (float)theClusterParam.theCluster->sizeX() * yerr;
0457 
0458       if (!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
0459         throw cms::Exception("PixelCPETemplateReco::localPosition :") << "A non-pixel detector type in here?";
0460 
0461       // Assign better errors based on the residuals for failed template cases
0462       if (GeomDetEnumerators::isBarrel(theDetParam.thePart)) {
0463         xerr = 55.0f * micronsToCm;
0464         yerr = 36.0f * micronsToCm;
0465       } else {
0466         xerr = 42.0f * micronsToCm;
0467         yerr = 39.0f * micronsToCm;
0468       }
0469 
0470       //cout << "xerr = " << xerr << endl;
0471       //cout << "yerr = " << yerr << endl;
0472 
0473       //return LocalError(xerr*xerr, 0, yerr*yerr);
0474     } else if (edgex || edgey) {
0475       // for edge pixels assign errors according to observed residual RMS
0476       if (edgex && !edgey) {
0477         xerr = xEdgeXError_ * micronsToCm;
0478         yerr = xEdgeYError_ * micronsToCm;
0479       } else if (!edgex && edgey) {
0480         xerr = yEdgeXError_ * micronsToCm;
0481         yerr = yEdgeYError_ * micronsToCm;
0482       } else if (edgex && edgey) {
0483         xerr = bothEdgeXError_ * micronsToCm;
0484         yerr = bothEdgeYError_ * micronsToCm;
0485       } else {
0486         throw cms::Exception(" PixelCPETemplateReco::localError: Something wrong with pixel edge flag !!!");
0487       }
0488 
0489       //cout << "xerr = " << xerr << endl;
0490       //cout << "yerr = " << yerr << endl;
0491     } else {
0492       // &&& need a class const
0493       //const float micronsToCm = 1.0e-4;
0494 
0495       xerr = theClusterParam.templSigmaX_ * micronsToCm;
0496       yerr = theClusterParam.templSigmaY_ * micronsToCm;
0497 
0498       //cout << "xerr = " << xerr << endl;
0499       //cout << "yerr = " << yerr << endl;
0500 
0501       // &&& should also check ierr (saved as class variable) and return
0502       // &&& nonsense (another class static) if the template fit failed.
0503     }
0504 
0505     if (theVerboseLevel > 9) {
0506       LogDebug("PixelCPETemplateReco") << " Sizex = " << theClusterParam.theCluster->sizeX()
0507                                        << " Sizey = " << theClusterParam.theCluster->sizeY() << " Edgex = " << edgex
0508                                        << " Edgey = " << edgey << " ErrX  = " << xerr << " ErrY  = " << yerr;
0509     }
0510 
0511   }  // else
0512 
0513   if (!(xerr > 0.0f))
0514     throw cms::Exception("PixelCPETemplateReco::localError")
0515         << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
0516 
0517   if (!(yerr > 0.0f))
0518     throw cms::Exception("PixelCPETemplateReco::localError")
0519         << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
0520 
0521   //cout << "Final errors set to: " << endl;
0522   //cout << "xerr = " << xerr << endl;
0523   //cout << "yerr = " << yerr << endl;
0524   //cout << "Out of PixelCPETemplateREco..........................................................................." << endl;
0525   //cout << endl;
0526 
0527   return LocalError(xerr * xerr, 0, yerr * yerr);
0528 }
0529 
0530 void PixelCPETemplateReco::fillPSetDescription(edm::ParameterSetDescription& desc) {
0531   desc.add<int>("barrelTemplateID", 0);
0532   desc.add<int>("forwardTemplateID", 0);
0533   desc.add<int>("directoryWithTemplates", 0);
0534   desc.add<int>("speed", -2);
0535   desc.add<bool>("UseClusterSplitter", false);
0536 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!