CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoLocalTracker/SiPixelRecHits/src/PixelCPETemplateReco.cc

Go to the documentation of this file.
00001 // Include our own header first
00002 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPETemplateReco.h"
00003 
00004 // Geometry services
00005 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00006 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
00007 
00008 //#define DEBUG
00009 
00010 // MessageLogger
00011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00012 
00013 // Magnetic field
00014 #include "MagneticField/Engine/interface/MagneticField.h"
00015 
00016 // The template header files
00017 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateReco.h"
00018 
00019 #include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplateSplit.h"
00020 
00021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00022 
00023 #include <vector>
00024 #include "boost/multi_array.hpp"
00025 
00026 #include <iostream>
00027 
00028 using namespace SiPixelTemplateReco;
00029 //using namespace SiPixelTemplateSplit;
00030 using namespace std;
00031 
00032 const float PI = 3.141593;
00033 const float HALFPI = PI * 0.5;
00034 const float degsPerRad = 57.29578;
00035 
00036 // &&& need a class const
00037 const float micronsToCm = 1.0e-4;
00038 
00039 const int cluster_matrix_size_x = 13;
00040 const int cluster_matrix_size_y = 21;
00041 
00042 //-----------------------------------------------------------------------------
00043 //  Constructor.  All detUnit-dependent quantities will be initialized later,
00044 //  in setTheDet().  Here we only load the templates into the template store templ_ .
00045 //-----------------------------------------------------------------------------
00046 PixelCPETemplateReco::PixelCPETemplateReco(edm::ParameterSet const & conf, 
00047                                            const MagneticField * mag, const SiPixelLorentzAngle * lorentzAngle, 
00048                                            const SiPixelTemplateDBObject * templateDBobject) 
00049   : PixelCPEBase(conf, mag, lorentzAngle, 0, templateDBobject)
00050 {
00051   //cout << endl;
00052   //cout << "Constructing PixelCPETemplateReco::PixelCPETemplateReco(...)................................................." << endl;
00053   //cout << endl;
00054 
00055   // &&& initialize the templates, etc.
00056   
00057   //-- Use Magnetic field at (0,0,0) to select a template ID [Morris, 6/25/08] (temporary until we implement DB access)
00058   
00059   DoCosmics_ = conf.getParameter<bool>("DoCosmics");
00060   LoadTemplatesFromDB_ = conf.getParameter<bool>("LoadTemplatesFromDB");
00061   //cout << " PixelCPETemplateReco : (int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
00062   //cout << "field_magnitude = " << field_magnitude << endl;
00063   
00064   // ggiurgiu@fnal.gov, 12/17/2008: use configuration parameter to decide between DB or text file template access
00065   if ( LoadTemplatesFromDB_ )
00066     {
00067       //cout << "PixelCPETemplateReco: Loading templates from database (DB) ------------------------------- " << endl;
00068       
00069       // Initialize template store to the selected ID [Morris, 6/25/08]  
00070       if ( !templ_.pushfile( *templateDBobject_) )
00071         throw cms::Exception("PixelCPETemplateReco") 
00072           << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version " 
00073           << (*templateDBobject_).version() << "\n\n";
00074     }
00075   else 
00076     {
00077       //cout << "PixelCPETemplateReco : Loading templates 40 and 41 from ASCII files ------------------------" << endl;
00078 
00079       if ( !templ_.pushfile( 40 ) )
00080         throw cms::Exception("PixelCPETemplateReco") 
00081           << "\nERROR: Templates 40 not loaded correctly from text file. Reconstruction will fail.\n\n";
00082     
00083       if ( !templ_.pushfile( 41 ) )
00084         throw cms::Exception("PixelCPETemplateReco") 
00085           << "\nERROR: Templates 41 not loaded correctly from text file. Reconstruction will fail.\n\n";
00086     }
00087   
00088   speed_ = conf.getParameter<int>( "speed");
00089   LogDebug("PixelCPETemplateReco::PixelCPETemplateReco:") <<
00090     "Template speed = " << speed_ << "\n";
00091   
00092   UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
00093   
00094 }
00095 
00096 //-----------------------------------------------------------------------------
00097 //  Clean up.
00098 //-----------------------------------------------------------------------------
00099 PixelCPETemplateReco::~PixelCPETemplateReco()
00100 {
00101   // &&& delete template store?
00102 }
00103 
00104 
00105 MeasurementPoint 
00106 PixelCPETemplateReco::measurementPosition(const SiPixelCluster& cluster, 
00107                                           const GeomDetUnit & det) const
00108 {
00109   LocalPoint lp = localPosition(cluster,det);
00110 
00111   // ggiurgiu@jhu.edu 12/09/2010 : trk angles needed to correct for bows/kinks
00112   if ( with_track_angle )
00113     return theTopol->measurementPosition(lp, Topology::LocalTrackAngles( loc_traj_param_.dxdz(), loc_traj_param_.dydz() ) );
00114   else
00115     {
00116       edm::LogError("PixelCPETemplateReco") 
00117         << "@SUB = PixelCPETemplateReco::measurementPosition" 
00118         << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
00119       
00120       return theTopol->measurementPosition( lp ); 
00121     }
00122 
00123 }
00124 
00125 
00126 //------------------------------------------------------------------
00127 //  Public methods mandated by the base class.
00128 //------------------------------------------------------------------
00129 
00130 //------------------------------------------------------------------
00131 //  The main call to the template code.
00132 //------------------------------------------------------------------
00133 LocalPoint
00134 PixelCPETemplateReco::localPosition(const SiPixelCluster& cluster, const GeomDetUnit & det) const 
00135 {
00136   setTheDet( det, cluster );
00137   
00138   bool fpix;  //  barrel(false) or forward(true)
00139   if ( thePart == GeomDetEnumerators::PixelBarrel )   
00140     fpix = false;    // no, it's not forward -- it's barrel
00141   else                                              
00142     fpix = true;     // yes, it's forward
00143   
00144   int ID = -9999;
00145 
00146   if ( LoadTemplatesFromDB_ )
00147     {
00148       ID = templateDBobject_->getTemplateID(theDet->geographicalId());
00149     }
00150   else
00151     {
00152       if ( !fpix )
00153         ID = 40; // barrel
00154       else 
00155         ID = 41; // endcap
00156     }
00157   
00158   //cout << "PixelCPETemplateReco : ID = " << ID << endl;
00159   
00160 
00161   // Make from cluster (a SiPixelCluster) a boost multi_array_2d called 
00162   // clust_array_2d.
00163   boost::multi_array<float, 2> clust_array_2d(boost::extents[cluster_matrix_size_x][cluster_matrix_size_y]);
00164   
00165   // Preparing to retrieve ADC counts from the SiPixelCluster.  In the cluster,
00166   // we have the following:
00167   //   int minPixelRow(); // Minimum pixel index in the x direction (low edge).
00168   //   int maxPixelRow(); // Maximum pixel index in the x direction (top edge).
00169   //   int minPixelCol(); // Minimum pixel index in the y direction (left edge).
00170   //   int maxPixelCol(); // Maximum pixel index in the y direction (right edge).
00171   // So the pixels from minPixelRow() will go into clust_array_2d[0][*],
00172   // and the pixels from minPixelCol() will go into clust_array_2d[*][0].
00173   int row_offset = cluster.minPixelRow();
00174   int col_offset = cluster.minPixelCol();
00175   
00176   // Store the coordinates of the center of the (0,0) pixel of the array that 
00177   // gets passed to PixelTempReco2D
00178   // Will add these values to the output of  PixelTempReco2D
00179   float tmp_x = float(cluster.minPixelRow()) + 0.5;
00180   float tmp_y = float(cluster.minPixelCol()) + 0.5;
00181   
00182   // Store these offsets (to be added later) in a LocalPoint after tranforming 
00183   // them from measurement units (pixel units) to local coordinates (cm)
00184 
00185   // ggiurgiu@jhu.edu 12/09/2010 : update call with trk angles needed for bow/kink corrections
00186   LocalPoint lp;
00187   
00188   if ( with_track_angle )
00189     lp = theTopol->localPosition( MeasurementPoint(tmp_x, tmp_y), loc_trk_pred_ );
00190   else
00191     {
00192       edm::LogError("PixelCPETemplateReco") 
00193         << "@SUB = PixelCPETemplateReco::localPosition" 
00194         << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
00195       
00196       lp = theTopol->localPosition( MeasurementPoint(tmp_x, tmp_y) );
00197     }
00198     
00199   const std::vector<SiPixelCluster::Pixel> & pixVec = cluster.pixels();
00200   std::vector<SiPixelCluster::Pixel>::const_iterator 
00201     pixIter = pixVec.begin(), pixEnd = pixVec.end();
00202   
00203   // Visualize large clusters ---------------------------------------------------------
00204   // From Petar: maybe this should be moved into a method in the base class?
00205   /*
00206     char cluster_matrix[100][100];
00207     for (int i=0; i<100; i++)
00208     for (int j=0; j<100; j++)
00209     cluster_matrix[i][j] = '.';
00210     
00211     if ( cluster.sizeX()>cluster_matrix_size_x || cluster.sizeY()>cluster_matrix_size_y )
00212     {           
00213     cout << "cluster.size()  = " << cluster.size()  << endl;
00214     cout << "cluster.sizeX() = " << cluster.sizeX() << endl;
00215     cout << "cluster.sizeY() = " << cluster.sizeY() << endl;
00216     
00217     for ( std::vector<SiPixelCluster::Pixel>::const_iterator pix = pixVec.begin(); pix != pixVec.end(); ++pix )
00218     {
00219     int i = (int)(pix->x) - row_offset;
00220     int j = (int)(pix->y) - col_offset;
00221     cluster_matrix[i][j] = '*';
00222     }
00223     
00224     for (int i=0; i<(int)cluster.sizeX()+2; i++)
00225     {
00226     for (int j=0; j<(int)cluster.sizeY()+2; j++)
00227     cout << cluster_matrix[i][j];
00228     cout << endl;
00229     }
00230     } // if ( cluster.sizeX()>cluster_matrix_size_x || cluster.sizeY()>cluster_matrix_size_y )
00231   */
00232   // End Visualize clusters ---------------------------------------------------------
00233   
00234 
00235   // Copy clust's pixels (calibrated in electrons) into clust_array_2d;
00236   for ( ; pixIter != pixEnd; ++pixIter ) 
00237     {
00238       // *pixIter dereferences to Pixel struct, with public vars x, y, adc (all float)
00239       // 02/13/2008 ggiurgiu@fnal.gov: type of x, y and adc has been changed to unsigned char, unsigned short, unsigned short
00240       // in DataFormats/SiPixelCluster/interface/SiPixelCluster.h so the type cast to int is redundant. Leave it there, it 
00241       // won't hurt. 
00242       int irow = int(pixIter->x) - row_offset;   // &&& do we need +0.5 ???
00243       int icol = int(pixIter->y) - col_offset;   // &&& do we need +0.5 ???
00244       
00245       // Gavril : what do we do here if the row/column is larger than cluster_matrix_size_x/cluster_matrix_size_y = 7/21 ?
00246       // Ignore them for the moment...
00247       if ( irow<cluster_matrix_size_x && icol<cluster_matrix_size_y )
00248         // 02/13/2008 ggiurgiu@fnal.gov typecast pixIter->adc to float
00249         clust_array_2d[irow][icol] = (float)pixIter->adc;
00250     }
00251   
00252   // Make and fill the bool arrays flagging double pixels
00253   std::vector<bool> ydouble(cluster_matrix_size_y), xdouble(cluster_matrix_size_x);
00254   // x directions (shorter), rows
00255   for (int irow = 0; irow < cluster_matrix_size_x; ++irow)
00256     {
00257       xdouble[irow] = theTopol->isItBigPixelInX( irow+row_offset );
00258     }
00259       
00260   // y directions (longer), columns
00261   for (int icol = 0; icol < cluster_matrix_size_y; ++icol) 
00262     {
00263       ydouble[icol] = theTopol->isItBigPixelInY( icol+col_offset );
00264     }
00265 
00266   // Output:
00267   float nonsense = -99999.9; // nonsense init value
00268   templXrec_ = templYrec_ = templSigmaX_ = templSigmaY_ = nonsense;
00269         // If the template recontruction fails, we want to return 1.0 for now
00270         templProbY_ = templProbX_ = templProbQ_ = 1.0;
00271         templQbin_ = 0;
00272         // We have a boolean denoting whether the reco failed or not
00273         hasFilledProb_ = false;
00274         
00275   float templYrec1_ = nonsense;
00276   float templXrec1_ = nonsense;
00277   float templYrec2_ = nonsense;
00278   float templXrec2_ = nonsense;
00279 
00280   // ******************************************************************
00281   // Do it! Use cotalpha_ and cotbeta_ calculated in PixelCPEBase
00282 
00283   GlobalVector bfield = magfield_->inTesla( theDet->surface().position() ); 
00284   
00285   Frame detFrame( theDet->surface().position(), theDet->surface().rotation() );
00286   LocalVector Bfield = detFrame.toLocal( bfield );
00287   float locBz = Bfield.z();
00288     
00289   ierr =
00290     PixelTempReco2D( ID, cotalpha_, cotbeta_,
00291                      locBz, 
00292                      clust_array_2d, ydouble, xdouble,
00293                      templ_,
00294                      templYrec_, templSigmaY_, templProbY_,
00295                      templXrec_, templSigmaX_, templProbX_, 
00296                      templQbin_, 
00297                      speed_,
00298                      templProbQ_
00299                      );
00300 
00301   // ******************************************************************
00302 
00303   // Check exit status
00304   if ( ierr != 0 ) 
00305     {
00306       LogDebug("PixelCPETemplateReco::localPosition") <<
00307         "reconstruction failed with error " << ierr << "\n";
00308 
00309       // Gavril: what do we do in this case ? For now, just return the cluster center of gravity in microns
00310       // In the x case, apply a rough Lorentz drift average correction
00311       // To do: call PixelCPEGeneric whenever PixelTempReco2D fails
00312       double lorentz_drift = -999.9;
00313       if ( thePart == GeomDetEnumerators::PixelBarrel )
00314         lorentz_drift = 60.0; // in microns
00315       else if ( thePart == GeomDetEnumerators::PixelEndcap )
00316         lorentz_drift = 10.0; // in microns
00317       else 
00318         throw cms::Exception("PixelCPETemplateReco::localPosition :") 
00319           << "A non-pixel detector type in here?" << "\n";
00320       // ggiurgiu@jhu.edu, 21/09/2010 : trk angles needed to correct for bows/kinks
00321       if ( with_track_angle )
00322         {
00323           templXrec_ = theTopol->localX( cluster.x(), loc_trk_pred_ ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
00324           templYrec_ = theTopol->localY( cluster.y(), loc_trk_pred_ ); 
00325         }
00326       else
00327         {
00328           edm::LogError("PixelCPETemplateReco") 
00329             << "@SUB = PixelCPETemplateReco::localPosition" 
00330             << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
00331           
00332           templXrec_ = theTopol->localX( cluster.x() ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
00333           templYrec_ = theTopol->localY( cluster.y() );
00334         }
00335     }
00336   else if ( UseClusterSplitter_ && templQbin_ == 0 )
00337     {
00338       cout << " PixelCPETemplateReco : We should never be here !!!!!!!!!!!!!!!!!!!!!!" << endl;
00339       cout << "                 (int)UseClusterSplitter_ = " << (int)UseClusterSplitter_ << endl;
00340 
00341       //ierr = 
00342       //PixelTempSplit( ID, fpix, cotalpha_, cotbeta_, 
00343       //                clust_array_2d, ydouble, xdouble, 
00344       //                templ_, 
00345       //                templYrec1_, templYrec2_, templSigmaY_, templProbY_,
00346       //                templXrec1_, templXrec2_, templSigmaX_, templProbX_, 
00347       //                templQbin_ );
00348 
00349 
00350       float dchisq;
00351       float templProbQ_;
00352       SiPixelTemplate2D templ2D_;
00353       templ2D_.pushfile(ID);
00354         
00355       ierr =
00356         SiPixelTemplateSplit::PixelTempSplit( ID, cotalpha_, cotbeta_,
00357                                               clust_array_2d, 
00358                                               ydouble, xdouble,
00359                                               templ_,
00360                                               templYrec1_, templYrec2_, templSigmaY_, templProbY_,
00361                                               templXrec1_, templXrec2_, templSigmaX_, templProbX_,
00362                                               templQbin_, 
00363                                               templProbQ_,
00364                                               true, 
00365                                               dchisq, 
00366                                               templ2D_ );
00367       
00368 
00369       if ( ierr != 0 )
00370         {
00371           LogDebug("PixelCPETemplateReco::localPosition") <<
00372             "reconstruction failed with error " << ierr << "\n";
00373           
00374           // Gavril: what do we do in this case ? For now, just return the cluster center of gravity in microns
00375           // In the x case, apply a rough Lorentz drift average correction
00376           // To do: call PixelCPEGeneric whenever PixelTempReco2D fails
00377           double lorentz_drift = -999.9;
00378           if ( thePart == GeomDetEnumerators::PixelBarrel )
00379             lorentz_drift = 60.0; // in microns
00380           else if ( thePart == GeomDetEnumerators::PixelEndcap )
00381             lorentz_drift = 10.0; // in microns
00382           else 
00383             throw cms::Exception("PixelCPETemplateReco::localPosition :") 
00384               << "A non-pixel detector type in here?" << "\n";
00385 
00386           // ggiurgiu@jhu.edu, 12/09/2010 : trk angles needed to correct for bows/kinks
00387           if ( with_track_angle )
00388             {
00389               templXrec_ = theTopol->localX( cluster.x(),loc_trk_pred_ ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
00390               templYrec_ = theTopol->localY( cluster.y(),loc_trk_pred_ );
00391             }
00392           else
00393             {
00394               edm::LogError("PixelCPETemplateReco") 
00395                 << "@SUB = PixelCPETemplateReco::localPosition" 
00396                 << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
00397               
00398               templXrec_ = theTopol->localX( cluster.x() ) - lorentz_drift * micronsToCm; // very rough Lorentz drift correction
00399               templYrec_ = theTopol->localY( cluster.y() );
00400               
00401             }
00402         }
00403       else
00404         {
00405           // go from micrometer to centimeter      
00406           templXrec1_ *= micronsToCm;
00407           templYrec1_ *= micronsToCm;     
00408           templXrec2_ *= micronsToCm;
00409           templYrec2_ *= micronsToCm;
00410       
00411           // go back to the module coordinate system   
00412           templXrec1_ += lp.x();
00413           templYrec1_ += lp.y();
00414           templXrec2_ += lp.x();
00415           templYrec2_ += lp.y();
00416       
00417           // calculate distance from each hit to the track and choose the 
00418           // hit closest to the track
00419           float distance11 = sqrt( (templXrec1_ - trk_lp_x)*(templXrec1_ - trk_lp_x) + 
00420                                    (templYrec1_ - trk_lp_y)*(templYrec1_ - trk_lp_y) );
00421           
00422           float distance12 = sqrt( (templXrec1_ - trk_lp_x)*(templXrec1_ - trk_lp_x) + 
00423                                    (templYrec2_ - trk_lp_y)*(templYrec2_ - trk_lp_y) );
00424           
00425           float distance21 = sqrt( (templXrec2_ - trk_lp_x)*(templXrec2_ - trk_lp_x) + 
00426                                    (templYrec1_ - trk_lp_y)*(templYrec1_ - trk_lp_y) );
00427           
00428           float distance22 = sqrt( (templXrec2_ - trk_lp_x)*(templXrec2_ - trk_lp_x) + 
00429                                    (templYrec2_ - trk_lp_y)*(templYrec2_ - trk_lp_y) );
00430           
00431           float min_templXrec_ = -999.9;
00432           float min_templYrec_ = -999.9;
00433           float distance_min = 9999999999.9;
00434           if ( distance11 < distance_min )
00435             {
00436               distance_min = distance11;
00437               min_templXrec_ = templXrec1_;
00438               min_templYrec_ = templYrec1_;
00439             }
00440           if ( distance12 < distance_min )
00441             {
00442               distance_min = distance12;
00443               min_templXrec_ = templXrec1_;
00444               min_templYrec_ = templYrec2_;
00445             }
00446           if ( distance21 < distance_min )
00447             {
00448               distance_min = distance21;
00449               min_templXrec_ = templXrec2_;
00450               min_templYrec_ = templYrec1_;
00451             }
00452           if ( distance22 < distance_min )
00453             {
00454               distance_min = distance22;
00455               min_templXrec_ = templXrec2_;
00456               min_templYrec_ = templYrec2_;
00457             }
00458                   
00459           templXrec_ = min_templXrec_;
00460           templYrec_ = min_templYrec_;
00461         }
00462     } // else if ( UseClusterSplitter_ && templQbin_ == 0 )
00463   else 
00464     {
00465       // go from micrometer to centimeter      
00466       templXrec_ *= micronsToCm;
00467       templYrec_ *= micronsToCm;
00468       
00469       // go back to the module coordinate system 
00470       templXrec_ += lp.x();
00471       templYrec_ += lp.y();
00472     }
00473     
00474   // Save probabilities and qBin in the quantities given to us by the base class
00475   // (for which there are also inline getters).  &&& templProbX_ etc. should be retired...
00476   probabilityX_  = templProbX_;
00477   probabilityY_  = templProbY_;
00478   probabilityQ_  = templProbQ_;
00479   qBin_          = templQbin_;
00480   
00481   if ( ierr == 0 ) 
00482     hasFilledProb_ = true;
00483   
00484   LocalPoint template_lp = LocalPoint( nonsense, nonsense );
00485   template_lp = LocalPoint( templXrec_, templYrec_ );      
00486   
00487   return template_lp;
00488   
00489 }
00490   
00491 //------------------------------------------------------------------
00492 //  localError() relies on localPosition() being called FIRST!!!
00493 //------------------------------------------------------------------
00494 LocalError  
00495 PixelCPETemplateReco::localError( const SiPixelCluster& cluster, 
00496                                   const GeomDetUnit& det ) const 
00497 {
00498   //cout << endl;
00499   //cout << "Set PixelCPETemplate errors .............................................." << endl;
00500 
00501   //cout << "CPETemplate : " << endl;
00502 
00503   //--- Default is the maximum error used for edge clusters.
00504   float xerr = thePitchX / sqrt(12.0);
00505   float yerr = thePitchY / sqrt(12.0);
00506   
00507   // Check if the errors were already set at the clusters splitting level
00508   if ( cluster.getSplitClusterErrorX() > 0.0 && cluster.getSplitClusterErrorX() < 7777.7 && 
00509        cluster.getSplitClusterErrorY() > 0.0 && cluster.getSplitClusterErrorY() < 7777.7 )
00510     {
00511       xerr = cluster.getSplitClusterErrorX() * micronsToCm;
00512       yerr = cluster.getSplitClusterErrorY() * micronsToCm;
00513 
00514       //cout << "Errors set at cluster splitting level : " << endl;
00515       //cout << "xerr = " << xerr << endl;
00516       //cout << "yerr = " << yerr << endl;
00517     }
00518   else
00519     {
00520       // If errors are not split at the cluster splitting level, set the errors here
00521 
00522       //cout  << "Errors are not split at the cluster splitting level, set the errors here : " << endl; 
00523   
00524       setTheDet( det, cluster );
00525       
00526       int maxPixelCol = cluster.maxPixelCol();
00527       int maxPixelRow = cluster.maxPixelRow();
00528       int minPixelCol = cluster.minPixelCol();
00529       int minPixelRow = cluster.minPixelRow();
00530       
00531       //--- Are we near either of the edges?
00532       bool edgex = ( theTopol->isItEdgePixelInX( minPixelRow ) || theTopol->isItEdgePixelInX( maxPixelRow ) );
00533       bool edgey = ( theTopol->isItEdgePixelInY( minPixelCol ) || theTopol->isItEdgePixelInY( maxPixelCol ) );
00534       
00535       if ( ierr !=0 ) 
00536         {
00537           // If reconstruction fails the hit position is calculated from cluster center of gravity 
00538           // corrected in x by average Lorentz drift. Assign huge errors.
00539           //xerr = 10.0 * (float)cluster.sizeX() * xerr;
00540           //yerr = 10.0 * (float)cluster.sizeX() * yerr;
00541           
00542           // Assign better errors based on the residuals for failed template cases
00543           if ( thePart == GeomDetEnumerators::PixelBarrel )
00544             {
00545               xerr = 55.0 * micronsToCm;
00546               yerr = 36.0 * micronsToCm;
00547             }
00548           else if ( thePart == GeomDetEnumerators::PixelEndcap )
00549             {
00550               xerr = 42.0 * micronsToCm;
00551               yerr = 39.0 * micronsToCm;
00552             }
00553           else 
00554             throw cms::Exception("PixelCPETemplateReco::localError :") << "A non-pixel detector type in here?" ;
00555 
00556           //cout << "xerr = " << xerr << endl;
00557           //cout << "yerr = " << yerr << endl;
00558           
00559           //return LocalError(xerr*xerr, 0, yerr*yerr);
00560         }
00561       else if ( edgex || edgey )
00562         {
00563           // for edge pixels assign errors according to observed residual RMS 
00564           if      ( edgex && !edgey )
00565             {
00566               xerr = 23.0 * micronsToCm;
00567               yerr = 39.0 * micronsToCm;
00568             }
00569           else if ( !edgex && edgey )
00570             {
00571               xerr = 24.0 * micronsToCm;
00572               yerr = 96.0 * micronsToCm;
00573             }
00574           else if ( edgex && edgey )
00575             {
00576               xerr = 31.0 * micronsToCm;
00577               yerr = 90.0 * micronsToCm;
00578             }
00579           else
00580             {
00581               throw cms::Exception(" PixelCPETemplateReco::localError: Something wrong with pixel edge flag !!!");
00582             }
00583 
00584           //cout << "xerr = " << xerr << endl;
00585           //cout << "yerr = " << yerr << endl;
00586         }
00587       else 
00588         {
00589           // &&& need a class const
00590           //const float micronsToCm = 1.0e-4;
00591           
00592           xerr = templSigmaX_ * micronsToCm;
00593           yerr = templSigmaY_ * micronsToCm;
00594           
00595           //cout << "xerr = " << xerr << endl;
00596           //cout << "yerr = " << yerr << endl;
00597 
00598           // &&& should also check ierr (saved as class variable) and return
00599           // &&& nonsense (another class static) if the template fit failed.
00600         }       
00601       
00602       if (theVerboseLevel > 9) 
00603         {
00604           LogDebug("PixelCPETemplateReco") <<
00605             " Sizex = " << cluster.sizeX() << " Sizey = " << cluster.sizeY() << " Edgex = " << edgex << " Edgey = " << edgey << 
00606             " ErrX  = " << xerr            << " ErrY  = " << yerr;
00607         }
00608 
00609     } // else
00610   
00611   if ( !(xerr > 0.0) )
00612     throw cms::Exception("PixelCPETemplateReco::localError") 
00613       << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
00614   
00615   if ( !(yerr > 0.0) )
00616     throw cms::Exception("PixelCPETemplateReco::localError") 
00617       << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
00618 
00619   //cout << "Final errors set to: " << endl;
00620   //cout << "xerr = " << xerr << endl;
00621   //cout << "yerr = " << yerr << endl;
00622   //cout << "Out of PixelCPETemplateREco..........................................................................." << endl;
00623   //cout << endl;
00624 
00625   return LocalError(xerr*xerr, 0, yerr*yerr);
00626 }
00627