#include <PixelCPETemplateReco.h>

Inheritance diagram for PixelCPETemplateReco:

Classes
struct	ClusterParamTemplate

Public Member Functions
	PixelCPETemplateReco (edm::ParameterSet const &conf, const MagneticField , const TrackerGeometry &, const TrackerTopology &, const SiPixelLorentzAngle , const SiPixelTemplateDBObject *)

	~PixelCPETemplateReco () override

Public Member Functions inherited from PixelCPEBase
ReturnType	getParameters (const SiPixelCluster &cl, const GeomDetUnit &det) const override

ReturnType	getParameters (const SiPixelCluster &cl, const GeomDetUnit &det, const LocalTrajectoryParameters &ltp) const override

	PixelCPEBase (edm::ParameterSet const &conf, const MagneticField mag, const TrackerGeometry &geom, const TrackerTopology &ttopo, const SiPixelLorentzAngle lorentzAngle, const SiPixelGenErrorDBObject genErrorDBObject, const SiPixelTemplateDBObject templateDBobject, const SiPixelLorentzAngle *lorentzAngleWidth, int flag=0)

Public Member Functions inherited from PixelClusterParameterEstimator
unsigned int	clusterProbComputationFlag () const

virtual ReturnType	getParameters (const SiPixelCluster &cl, const GeomDetUnit &det, const TrajectoryStateOnSurface &tsos) const

virtual VLocalValues	localParametersV (const SiPixelCluster &cluster, const GeomDetUnit &gd) const

virtual VLocalValues	localParametersV (const SiPixelCluster &cluster, const GeomDetUnit &gd, TrajectoryStateOnSurface &tsos) const

	PixelClusterParameterEstimator ()

virtual	~PixelClusterParameterEstimator ()

Private Member Functions
ClusterParam *	createClusterParam (const SiPixelCluster &cl) const override

LocalError	localError (DetParam const &theDetParam, ClusterParam &theClusterParam) const override

LocalPoint	localPosition (DetParam const &theDetParam, ClusterParam &theClusterParam) const override

Private Attributes
int	barrelTemplateID_

int	forwardTemplateID_

int	speed_

std::string	templateDir_

std::vector< SiPixelTemplateStore >	thePixelTemp_

bool	UseClusterSplitter_

Additional Inherited Members
Public Types inherited from PixelClusterParameterEstimator
typedef std::pair< LocalPoint, LocalError >	LocalValues

using	ReturnType = std::tuple< LocalPoint, LocalError, SiPixelRecHitQuality::QualWordType >

typedef std::vector< LocalValues >	VLocalValues

Protected Types inherited from PixelCPEBase
using	DetParams = std::vector< DetParam >

typedef GloballyPositioned< double >	Frame

Protected Member Functions inherited from PixelCPEBase
void	computeAnglesFromDetPosition (DetParam const &theDetParam, ClusterParam &theClusterParam) const

void	computeAnglesFromTrajectory (DetParam const &theDetParam, ClusterParam &theClusterParam, const LocalTrajectoryParameters &ltp) const

void	computeLorentzShifts (DetParam &) const

DetParam const &	detParam (const GeomDetUnit &det) const

LocalVector	driftDirection (DetParam &theDetParam, GlobalVector bfield) const

LocalVector	driftDirection (DetParam &theDetParam, LocalVector bfield) const

void	setTheClu (DetParam const &, ClusterParam &theClusterParam) const

Protected Attributes inherited from PixelCPEBase
bool	alpha2Order

bool	DoLorentz_

const SiPixelGenErrorDBObject *	genErrorDBObject_

const TrackerGeometry &	geom_

float	lAOffset_

float	lAWidthBPix_

float	lAWidthFPix_

bool	LoadTemplatesFromDB_

const SiPixelLorentzAngle *	lorentzAngle_

const SiPixelLorentzAngle *	lorentzAngleWidth_

DetParams	m_DetParams =DetParams(1440)

const MagneticField *	magfield_

const SiPixelTemplateDBObject *	templateDBobject_

int	theFlag_

int	theVerboseLevel

const TrackerTopology &	ttopo_

bool	useLAOffsetFromConfig_

bool	useLAWidthFromConfig_

bool	useLAWidthFromDB_

Protected Attributes inherited from PixelClusterParameterEstimator
unsigned int	clusterProbComputationFlag_

Detailed Description

Definition at line 34 of file PixelCPETemplateReco.h.

Constructor & Destructor Documentation

PixelCPETemplateReco::PixelCPETemplateReco	(	edm::ParameterSet const &	conf,
		const MagneticField *	mag,
		const TrackerGeometry &	geom,
		const TrackerTopology &	ttopo,
		const SiPixelLorentzAngle *	lorentzAngle,
		const SiPixelTemplateDBObject *	templateDBobject
	)

Definition at line 43 of file PixelCPETemplateReco.cc.

References barrelTemplateID_, Exception, forwardTemplateID_, edm::ParameterSet::getParameter(), PixelCPEBase::LoadTemplatesFromDB_, LogDebug, SiPixelTemplate::pushfile(), speed_, PixelCPEBase::templateDBobject_, templateDir_, thePixelTemp_, and UseClusterSplitter_.

 : PixelCPEBase(conf, mag, geom, ttopo, lorentzAngle, nullptr, templateDBobject, nullptr,1)
 {
    //cout << endl;
    //cout << "Constructing PixelCPETemplateReco::PixelCPETemplateReco(...)................................................." << endl;
    //cout << endl;
    
    // Configurable parameters
    //DoCosmics_ = conf.getParameter<bool>("DoCosmics"); // Not used in templates
    //LoadTemplatesFromDB_ = conf.getParameter<bool>("LoadTemplatesFromDB"); // Moved to Base
    
    //cout << " PixelCPETemplateReco : (int)LoadTemplatesFromDB_ = " << (int)LoadTemplatesFromDB_ << endl;
    //cout << "field_magnitude = " << field_magnitude << endl;
    
    // configuration parameter to decide between DB or text file template access
    
    if ( LoadTemplatesFromDB_ )
    {
       //cout << "PixelCPETemplateReco: Loading templates from database (DB) --------- " << endl;
       
       // Initialize template store to the selected ID [Morris, 6/25/08]
       if ( !SiPixelTemplate::pushfile( *templateDBobject_, thePixelTemp_) )
          throw cms::Exception("PixelCPETemplateReco")
          << "\nERROR: Templates not filled correctly. Check the sqlite file. Using SiPixelTemplateDBObject version "
          << (*templateDBobject_).version() << "\n\n";
    }
    else
    {
       //cout << "PixelCPETemplateReco : Loading templates for barrel and forward from ASCII files ----------" << endl;
       barrelTemplateID_  = conf.getParameter<int>( "barrelTemplateID" );
       forwardTemplateID_ = conf.getParameter<int>( "forwardTemplateID" );
       templateDir_       = conf.getParameter<int>( "directoryWithTemplates" );
      
       if ( !SiPixelTemplate::pushfile( barrelTemplateID_  , thePixelTemp_ , templateDir_ ) )
          throw cms::Exception("PixelCPETemplateReco")
      << "\nERROR: Template ID " << barrelTemplateID_ << " not loaded correctly from text file. Reconstruction will fail.\n\n";
       
       if ( !SiPixelTemplate::pushfile( forwardTemplateID_ , thePixelTemp_ , templateDir_ ) )
          throw cms::Exception("PixelCPETemplateReco")
      << "\nERROR: Template ID " << forwardTemplateID_ << " not loaded correctly from text file. Reconstruction will fail.\n\n";
    }
    
    speed_ = conf.getParameter<int>( "speed");
    LogDebug("PixelCPETemplateReco::PixelCPETemplateReco:") <<
    "Template speed = " << speed_ << "\n";
    
    UseClusterSplitter_ = conf.getParameter<bool>("UseClusterSplitter");
    
 }

PixelCPETemplateReco::~PixelCPETemplateReco ( )

override

Definition at line 101 of file PixelCPETemplateReco.cc.

References thePixelTemp_, and x.

 {
    for(auto x : thePixelTemp_) x.destroy();
 }

Member Function Documentation

PixelCPEBase::ClusterParam * PixelCPETemplateReco::createClusterParam ( const SiPixelCluster & cl ) const

overrideprivatevirtual

Implements PixelCPEBase.

Definition at line 106 of file PixelCPETemplateReco.cc.

 {
    return new ClusterParamTemplate(cl);
 }

LocalError PixelCPETemplateReco::localError	(	DetParam const &	theDetParam,
		ClusterParam &	theClusterParam
	)		const

overrideprivatevirtual

Implements PixelCPEBase.

Definition at line 436 of file PixelCPETemplateReco.cc.

References Exception, f, PixelCPETemplateReco::ClusterParamTemplate::ierr, GeomDetEnumerators::isBarrel(), RectangularPixelTopology::isItEdgePixelInX(), RectangularPixelTopology::isItEdgePixelInY(), GeomDetEnumerators::isTrackerPixel(), LogDebug, PixelCPETemplateReco::ClusterParamTemplate::templSigmaX_, PixelCPETemplateReco::ClusterParamTemplate::templSigmaY_, PixelCPEBase::DetParam::thePart, PixelCPEBase::DetParam::theRecTopol, and PixelCPEBase::theVerboseLevel.

 {
    
    ClusterParamTemplate & theClusterParam = static_cast<ClusterParamTemplate &>(theClusterParamBase);
    
    //cout << endl;
    //cout << "Set PixelCPETemplate errors .............................................." << endl;
    
    //cout << "CPETemplate : " << endl;
    
    //--- Default is the maximum error used for edge clusters.
    //--- (never used, in fact: let comment it out, shut up the complains of the static analyzer, and save a few CPU cycles)
    //   const float sig12 = 1./sqrt(12.0);
    //   float xerr = theDetParam.thePitchX *sig12;
    //   float yerr = theDetParam.thePitchY *sig12;
    float xerr, yerr;
    
    // Check if the errors were already set at the clusters splitting level
    if ( theClusterParam.theCluster->getSplitClusterErrorX() > 0.0f && theClusterParam.theCluster->getSplitClusterErrorX() < 7777.7f &&
        theClusterParam.theCluster->getSplitClusterErrorY() > 0.0f && theClusterParam.theCluster->getSplitClusterErrorY() < 7777.7f )
    {
       xerr = theClusterParam.theCluster->getSplitClusterErrorX() * micronsToCm;
       yerr = theClusterParam.theCluster->getSplitClusterErrorY() * micronsToCm;
       
       //cout << "Errors set at cluster splitting level : " << endl;
       //cout << "xerr = " << xerr << endl;
       //cout << "yerr = " << yerr << endl;
    }
    else
    {
       // If errors are not split at the cluster splitting level, set the errors here
       
       //cout  << "Errors are not split at the cluster splitting level, set the errors here : " << endl;
       
       int maxPixelCol = theClusterParam.theCluster->maxPixelCol();
       int maxPixelRow = theClusterParam.theCluster->maxPixelRow();
       int minPixelCol = theClusterParam.theCluster->minPixelCol();
       int minPixelRow = theClusterParam.theCluster->minPixelRow();
       
       //--- Are we near either of the edges?
       bool edgex = ( theDetParam.theRecTopol->isItEdgePixelInX( minPixelRow ) || theDetParam.theRecTopol->isItEdgePixelInX( maxPixelRow ) );
       bool edgey = ( theDetParam.theRecTopol->isItEdgePixelInY( minPixelCol ) || theDetParam.theRecTopol->isItEdgePixelInY( maxPixelCol ) );
       
       if ( theClusterParam.ierr !=0 )
       {
          // If reconstruction fails the hit position is calculated from cluster center of gravity
          // corrected in x by average Lorentz drift. Assign huge errors.
          //xerr = 10.0 * (float)theClusterParam.theCluster->sizeX() * xerr;
          //yerr = 10.0 * (float)theClusterParam.theCluster->sizeX() * yerr;
          
          if(!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
             throw cms::Exception("PixelCPETemplateReco::localPosition :")
             << "A non-pixel detector type in here?";
          
          // Assign better errors based on the residuals for failed template cases
          if ( GeomDetEnumerators::isBarrel(theDetParam.thePart) )
          {
             xerr = 55.0f * micronsToCm;
             yerr = 36.0f * micronsToCm;
          }
          else
          {
             xerr = 42.0f * micronsToCm;
             yerr = 39.0f * micronsToCm;
          }
          
          //cout << "xerr = " << xerr << endl;
          //cout << "yerr = " << yerr << endl;
          
          //return LocalError(xerr*xerr, 0, yerr*yerr);
       }
       else if ( edgex || edgey )
       {
          // for edge pixels assign errors according to observed residual RMS
          if      ( edgex && !edgey )
          {
             xerr = 23.0f * micronsToCm;
             yerr = 39.0f * micronsToCm;
          }
          else if ( !edgex && edgey )
          {
             xerr = 24.0f * micronsToCm;
             yerr = 96.0f * micronsToCm;
          }
          else if ( edgex && edgey )
          {
             xerr = 31.0f * micronsToCm;
             yerr = 90.0f * micronsToCm;
          }
          else
          {
             throw cms::Exception(" PixelCPETemplateReco::localError: Something wrong with pixel edge flag !!!");
          }
          
          //cout << "xerr = " << xerr << endl;
          //cout << "yerr = " << yerr << endl;
       }
       else 
       {
          // &&& need a class const
          //const float micronsToCm = 1.0e-4;
          
          xerr = theClusterParam.templSigmaX_ * micronsToCm;
          yerr = theClusterParam.templSigmaY_ * micronsToCm;
          
          //cout << "xerr = " << xerr << endl;
          //cout << "yerr = " << yerr << endl;
          
          // &&& should also check ierr (saved as class variable) and return
          // &&& nonsense (another class static) if the template fit failed.
       }       
       
       if (theVerboseLevel > 9) 
       {
          LogDebug("PixelCPETemplateReco") <<
          " Sizex = " << theClusterParam.theCluster->sizeX() << " Sizey = " << theClusterParam.theCluster->sizeY() << " Edgex = " << edgex << " Edgey = " << edgey << 
          " ErrX  = " << xerr            << " ErrY  = " << yerr;
       }
       
    } // else
    
    if ( !(xerr > 0.0f) )
       throw cms::Exception("PixelCPETemplateReco::localError") 
       << "\nERROR: Negative pixel error xerr = " << xerr << "\n\n";
    
    if ( !(yerr > 0.0f) )
       throw cms::Exception("PixelCPETemplateReco::localError") 
       << "\nERROR: Negative pixel error yerr = " << yerr << "\n\n";
    
    //cout << "Final errors set to: " << endl;
    //cout << "xerr = " << xerr << endl;
    //cout << "yerr = " << yerr << endl;
    //cout << "Out of PixelCPETemplateREco..........................................................................." << endl;
    //cout << endl;
    
    return LocalError(xerr*xerr, 0, yerr*yerr);
 }

LocalPoint PixelCPETemplateReco::localPosition	(	DetParam const &	theDetParam,
		ClusterParam &	theClusterParam
	)		const

overrideprivatevirtual

std::vector< SiPixelTemplateStore2D > thePixelTemp2D_; SiPixelTemplate2D::pushfile(ID, thePixelTemp2D_); SiPixelTemplate2D templ2D_(thePixelTemp2D_);

Implements PixelCPEBase.

Definition at line 120 of file PixelCPETemplateReco.cc.

 {
    
    ClusterParamTemplate & theClusterParam = static_cast<ClusterParamTemplate &>(theClusterParamBase);
    
    if(!GeomDetEnumerators::isTrackerPixel(theDetParam.thePart))
       throw cms::Exception("PixelCPETemplateReco::localPosition :")
       << "A non-pixel detector type in here?";
    //  barrel(false) or forward(true)
    const bool fpix = GeomDetEnumerators::isEndcap(theDetParam.thePart);
    
    int ID = -9999;
    if ( LoadTemplatesFromDB_ ) {
       int ID0 = templateDBobject_->getTemplateID(theDetParam.theDet->geographicalId()); // just to comapre
       ID = theDetParam.detTemplateId;
       if(ID0!=ID) edm::LogError("PixelCPETemplateReco") <<" different id"<< ID<<" "<<ID0<<endl;
    } else { // from asci file
       if ( !fpix )
     ID = barrelTemplateID_  ; // barrel
       else
     ID = forwardTemplateID_ ; // forward
    }
    //cout << "PixelCPETemplateReco : ID = " << ID << endl;
    
    SiPixelTemplate templ(thePixelTemp_);
    
    // Preparing to retrieve ADC counts from the SiPixeltheClusterParam.theCluster->  In the cluster,
    // we have the following:
    //   int minPixelRow(); // Minimum pixel index in the x direction (low edge).
    //   int maxPixelRow(); // Maximum pixel index in the x direction (top edge).
    //   int minPixelCol(); // Minimum pixel index in the y direction (left edge).
    //   int maxPixelCol(); // Maximum pixel index in the y direction (right edge).
    // So the pixels from minPixelRow() will go into clust_array_2d[0][*],
    // and the pixels from minPixelCol() will go into clust_array_2d[*][0].
    
    int row_offset = theClusterParam.theCluster->minPixelRow();
    int col_offset = theClusterParam.theCluster->minPixelCol();
    
    // Store the coordinates of the center of the (0,0) pixel of the array that
    // gets passed to PixelTempReco1D
    // Will add these values to the output of  PixelTempReco1D
    float tmp_x = float(row_offset) + 0.5f;
    float tmp_y = float(col_offset) + 0.5f;
    
    // Store these offsets (to be added later) in a LocalPoint after tranforming
    // them from measurement units (pixel units) to local coordinates (cm)
    
    // ggiurgiu@jhu.edu 12/09/2010 : update call with trk angles needed for bow/kink corrections
    LocalPoint lp;
    
    if ( theClusterParam.with_track_angle )
       lp = theDetParam.theTopol->localPosition( MeasurementPoint(tmp_x, tmp_y), theClusterParam.loc_trk_pred );
    else
    {
       edm::LogError("PixelCPETemplateReco")
       << "@SUB = PixelCPETemplateReco::localPosition"
       << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
       
       lp = theDetParam.theTopol->localPosition( MeasurementPoint(tmp_x, tmp_y) );
    }
    
    // first compute matrix size
    int mrow=0, mcol=0;
    for (int i=0 ; i!=theClusterParam.theCluster->size(); ++i )
    {
       auto pix = theClusterParam.theCluster->pixel(i);
       int irow = int(pix.x);
       int icol = int(pix.y);
       mrow = std::max(mrow,irow);
       mcol = std::max(mcol,icol);
    }
    mrow -= row_offset; mrow+=1; mrow = std::min(mrow,cluster_matrix_size_x);
    mcol -= col_offset; mcol+=1; mcol = std::min(mcol,cluster_matrix_size_y);
    assert(mrow>0); assert(mcol>0);
    
    float clustMatrix[mrow][mcol];
    memset(clustMatrix,0,sizeof(float)*mrow*mcol);
    
    // Copy clust's pixels (calibrated in electrons) into clusMatrix;
    for (int i=0 ; i!=theClusterParam.theCluster->size(); ++i )
    {
       auto pix = theClusterParam.theCluster->pixel(i);
       int irow = int(pix.x) - row_offset;
       int icol = int(pix.y) - col_offset;
       
       // Gavril : what do we do here if the row/column is larger than cluster_matrix_size_x/cluster_matrix_size_y  ?
       // Ignore them for the moment...
       if ( (irow<mrow) & (icol<mcol) ) clustMatrix[irow][icol] =  float(pix.adc);
       
    }
    
    
    // Make and fill the bool arrays flagging double pixels
    bool xdouble[mrow], ydouble[mcol];
    // x directions (shorter), rows
    for (int irow = 0; irow < mrow; ++irow)
       xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX( irow+row_offset );
    
    // y directions (longer), columns
    for (int icol = 0; icol < mcol; ++icol)
       ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY( icol+col_offset );
    
    SiPixelTemplateReco::ClusMatrix clusterPayload{&clustMatrix[0][0], xdouble, ydouble, mrow,mcol};
    
    // Output:
    float nonsense = -99999.9f; // nonsense init value
    theClusterParam.templXrec_ = theClusterParam.templYrec_ = theClusterParam.templSigmaX_ = theClusterParam.templSigmaY_ = nonsense;
    // If the template recontruction fails, we want to return 1.0 for now
    theClusterParam.templProbY_ = theClusterParam.templProbX_ = theClusterParam.templProbQ_ = 1.0f;
    theClusterParam.templQbin_ = 0;
    // We have a boolean denoting whether the reco failed or not
    theClusterParam.hasFilledProb_ = false;
    
    float templYrec1_ = nonsense;
    float templXrec1_ = nonsense;
    float templYrec2_ = nonsense;
    float templXrec2_ = nonsense;
    
    // ******************************************************************
    // Do it! Use cotalpha_ and cotbeta_ calculated in PixelCPEBase
    
    
    float locBz = theDetParam.bz;
    float locBx = theDetParam.bx;
    
    theClusterParam.ierr =
    PixelTempReco1D( ID, theClusterParam.cotalpha, theClusterParam.cotbeta,
                    locBz, locBx,
                    clusterPayload,
                    templ,
                    theClusterParam.templYrec_, theClusterParam.templSigmaY_, theClusterParam.templProbY_,
                    theClusterParam.templXrec_, theClusterParam.templSigmaX_, theClusterParam.templProbX_,
                    theClusterParam.templQbin_,
                    speed_,
                    theClusterParam.templProbQ_
                    );
    
    // ******************************************************************
    
    // Check exit status
    if UNLIKELY( theClusterParam.ierr != 0 )
    {
       LogDebug("PixelCPETemplateReco::localPosition") <<
       "reconstruction failed with error " << theClusterParam.ierr << "\n";
       
       // Gavril: what do we do in this case ? For now, just return the cluster center of gravity in microns
       // In the x case, apply a rough Lorentz drift average correction
       // To do: call PixelCPEGeneric whenever PixelTempReco1D fails
       float lorentz_drift = -999.9;
       if ( !fpix )
          lorentz_drift = 60.0f; // in microns
       else
          lorentz_drift = 10.0f; // in microns
       // ggiurgiu@jhu.edu, 21/09/2010 : trk angles needed to correct for bows/kinks
       if ( theClusterParam.with_track_angle )
       {
          theClusterParam.templXrec_ = theDetParam.theTopol->localX( theClusterParam.theCluster->x(), theClusterParam.loc_trk_pred ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
          theClusterParam.templYrec_ = theDetParam.theTopol->localY( theClusterParam.theCluster->y(), theClusterParam.loc_trk_pred );
       }
       else
       {
          edm::LogError("PixelCPETemplateReco")
          << "@SUB = PixelCPETemplateReco::localPosition"
          << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
          
          theClusterParam.templXrec_ = theDetParam.theTopol->localX( theClusterParam.theCluster->x() ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
          theClusterParam.templYrec_ = theDetParam.theTopol->localY( theClusterParam.theCluster->y() );
       }
    }
    else if UNLIKELY( UseClusterSplitter_ && theClusterParam.templQbin_ == 0 )
    {
       cout << " PixelCPETemplateReco : We should never be here !!!!!!!!!!!!!!!!!!!!!!" << endl;
       cout << "                 (int)UseClusterSplitter_ = " << (int)UseClusterSplitter_ << endl;
       
       //ierr =
       //PixelTempSplit( ID, fpix, cotalpha_, cotbeta_,
       //        clust_array_2d, ydouble, xdouble,
       //        templ,
       //        templYrec1_, templYrec2_, templSigmaY_, templProbY_,
       //        templXrec1_, templXrec2_, templSigmaX_, templProbX_,
       //        templQbin_ );
       
 
       // Commented for now (3/10/17) until we figure out how to resuscitate 2D template splitter
       
       theClusterParam.ierr = -123;
       /*
        float dchisq;
        float templProbQ_;
        SiPixelTemplateSplit::PixelTempSplit( ID, theClusterParam.cotalpha, theClusterParam.cotbeta,
        clust_array_2d,
        ydouble, xdouble,
        templ,
        templYrec1_, templYrec2_, theClusterParam.templSigmaY_, theClusterParam.templProbY_,
        templXrec1_, templXrec2_, theClusterParam.templSigmaX_, theClusterParam.templProbX_,
        theClusterParam.templQbin_,
        templProbQ_,
        true,
        dchisq,
        templ2D_ );
        
        */
       if ( theClusterParam.ierr != 0 )
       {
          LogDebug("PixelCPETemplateReco::localPosition") <<
          "reconstruction failed with error " << theClusterParam.ierr << "\n";
          
          // Gavril: what do we do in this case ? For now, just return the cluster center of gravity in microns
          // In the x case, apply a rough Lorentz drift average correction
          // To do: call PixelCPEGeneric whenever PixelTempReco1D fails
          float lorentz_drift = -999.9f;
          if ( !fpix )
             lorentz_drift = 60.0f; // in microns
          else
             lorentz_drift = 10.0f; // in microns
          
          // ggiurgiu@jhu.edu, 12/09/2010 : trk angles needed to correct for bows/kinks
          if ( theClusterParam.with_track_angle )
          {
             theClusterParam.templXrec_ = theDetParam.theTopol->localX( theClusterParam.theCluster->x(),theClusterParam.loc_trk_pred ) - lorentz_drift * micronsToCm; // rough Lorentz drift correction
             theClusterParam.templYrec_ = theDetParam.theTopol->localY( theClusterParam.theCluster->y(),theClusterParam.loc_trk_pred );
          }
          else
          {
             edm::LogError("PixelCPETemplateReco")
             << "@SUB = PixelCPETemplateReco::localPosition"
             << "Should never be here. PixelCPETemplateReco should always be called with track angles. This is a bad error !!! ";
             
             theClusterParam.templXrec_ = theDetParam.theTopol->localX( theClusterParam.theCluster->x() ) - lorentz_drift * micronsToCm; // very rough Lorentz drift correction
             theClusterParam.templYrec_ = theDetParam.theTopol->localY( theClusterParam.theCluster->y() );
             
          }
       }
       else
       {
          // go from micrometer to centimeter
          templXrec1_ *= micronsToCm;
          templYrec1_ *= micronsToCm;
          templXrec2_ *= micronsToCm;
          templYrec2_ *= micronsToCm;
          
          // go back to the module coordinate system
          templXrec1_ += lp.x();
          templYrec1_ += lp.y();
          templXrec2_ += lp.x();
          templYrec2_ += lp.y();
          
          // calculate distance from each hit to the track and choose the hit closest to the track
          float distX1 = std::abs (templXrec1_ - theClusterParam.trk_lp_x);
          float distX2 = std::abs (templXrec2_ - theClusterParam.trk_lp_x);
          float distY1 = std::abs (templYrec1_ - theClusterParam.trk_lp_y);
          float distY2 = std::abs (templYrec2_ - theClusterParam.trk_lp_y);
          theClusterParam.templXrec_ = (distX1<distX2? templXrec1_ : templXrec2_);
          theClusterParam.templYrec_ = (distY1<distY2? templYrec1_ : templYrec2_);
       }
    } // else if ( UseClusterSplitter_ && templQbin_ == 0 )
    
    else // apparenly this is the good one!
    {
       // go from micrometer to centimeter
       theClusterParam.templXrec_ *= micronsToCm;
       theClusterParam.templYrec_ *= micronsToCm;
       
       // go back to the module coordinate system
       theClusterParam.templXrec_ += lp.x();
       theClusterParam.templYrec_ += lp.y();
       
       // Compute the Alignment Group Corrections [template ID should already be selected from call to reco procedure]
       if ( DoLorentz_ ) {
          // Do only if the lotentzshift has meaningfull numbers
          if( theDetParam.lorentzShiftInCmX!= 0.0 ||  theDetParam.lorentzShiftInCmY!= 0.0 ) {
             // the LA width/shift returned by templates use (+)
             // the LA width/shift produced by PixelCPEBase for positive LA is (-)
             // correct this by iserting (-)
             //float temp1 = -micronsToCm*templ.lorxwidth();  // old
             //float temp2 = -micronsToCm*templ.lorywidth();  // does not incl 1/2
             float templateLorbiasCmX = -micronsToCm*templ.lorxbias();  // new
             float templateLorbiasCmY = -micronsToCm*templ.lorybias(); //incl. 1/2
             // now, correctly, we can use the difference of shifts
             //theClusterParam.templXrec_ += 0.5*(theDetParam.lorentzShiftInCmX - templateLorbiasCmX);
             //theClusterParam.templYrec_ += 0.5*(theDetParam.lorentzShiftInCmY - templateLorbiasCmY);
             theClusterParam.templXrec_ += (0.5*(theDetParam.lorentzShiftInCmX) - templateLorbiasCmX);
             theClusterParam.templYrec_ += (0.5*(theDetParam.lorentzShiftInCmY) - templateLorbiasCmY);
             //cout << "Templates: la lorentz offset = "
             //   <<(0.5*(theDetParam.lorentzShiftInCmX)-templateLorbiasCmX)
             //   <<" "<<templateLorbiasCmX<<" "<<templateLorbiasCmY
             //   <<" "<<temp1<<" "<<temp2
             //   <<" "<<theDetParam.lorentzShiftInCmX
             //   <<" "<<theDetParam.lorentzShiftInCmY
             //   << endl; //dk
          } //else {cout<<" LA is 0, disable offset corrections "<<endl;} //dk
       } //else {cout<<" Do not do LA offset correction "<<endl;} //dk
       
    }
    
    // Save probabilities and qBin in the quantities given to us by the base class
    // (for which there are also inline getters).  &&& templProbX_ etc. should be retired...
    theClusterParam.probabilityX_  = theClusterParam.templProbX_;
    theClusterParam.probabilityY_  = theClusterParam.templProbY_;
    theClusterParam.probabilityQ_  = theClusterParam.templProbQ_;
    theClusterParam.qBin_          = theClusterParam.templQbin_;
    
    if ( theClusterParam.ierr == 0 ) // always true here
       theClusterParam.hasFilledProb_ = true;
    
    return LocalPoint( theClusterParam.templXrec_, theClusterParam.templYrec_ );
    
 }