#include <PixelCPETemplateReco.h>

Inheritance diagram for PixelCPETemplateReco:

Classes
struct	ClusterParamTemplate

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

	~PixelCPETemplateReco ()

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

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

	PixelCPEBase (edm::ParameterSet const &conf, const MagneticField mag, const TrackerGeometry &geom, 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

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

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

Private Attributes
int	speed_

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
typedef GloballyPositioned < double >	Frame

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_

const MagneticField *	magfield_

const SiPixelTemplateDBObject *	templateDBobject_

int	theFlag_

int	theVerboseLevel

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 SiPixelLorentzAngle *	lorentzAngle,
		const SiPixelTemplateDBObject *	templateDBobject
	)

Definition at line 42 of file PixelCPETemplateReco.cc.

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

   : PixelCPEBase(conf, mag, geom, lorentzAngle, 0, templateDBobject, 0,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 40 and 41 from ASCII files ------------------------" << endl;
 
       if ( !SiPixelTemplate::pushfile( 40, thePixelTemp_ ) )
     throw cms::Exception("PixelCPETemplateReco") 
       << "\nERROR: Templates 40 not loaded correctly from text file. Reconstruction will fail.\n\n";
     
       if ( !SiPixelTemplate::pushfile( 41, thePixelTemp_ ) )
     throw cms::Exception("PixelCPETemplateReco") 
       << "\nERROR: Templates 41 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 ( )

Definition at line 96 of file PixelCPETemplateReco.cc.

 {
   // &&& delete template store?
 }

Member Function Documentation

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

privatevirtual

Implements PixelCPEBase.

Definition at line 101 of file PixelCPETemplateReco.cc.

 {
    return new ClusterParamTemplate(cl);
 }

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

privatevirtual

Implements PixelCPEBase.

Definition at line 462 of file PixelCPETemplateReco.cc.

References edm::hlt::Exception, f, PixelCPETemplateReco::ClusterParamTemplate::ierr, RectangularPixelTopology::isItEdgePixelInX(), RectangularPixelTopology::isItEdgePixelInY(), LogDebug, GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, mathSSE::sqrt(), PixelCPETemplateReco::ClusterParamTemplate::templSigmaX_, PixelCPETemplateReco::ClusterParamTemplate::templSigmaY_, PixelCPEBase::DetParam::thePart, PixelCPEBase::DetParam::thePitchX, PixelCPEBase::DetParam::thePitchY, 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.
   const float sig12 = 1./sqrt(12.0);
   float xerr = theDetParam.thePitchX *sig12;
   float yerr = theDetParam.thePitchY *sig12;
   
   // 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;
       
       // Assign better errors based on the residuals for failed template cases
       if ( theDetParam.thePart == GeomDetEnumerators::PixelBarrel )
         {
           xerr = 55.0f * micronsToCm;
           yerr = 36.0f * micronsToCm;
         }
       else if ( theDetParam.thePart == GeomDetEnumerators::PixelEndcap )
         {
           xerr = 42.0f * micronsToCm;
           yerr = 39.0f * micronsToCm;
         }
       else 
         throw cms::Exception("PixelCPETemplateReco::localError :") << "A non-pixel detector type in here?" ;
 
       //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

privatevirtual

Implements PixelCPEBase.

Definition at line 115 of file PixelCPETemplateReco.cc.

References PixelCPEBase::DetParam::bz, gather_cfg::cout, PixelCPEBase::DetParam::detTemplateId, PixelCPEBase::DoLorentz_, edm::hlt::Exception, GeomDet::geographicalId(), SiPixelTemplateDBObject::getTemplateID(), i, PixelCPETemplateReco::ClusterParamTemplate::ierr, RectangularPixelTopology::isItBigPixelInX(), RectangularPixelTopology::isItBigPixelInY(), PixelCPEBase::LoadTemplatesFromDB_, Topology::localPosition(), PixelTopology::localX(), PixelTopology::localY(), LogDebug, PixelCPEBase::DetParam::lorentzShiftInCmX, PixelCPEBase::DetParam::lorentzShiftInCmY, SiPixelTemplate::lorxbias(), SiPixelTemplate::lorybias(), GeomDetEnumerators::PixelBarrel, GeomDetEnumerators::PixelEndcap, SiPixelTemplateReco::PixelTempReco2D(), SiPixelTemplateSplit::PixelTempSplit(), SiPixelTemplate2D::pushfile(), speed_, mathSSE::sqrt(), PixelCPEBase::templateDBobject_, PixelCPETemplateReco::ClusterParamTemplate::templProbQ_, PixelCPETemplateReco::ClusterParamTemplate::templProbX_, PixelCPETemplateReco::ClusterParamTemplate::templProbY_, PixelCPETemplateReco::ClusterParamTemplate::templQbin_, PixelCPETemplateReco::ClusterParamTemplate::templSigmaX_, PixelCPETemplateReco::ClusterParamTemplate::templSigmaY_, PixelCPETemplateReco::ClusterParamTemplate::templXrec_, PixelCPETemplateReco::ClusterParamTemplate::templYrec_, PixelCPEBase::DetParam::theDet, PixelCPEBase::DetParam::thePart, thePixelTemp_, PixelCPEBase::DetParam::theRecTopol, PixelCPEBase::DetParam::theTopol, create_public_lumi_plots::tmp_x, create_public_lumi_plots::tmp_y, unlikely, UseClusterSplitter_, PV3DBase< T, PVType, FrameType >::x(), and PV3DBase< T, PVType, FrameType >::y().

 {
 
   ClusterParamTemplate & theClusterParam = static_cast<ClusterParamTemplate &>(theClusterParamBase);
 
   bool fpix;  //  barrel(false) or forward(true)
   if ( theDetParam.thePart == GeomDetEnumerators::PixelBarrel )   
     fpix = false;    // no, it's not forward -- it's barrel
   else                                              
     fpix = true;     // yes, it's forward
   
   int ID = -9999;
   if ( LoadTemplatesFromDB_ ) {
     int ID0 = templateDBobject_->getTemplateID(theDetParam.theDet->geographicalId()); // just to comapre 
       ID = theDetParam.detTemplateId;
       if(ID0!=ID) cout<<" different id"<< ID<<" "<<ID0<<endl;
   } else { // from asci file 
     if ( !fpix ) ID = 40; // barrel
     else ID = 41; // endcap
   }
   //cout << "PixelCPETemplateReco : ID = " << ID << endl;
   
   SiPixelTemplate templ(thePixelTemp_);
 
   // Make from cluster (a SiPixelCluster) a boost multi_array_2d called 
   // clust_array_2d.
   boost::multi_array<float, 2> clust_array_2d(boost::extents[cluster_matrix_size_x][cluster_matrix_size_y]);
   
   // 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 PixelTempReco2D
   // Will add these values to the output of  PixelTempReco2D
   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) );
     }
     
   // Copy clust's pixels (calibrated in electrons) into clust_array_2d;
   for (int i=0 ; i!=theClusterParam.theCluster->size(); ++i ) 
     {
       auto pix = theClusterParam.theCluster->pixel(i);
       // *pixIter dereferences to Pixel struct, with public vars x, y, adc (all float)
       // 02/13/2008 ggiurgiu@fnal.gov: type of x, y and adc has been changed to unsigned char, unsigned short, unsigned short
       // in DataFormats/SiPixelCluster/interface/SiPixeltheClusterParam.theCluster->h so the type cast to int is redundant. Leave it there, it 
       // won't hurt. 
       int irow = int(pix.x) - row_offset;   // &&& do we need +0.5 ???
       int icol = int(pix.y) - col_offset;   // &&& do we need +0.5 ???
       
       // Gavril : what do we do here if the row/column is larger than cluster_matrix_size_x/cluster_matrix_size_y = 7/21 ?
       // Ignore them for the moment...
       if ( irow<cluster_matrix_size_x && icol<cluster_matrix_size_y )
     // 02/13/2008 ggiurgiu@fnal.gov typecast pixIter->adc to float
     clust_array_2d[irow][icol] = (float)pix.adc;
     }
   
   // Make and fill the bool arrays flagging double pixels
   std::vector<bool> ydouble(cluster_matrix_size_y), xdouble(cluster_matrix_size_x);
   // x directions (shorter), rows
   for (int irow = 0; irow < cluster_matrix_size_x; ++irow)
     {
       xdouble[irow] = theDetParam.theRecTopol->isItBigPixelInX( irow+row_offset );
     }
       
   // y directions (longer), columns
   for (int icol = 0; icol < cluster_matrix_size_y; ++icol) 
     {
       ydouble[icol] = theDetParam.theRecTopol->isItBigPixelInY( icol+col_offset );
     }
 
   // 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;
     
   theClusterParam.ierr =
     PixelTempReco2D( ID, theClusterParam.cotalpha, theClusterParam.cotbeta,
              locBz, 
              clust_array_2d, ydouble, xdouble,
              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 PixelTempReco2D fails
       float lorentz_drift = -999.9;
       if ( theDetParam.thePart == GeomDetEnumerators::PixelBarrel )
     lorentz_drift = 60.0f; // in microns
       else if ( theDetParam.thePart == GeomDetEnumerators::PixelEndcap )
     lorentz_drift = 10.0f; // in microns
       else 
     throw cms::Exception("PixelCPETemplateReco::localPosition :") 
       << "A non-pixel detector type in here?" << "\n";
       // 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_ );
 
 
       float dchisq;
       float templProbQ_;
       std::vector< SiPixelTemplateStore2D > thePixelTemp2D_;
       SiPixelTemplate2D::pushfile(ID, thePixelTemp2D_);
       SiPixelTemplate2D templ2D_(thePixelTemp2D_);
         
       theClusterParam.ierr =
     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 PixelTempReco2D fails
       float lorentz_drift = -999.9f;
       if ( theDetParam.thePart == GeomDetEnumerators::PixelBarrel )
         lorentz_drift = 60.0f; // in microns
       else if ( theDetParam.thePart == GeomDetEnumerators::PixelEndcap )
         lorentz_drift = 10.0f; // in microns
       else 
         throw cms::Exception("PixelCPETemplateReco::localPosition :") 
           << "A non-pixel detector type in here?" << "\n";
 
       // 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 distance11 = sqrt( (templXrec1_ - theClusterParam.trk_lp_x)*(templXrec1_ - theClusterParam.trk_lp_x) + 
                    (templYrec1_ - theClusterParam.trk_lp_y)*(templYrec1_ - theClusterParam.trk_lp_y) );
       
       float distance12 = sqrt( (templXrec1_ - theClusterParam.trk_lp_x)*(templXrec1_ - theClusterParam.trk_lp_x) + 
                    (templYrec2_ - theClusterParam.trk_lp_y)*(templYrec2_ - theClusterParam.trk_lp_y) );
       
       float distance21 = sqrt( (templXrec2_ - theClusterParam.trk_lp_x)*(templXrec2_ - theClusterParam.trk_lp_x) + 
                    (templYrec1_ - theClusterParam.trk_lp_y)*(templYrec1_ - theClusterParam.trk_lp_y) );
       
       float distance22 = sqrt( (templXrec2_ - theClusterParam.trk_lp_x)*(templXrec2_ - theClusterParam.trk_lp_x) + 
                    (templYrec2_ - theClusterParam.trk_lp_y)*(templYrec2_ - theClusterParam.trk_lp_y) );
       
       float min_templXrec_ = -999.9;
       float min_templYrec_ = -999.9;
       float distance_min = 9999999999.9;
       if ( distance11 < distance_min )
         {
           distance_min = distance11;
           min_templXrec_ = templXrec1_;
           min_templYrec_ = templYrec1_;
         }
       if ( distance12 < distance_min )
         {
           distance_min = distance12;
           min_templXrec_ = templXrec1_;
           min_templYrec_ = templYrec2_;
         }
       if ( distance21 < distance_min )
         {
           distance_min = distance21;
           min_templXrec_ = templXrec2_;
           min_templYrec_ = templYrec1_;
         }
       if ( distance22 < distance_min )
         {
           distance_min = distance22;
           min_templXrec_ = templXrec2_;
           min_templYrec_ = templYrec2_;
         }
           
       theClusterParam.templXrec_ = min_templXrec_;
       theClusterParam.templYrec_ = min_templYrec_;
     }
     } // else if ( UseClusterSplitter_ && templQbin_ == 0 )
 
   else // apparenly this is he 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_ );      
   
 }