CMS 3D CMS Logo

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

Go to the documentation of this file.
00001 // Move geomCorrection to the concrete class. d.k. 06/06.
00002 // Change drift direction. d.k. 06/06
00003 
00004 // G. Giurgiu (ggiurgiu@pha.jhu.edu), 12/01/06, implemented the function: 
00005 // computeAnglesFromDetPosition(const SiPixelCluster & cl, 
00006 //                              const GeomDetUnit    & det ) const
00007 //                                    09/09/07, replaced assert statements with throw cms::Exception 
00008 //                                              and fix an invalid pointer check in setTheDet function 
00009 //                                    09/21/07, implement caching of Lorentz drift direction
00010 //                                    01/24/09, use atan2 to get the alpha and beta angles
00011 // change to use Lorentz angle from DB Lotte Wilke, Jan. 31st, 2008
00012 // Change to use Generic error & Template calibration from DB - D.Fehling 11/08
00013 
00014 
00015 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00016 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
00017 
00018 #include "RecoLocalTracker/SiPixelRecHits/interface/PixelCPEBase.h"
00019 
00020 //#define TPDEBUG
00021 #define CORRECT_FOR_BIG_PIXELS
00022 
00023 // MessageLogger
00024 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00025 
00026 // Magnetic field
00027 #include "MagneticField/Engine/interface/MagneticField.h"
00028 
00029 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00030 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00031 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00032 
00033 #include <iostream>
00034 
00035 using namespace std;
00036 
00037 const float PI = 3.141593;
00038 const float degsPerRad = 57.29578;
00039 
00040 //-----------------------------------------------------------------------------
00041 //  A fairly boring constructor.  All quantities are DetUnit-dependent, and
00042 //  will be initialized in setTheDet().
00043 //-----------------------------------------------------------------------------
00044 PixelCPEBase::PixelCPEBase(edm::ParameterSet const & conf, const MagneticField *mag, const SiPixelLorentzAngle * lorentzAngle, 
00045                            const SiPixelCPEGenericErrorParm * genErrorParm, const SiPixelTemplateDBObject * templateDBobject)
00046   : theDet(0), nRecHitsTotal_(0), nRecHitsUsedEdge_(0),
00047     probabilityX_(0.0), probabilityY_(0.0),
00048     probabilityQ_(0.0), qBin_(0),
00049     isOnEdge_(0), hasBadPixels_(0),
00050     spansTwoROCs_(0),
00051     hasFilledProb_(0),
00052     loc_trk_pred_(0.0, 0.0, 0.0, 0.0)
00053 {
00054   //--- Lorentz angle tangent per Tesla
00055   //   theTanLorentzAnglePerTesla =
00056 
00057   //     conf.getParameter<double>("TanLorentzAnglePerTesla");
00058   lorentzAngle_ = lorentzAngle;
00059   /*    if(!lorentzAngle_)
00060     theTanLorentzAnglePerTesla =
00061     conf.getParameter<double>("TanLorentzAnglePerTesla");
00062   */
00063   
00064   //--- Algorithm's verbosity
00065   theVerboseLevel = 
00066     conf.getUntrackedParameter<int>("VerboseLevel",0);
00067   
00068   //-- Magnetic Field
00069   magfield_ = mag;
00070   
00071   //-- Error Parametriaztion from DB for CPE Generic
00072   genErrorParm_ = genErrorParm;
00073   
00074   //-- Template Calibration Object from DB
00075   templateDBobject_ = templateDBobject;
00076   
00077   //-- Switch on/off E.B 
00078   alpha2Order = conf.getParameter<bool>("Alpha2Order");
00079   
00080   //--- A flag that could be used to change the behavior of
00081   //--- clusterProbability() in TSiPixelRecHit (the *transient* one).  
00082   //--- The problem is that the transient hits are made after the CPE runs
00083   //--- and they don't get the access to the PSet, so we pass it via the
00084   //--- CPE itself...
00085   //
00086   clusterProbComputationFlag_ 
00087     = (unsigned int) conf.getParameter<int>("ClusterProbComputationFlag");
00088   
00089 }
00090 
00091 //-----------------------------------------------------------------------------
00092 //  One function to cache the variables common for one DetUnit.
00093 //-----------------------------------------------------------------------------
00094 void
00095 PixelCPEBase::setTheDet( const GeomDetUnit & det, const SiPixelCluster & cluster ) const 
00096 {
00097   if ( theDet == &det )
00098     return;       // we have already seen this det unit
00099   
00100   //--- This is a new det unit, so cache it
00101   theDet = dynamic_cast<const PixelGeomDetUnit*>( &det );
00102 
00103   if ( !theDet ) 
00104     {
00105       throw cms::Exception(" PixelCPEBase::setTheDet : ")
00106             << " Wrong pointer to PixelGeomDetUnit object !!!";
00107     }
00108 
00109   //--- theDet->type() returns a GeomDetType, which implements subDetector()
00110   thePart = theDet->type().subDetector();
00111   switch ( thePart ) 
00112     {
00113     case GeomDetEnumerators::PixelBarrel:
00114       // A barrel!  A barrel!
00115       break;
00116     case GeomDetEnumerators::PixelEndcap:
00117       // A forward!  A forward!
00118       break;
00119     default:
00120       throw cms::Exception("PixelCPEBase::setTheDet :")
00121         << "PixelCPEBase: A non-pixel detector type in here?" ;
00122     }
00123 
00124   //--- The location in of this DetUnit in a cyllindrical coord system (R,Z)
00125   //--- The call goes via BoundSurface, returned by theDet->surface(), but
00126   //--- position() is implemented in GloballyPositioned<> template
00127   //--- ( BoundSurface : Surface : GloballyPositioned<float> )
00128   theDetR = theDet->surface().position().perp();
00129   theDetZ = theDet->surface().position().z();
00130   //--- Define parameters for chargewidth calculation
00131 
00132   //--- bounds() is implemented in BoundSurface itself.
00133   theThickness = theDet->surface().bounds().thickness();
00134 
00135   //--- Cache the topology.
00136   // ggiurgiu@jhu.edu 12/09/2010 : no longer need to dynamyc cast to RectangularPixelTopology
00137   //theTopol
00138   //= dynamic_cast<const RectangularPixelTopology*>( & (theDet->specificTopology()) );
00139 
00140   theTopol = &(theDet->specificTopology());
00141 
00142   //---- The geometrical description of one module/plaquette
00143   theNumOfRow = theTopol->nrows();      // rows in x
00144   theNumOfCol = theTopol->ncolumns();   // cols in y
00145   std::pair<float,float> pitchxy = theTopol->pitch();
00146   thePitchX = pitchxy.first;            // pitch along x
00147   thePitchY = pitchxy.second;           // pitch along y
00148 
00149   //--- Find the offset
00150   // ggiurgiu@jhu.edu 12/09/2010 : this piece is deprecated
00151   //MeasurementPoint offset = 
00152   //theTopol->measurementPosition( LocalPoint(0.0, 0.0), 
00153   //                       Topology::LocalTrackAngles( LocalTrajectoryParameters::dxdz(), 
00154   //                                                   LocalTrajectoryParameters::dydz() ) );  
00155   //
00156   //theOffsetX = offset.x();
00157   //theOffsetY = offset.y();
00158 
00159   //--- Find if the E field is flipped: i.e. whether it points
00160   //--- from the beam, or towards the beam.  (The voltages are
00161   //--- applied backwards on every other module in barrel and
00162   //--- blade in forward.)
00163   theSign = isFlipped() ? -1 : 1;
00164 
00165   //--- The Lorentz shift.
00166   theLShiftX = lorentzShiftX();
00167 
00168   theLShiftY = lorentzShiftY();
00169 
00170   // testing 
00171   if(thePart == GeomDetEnumerators::PixelBarrel) {
00172     //cout<<" lorentz shift "<<theLShiftX<<" "<<theLShiftY<<endl;
00173     theLShiftY=0.;
00174   }
00175 
00176         //--- Geometric Quality Information
00177         int minInX,minInY,maxInX,maxInY=0;
00178         minInX = cluster.minPixelRow();
00179         minInY = cluster.minPixelCol();
00180         maxInX = cluster.maxPixelRow();
00181         maxInY = cluster.maxPixelCol();
00182 
00183         if(theTopol->isItEdgePixelInX(minInX) || theTopol->isItEdgePixelInX(maxInX) ||
00184                  theTopol->isItEdgePixelInY(minInY) || theTopol->isItEdgePixelInY(maxInY) )  {
00185                 isOnEdge_ = true;
00186         }
00187         else isOnEdge_ = false;
00188 
00189         // Bad Pixels have their charge set to 0 in the clusterizer
00190         hasBadPixels_ = false;
00191         for(unsigned int i=0; i<cluster.pixelADC().size(); ++i) {
00192                 if(cluster.pixelADC()[i] == 0) hasBadPixels_ = true;
00193         }
00194 
00195         if(theTopol->containsBigPixelInX(minInX,maxInX) ||
00196                  theTopol->containsBigPixelInY(minInY,maxInY) )  {
00197                 spansTwoROCs_ = true;
00198         }
00199         else spansTwoROCs_ = false;
00200         
00201 
00202         if (theVerboseLevel > 1) 
00203     {
00204       LogDebug("PixelCPEBase") << "***** PIXEL LAYOUT *****" 
00205                                << " thePart = " << thePart
00206                                << " theThickness = " << theThickness
00207                                << " thePitchX  = " << thePitchX 
00208                                << " thePitchY  = " << thePitchY 
00209         // << " theOffsetX = " << theOffsetX 
00210         // << " theOffsetY = " << theOffsetY 
00211                                << " theLShiftX  = " << theLShiftX;
00212     }
00213         
00214 }
00215 
00216 //-----------------------------------------------------------------------------
00217 //  Compute alpha_ and beta_ from the position of this DetUnit.
00218 //  &&& DOES NOT WORK FOR NOW. d.k. 6/06
00219 // The angles from dets are calculated internaly in the PixelCPEInitial class
00220 //-----------------------------------------------------------------------------
00221 // G. Giurgiu, 12/01/06 : implement the function
00222 void PixelCPEBase::
00223 computeAnglesFromDetPosition(const SiPixelCluster & cl, 
00224                              const GeomDetUnit    & det ) const
00225 {
00226   //--- This is a new det unit, so cache it
00227   theDet = dynamic_cast<const PixelGeomDetUnit*>( &det );
00228   if ( ! theDet ) 
00229     {
00230       throw cms::Exception("PixelCPEBase::computeAngleFromDetPosition")
00231         << " Wrong pointer to pixel detector !!!" << endl;
00232     
00233     }
00234 
00235   // get cluster center of gravity (of charge)
00236   float xcenter = cl.x();
00237   float ycenter = cl.y();
00238   
00239   // get the cluster position in local coordinates (cm)
00240 
00241   // ggiurgiu@jhu.edu 12/09/2010 : This function is called without track info, therefore there are no track 
00242   // angles to provide here. Call the default localPosition (without track info)
00243   LocalPoint lp = theTopol->localPosition( MeasurementPoint(xcenter, ycenter) );
00244 
00245   // get the cluster position in global coordinates (cm)
00246   GlobalPoint gp = theDet->surface().toGlobal( lp );
00247   float gp_mod = sqrt( gp.x()*gp.x() + gp.y()*gp.y() + gp.z()*gp.z() );
00248 
00249   // normalize
00250   float gpx = gp.x()/gp_mod;
00251   float gpy = gp.y()/gp_mod;
00252   float gpz = gp.z()/gp_mod;
00253 
00254   // make a global vector out of the global point; this vector will point from the 
00255   // origin of the detector to the cluster
00256   GlobalVector gv(gpx, gpy, gpz);
00257 
00258   // make local unit vector along local X axis
00259   const Local3DVector lvx(1.0, 0.0, 0.0);
00260 
00261   // get the unit X vector in global coordinates/
00262   GlobalVector gvx = theDet->surface().toGlobal( lvx );
00263 
00264   // make local unit vector along local Y axis
00265   const Local3DVector lvy(0.0, 1.0, 0.0);
00266 
00267   // get the unit Y vector in global coordinates
00268   GlobalVector gvy = theDet->surface().toGlobal( lvy );
00269    
00270   // make local unit vector along local Z axis
00271   const Local3DVector lvz(0.0, 0.0, 1.0);
00272 
00273   // get the unit Z vector in global coordinates
00274   GlobalVector gvz = theDet->surface().toGlobal( lvz );
00275     
00276   // calculate the components of gv (the unit vector pointing to the cluster) 
00277   // in the local coordinate system given by the basis {gvx, gvy, gvz}
00278   // note that both gv and the basis {gvx, gvy, gvz} are given in global coordinates
00279   float gv_dot_gvx = gv.x()*gvx.x() + gv.y()*gvx.y() + gv.z()*gvx.z();
00280   float gv_dot_gvy = gv.x()*gvy.x() + gv.y()*gvy.y() + gv.z()*gvy.z();
00281   float gv_dot_gvz = gv.x()*gvz.x() + gv.y()*gvz.y() + gv.z()*gvz.z();
00282 
00283   // calculate angles
00284   alpha_ = atan2( gv_dot_gvz, gv_dot_gvx );
00285   beta_  = atan2( gv_dot_gvz, gv_dot_gvy );
00286 
00287   cotalpha_ = gv_dot_gvx / gv_dot_gvz;
00288   cotbeta_  = gv_dot_gvy / gv_dot_gvz;
00289 
00290   with_track_angle = false;
00291 }
00292 
00293 //-----------------------------------------------------------------------------
00294 //  Compute alpha_ and beta_ from the LocalTrajectoryParameters.
00295 //  Note: should become const after both localParameters() become const.
00296 //-----------------------------------------------------------------------------
00297 void PixelCPEBase::
00298 computeAnglesFromTrajectory( const SiPixelCluster & cl,
00299                              const GeomDetUnit    & det, 
00300                              const LocalTrajectoryParameters & ltp) const
00301 {
00302   loc_traj_param_ = ltp;
00303 
00304   LocalVector localDir = ltp.momentum()/ltp.momentum().mag();
00305   
00306   // &&& Or, maybe we need to move to the local frame ???
00307   //  LocalVector localDir( theDet->toLocal(theState.globalDirection()));
00308   //thePart = theDet->type().part();
00309   
00310   float locx = localDir.x();
00311   float locy = localDir.y();
00312   float locz = localDir.z();
00313 
00314   /*
00315     // Danek's definition 
00316     alpha_ = acos(locx/sqrt(locx*locx+locz*locz));
00317     if ( isFlipped() )                    // &&& check for FPIX !!!
00318     alpha_ = PI - alpha_ ;
00319     beta_ = acos(locy/sqrt(locy*locy+locz*locz));
00320   */
00321 
00322   // &&& In the above, why not use atan2() ?
00323   // ggiurgiu@fnal.gov, 01/24/09 :  Use it now.
00324   alpha_ = atan2( locz, locx );
00325   beta_  = atan2( locz, locy );
00326 
00327   cotalpha_ = locx/locz;
00328   cotbeta_  = locy/locz;
00329 
00330   LocalPoint trk_lp = ltp.position();
00331   trk_lp_x = trk_lp.x();
00332   trk_lp_y = trk_lp.y();
00333   
00334   with_track_angle = true;
00335 
00336 
00337   // ggiurgiu@jhu.edu 12/09/2010 : needed to correct for bows/kinks
00338   AlgebraicVector5 vec_trk_parameters = ltp.mixedFormatVector();
00339   //loc_trk_pred = &Topology::LocalTrackPred( vec_trk_parameters );
00340   loc_trk_pred_ = Topology::LocalTrackPred( vec_trk_parameters );
00341   
00342 }
00343 
00344 //-----------------------------------------------------------------------------
00345 //  Estimate theAlpha for barrel, based on the det position.
00346 //  &&& Needs to be consolidated from the above.
00347 //-----------------------------------------------------------------------------
00348 //float 
00349 //PixelCPEBase::estimatedAlphaForBarrel(float centerx) const
00350 //{
00351 //  float tanalpha = theSign * (centerx-theOffsetX) * thePitchX / theDetR;
00352 //  return PI/2.0 - atan(tanalpha);
00353 //}
00354 
00355 //-----------------------------------------------------------------------------
00356 //  The local position.
00357 // Should do correctly the big pixels.
00358 //-----------------------------------------------------------------------------
00359 LocalPoint
00360 PixelCPEBase::localPosition( const SiPixelCluster& cluster, 
00361                              const GeomDetUnit & det) const {
00362   setTheDet( det, cluster );
00363   
00364   float lpx = xpos(cluster);
00365   float lpy = ypos(cluster);
00366   float lxshift = theLShiftX * thePitchX;  // shift in cm
00367   float lyshift = theLShiftY * thePitchY;
00368   LocalPoint cdfsfs(lpx-lxshift, lpy-lyshift);
00369   return cdfsfs;
00370 }
00371 
00372 //-----------------------------------------------------------------------------
00373 //  Seems never used?
00374 //-----------------------------------------------------------------------------
00375 MeasurementPoint 
00376 PixelCPEBase::measurementPosition( const SiPixelCluster& cluster, 
00377                                    const GeomDetUnit & det) const {
00378 
00379   LocalPoint lp = localPosition(cluster,det);
00380 
00381   // ggiurgiu@jhu.edu 12/09/2010 : trk angles needed for bow/kink correction
00382 
00383   if ( with_track_angle )
00384     return theTopol->measurementPosition( lp, Topology::LocalTrackAngles( loc_traj_param_.dxdz(), loc_traj_param_.dydz() ) );
00385   else 
00386     return theTopol->measurementPosition( lp );
00387 
00388 }
00389 
00390 
00391 //-----------------------------------------------------------------------------
00392 //  Once we have the position, feed it to the topology to give us
00393 //  the error.  
00394 //  &&& APPARENTLY THIS METHOD IS NOT BEING USED ??? (Delete it?)
00395 //-----------------------------------------------------------------------------
00396 MeasurementError  
00397 PixelCPEBase::measurementError( const SiPixelCluster& cluster, const GeomDetUnit & det) const 
00398 {
00399   LocalPoint lp( localPosition(cluster, det) );
00400   LocalError le( localError(   cluster, det) );
00401 
00402   // ggiurgiu@jhu.edu 12/09/2010 : trk angles needed for bow/kink correction
00403   if ( with_track_angle )
00404     return theTopol->measurementError( lp, le, Topology::LocalTrackAngles( loc_traj_param_.dxdz(), loc_traj_param_.dydz() ) );
00405   else 
00406     return theTopol->measurementError( lp, le );
00407 }
00408 
00409 //-----------------------------------------------------------------------------
00410 // The isFlipped() is a silly way to determine which detectors are inverted.
00411 // In the barrel for every 2nd ladder the E field direction is in the
00412 // global r direction (points outside from the z axis), every other
00413 // ladder has the E field inside. Something similar is in the 
00414 // forward disks (2 sides of the blade). This has to be recognised
00415 // because the charge sharing effect is different.
00416 //
00417 // The isFliped does it by looking and the relation of the local (z always
00418 // in the E direction) to global coordinates. There is probably a much 
00419 // better way.
00420 //-----------------------------------------------------------------------------
00421 bool PixelCPEBase::isFlipped() const 
00422 {
00423   // Check the relative position of the local +/- z in global coordinates.
00424   float tmp1 = theDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
00425   float tmp2 = theDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
00426   //cout << " 1: " << tmp1 << " 2: " << tmp2 << endl;
00427   if ( tmp2<tmp1 ) return true;
00428   else return false;    
00429 }
00430 
00431 //-----------------------------------------------------------------------------
00432 // HALF OF the Lorentz shift (so for the full shift multiply by 2), and
00433 // in the units of pitch.  (So note these are neither local nor measurement
00434 // units!)
00435 //-----------------------------------------------------------------------------
00436 float PixelCPEBase::lorentzShiftX() const 
00437 {
00438   LocalVector dir;
00439   Param & p = const_cast<PixelCPEBase*>(this)->m_Params[ theDet->geographicalId().rawId() ];
00440   if ( p.topology ) 
00441     {
00442       //cout << "--------------- old ----------------------" << endl;
00443       //cout << "p.topology = " << p.topology << endl;
00444       dir = p.drift;
00445       //cout << "same direction: dir = " << dir << endl;
00446     }
00447   else 
00448     {
00449       //cout << "--------------- new ----------------------" << endl;
00450       //cout << "p.topology = " << p.topology << endl;
00451       
00452       // ggiurgiu@jhu.edu 12/09/2010 : no longer need to cast
00453       //p.topology = (RectangularPixelTopology*)( & ( theDet->specificTopology() ) );          
00454       p.topology = &( theDet->specificTopology() );   
00455 
00456       p.drift = driftDirection(magfield_->inTesla(theDet->surface().position()) );
00457       dir = p.drift;
00458       //cout << "p.topology = " << p.topology << endl;
00459       //cout << "new direction: dir = " << dir << endl;
00460 
00461     }
00462   //LocalVector dir = driftDirection(magfield_->inTesla(theDet->surface().position()) );
00463   
00464   // max shift in cm 
00465   float xdrift = dir.x()/dir.z() * theThickness;  
00466   // express the shift in units of pitch, 
00467   // divide by 2 to get the average correction
00468   float lshift = xdrift / thePitchX / 2.; 
00469   
00470   //cout << "Lorentz Drift = " << lshift << endl;
00471   //cout << "X Drift = " << dir.x() << endl;
00472   //cout << "Z Drift = " << dir.z() << endl;
00473   
00474   return lshift;  
00475   
00476 
00477 }
00478 
00479 float PixelCPEBase::lorentzShiftY() const 
00480 {
00481   
00482   LocalVector dir;
00483   
00484   Param & p = const_cast<PixelCPEBase*>(this)->m_Params[ theDet->geographicalId().rawId() ];
00485   if ( p.topology ) 
00486     {
00487       //cout << "--------------- old y ----------------------" << endl;
00488       //cout << "p.topology y = " << p.topology << endl;
00489       dir = p.drift;
00490       //cout << "same direction y: dir = " << dir << endl;
00491     }
00492   else 
00493     {
00494       //cout << "--------------- new y ----------------------" << endl;
00495       //cout << "p.topology y = " << p.topology << endl;
00496       
00497       //p.topology = (RectangularPixelTopology*)( & ( theDet->specificTopology() ) );    
00498       p.topology = &( theDet->specificTopology() );
00499 
00500       p.drift = driftDirection(magfield_->inTesla(theDet->surface().position()) );
00501       dir = p.drift;
00502       //cout << "p.topology y = " << p.topology << endl;
00503       //cout << "new direction y: dir = " << dir << endl;
00504 
00505     }
00506 
00507   //LocalVector dir = driftDirection(magfield_->inTesla(theDet->surface().position()) );
00508   
00509   float ydrift = dir.y()/dir.z() * theThickness;
00510   float lshift = ydrift / thePitchY / 2.;
00511   return lshift; 
00512   
00513 
00514 }
00515 
00516 //-----------------------------------------------------------------------------
00517 //  Sum the pixels in the first and the last row, and the total.  Returns
00518 //  a vector of three elements with q_first, q_last and q_total.
00519 //  &&& Really need a simpler & cleaner way, this is very confusing...
00520 //-----------------------------------------------------------------------------
00521 void PixelCPEBase::xCharge(const vector<SiPixelCluster::Pixel>& pixelsVec, 
00522                            const int& imin, const int& imax,
00523                            float& q1, float& q2) const {
00524   //calculate charge in the first and last pixel in y
00525   // and the total cluster charge
00526   q1 = 0.0; 
00527   q2 = 0.0;
00528   float qm = 0.0;
00529   int isize = pixelsVec.size();
00530   for (int i=0;  i<isize; ++i) {
00531     if ( (pixelsVec[i].x) == imin )
00532       q1 += pixelsVec[i].adc;
00533     else if ( (pixelsVec[i].x) == imax) 
00534       q2 += float(pixelsVec[i].adc);
00535     else 
00536       qm += float(pixelsVec[i].adc);
00537   }
00538   return;
00539 } 
00540 //-----------------------------------------------------------------------------
00541 //  Sum the pixels in the first and the last column, and the total.  Returns
00542 //  a vector of three elements with q_first, q_last and q_total.
00543 //  &&& Really need a simpler & cleaner way, this is very confusing...
00544 //-----------------------------------------------------------------------------
00545 void PixelCPEBase::yCharge(const vector<SiPixelCluster::Pixel>& pixelsVec,
00546                            const int& imin, const int& imax,
00547                            float& q1, float& q2) const {
00548   
00549   //calculate charge in the first and last pixel in y
00550   // and the inner cluster charge
00551   q1 = 0;
00552   q2 = 0;
00553   float qm=0;
00554   int isize = pixelsVec.size();
00555   for (int i=0;  i<isize; ++i) {
00556     if ( (pixelsVec[i].y) == imin) 
00557       q1 += pixelsVec[i].adc;
00558     else if ( (pixelsVec[i].y) == imax) 
00559       q2 += pixelsVec[i].adc;
00560     //else if (pixelsVec[i].y < ymax && pixelsVec[i].y > ymin ) 
00561     else  
00562       qm += float(pixelsVec[i].adc);
00563   }
00564   return;
00565 } 
00566 //-----------------------------------------------------------------------------
00567 //  Drift direction.
00568 //  Works OK for barrel and forward.
00569 //  The formulas used for dir_x,y,z have to be exactly the same as the ones
00570 //  used in the digitizer (SiPixelDigitizerAlgorithm.cc).
00571 //  Assumption: setTheDet() has been called already.
00572 //
00573 //  Petar (2/23/07): uhm, actually, there is a bug in the sign for both X and Y!
00574 //  (The signs have been fixed in SiPixelDigitizer, but not in here.)
00575 //-----------------------------------------------------------------------------
00576 LocalVector 
00577 PixelCPEBase::driftDirection( GlobalVector bfield ) const {
00578 
00579         Frame detFrame(theDet->surface().position(), theDet->surface().rotation());
00580         LocalVector Bfield = detFrame.toLocal(bfield);
00581 
00582         if(lorentzAngle_ == 0){
00583                 throw cms::Exception("invalidPointer") << "[PixelCPEBase::driftDirection] zero pointer to lorentz angle record ";
00584         }
00585         double langle = lorentzAngle_->getLorentzAngle(theDet->geographicalId().rawId());
00586         float alpha2;
00587         if (alpha2Order) {
00588                 alpha2 = langle*langle;
00589         } else {
00590                 alpha2 = 0.0;
00591         }
00592         // &&& dir_x should have a "-" and dir_y a "+"
00593         float dir_x =  ( langle * Bfield.y() + alpha2* Bfield.z()* Bfield.x() );
00594         float dir_y = -( langle * Bfield.x() - alpha2* Bfield.z()* Bfield.y() );
00595         float dir_z = -( 1 + alpha2* Bfield.z()*Bfield.z() );
00596         float scale = (1 + alpha2* Bfield.z()*Bfield.z() );
00597         LocalVector theDriftDirection = LocalVector(dir_x/scale, dir_y/scale, dir_z/scale );
00598         if ( theVerboseLevel > 9 ) 
00599                 LogDebug("PixelCPEBase") << " The drift direction in local coordinate is " 
00600                                         << theDriftDirection    ;
00601         
00602         return theDriftDirection;
00603 }
00604 
00605 //-----------------------------------------------------------------------------
00606 //  One-shot computation of the driftDirection and both lorentz shifts
00607 //-----------------------------------------------------------------------------
00608 void
00609 PixelCPEBase::computeLorentzShifts() const 
00610 {
00611   Frame detFrame(theDet->surface().position(), theDet->surface().rotation());
00612   GlobalVector global_Bfield = magfield_->inTesla( theDet->surface().position() );
00613   LocalVector  Bfield        = detFrame.toLocal(global_Bfield);
00614   if(lorentzAngle_ == 0){
00615         throw cms::Exception("invalidPointer") << "[PixelCPEBase::computeLorentzShifts] zero pointer to lorentz angle record ";
00616         }
00617   double langle = lorentzAngle_->getLorentzAngle(theDet->geographicalId().rawId());
00618   double alpha2;
00619   if ( alpha2Order) {
00620     alpha2 = langle * langle;
00621   }
00622   else  {
00623     alpha2 = 0.0;
00624   }
00625 
00626   // **********************************************************************
00627   // Our convention is the following:
00628   // +x is defined by the direction of the Lorentz drift!
00629   // +z is defined by the direction of E field (so electrons always go into -z!)
00630   // +y is defined by +x and +z, and it turns out to be always opposite to the +B field.
00631   // **********************************************************************
00632       
00633   // Note correct signs for dir_x and dir_y!
00634   double dir_x = -( langle * Bfield.y() + alpha2* Bfield.z()* Bfield.x() );
00635   double dir_y =  ( langle * Bfield.x() - alpha2* Bfield.z()* Bfield.y() );
00636   double dir_z = -( 1                                       + alpha2* Bfield.z()* Bfield.z() );
00637 
00638   // &&& Why do we need to scale???
00639   //double scale = (1 + alpha2* Bfield.z()*Bfield.z() );
00640   double scale = fabs( dir_z );  // same as 1 + alpha2*Bfield.z()*Bfield.z()
00641   driftDirection_ = LocalVector(dir_x/scale, dir_y/scale, dir_z/scale );  // last is -1 !
00642 
00643   // Max shift (at the other side of the sensor) in cm 
00644   lorentzShiftInCmX_ = driftDirection_.x()/driftDirection_.z() * theThickness;  // &&& redundant
00645   // Express the shift in units of pitch, 
00646   lorentzShiftX_ = lorentzShiftInCmX_ / thePitchX ; 
00647    
00648   // Max shift (at the other side of the sensor) in cm 
00649   lorentzShiftInCmY_ = driftDirection_.y()/driftDirection_.z() * theThickness;  // &&& redundant
00650   // Express the shift in units of pitch, 
00651   lorentzShiftY_ = lorentzShiftInCmY_ / thePitchY;
00652 
00653 
00654   if ( theVerboseLevel > 9 ) {
00655     LogDebug("PixelCPEBase") << " The drift direction in local coordinate is " 
00656                              << driftDirection_    ;
00657     
00658 //     cout << "Lorentz Drift (in cm) along X = " << lorentzShiftInCmX_ << endl;
00659 //     cout << "Lorentz Drift (in cm) along Y = " << lorentzShiftInCmY_ << endl;
00660   }
00661 }
00662 
00663 //-----------------------------------------------------------------------------
00668 //-----------------------------------------------------------------------------
00669 SiPixelRecHitQuality::QualWordType 
00670 PixelCPEBase::rawQualityWord() const
00671 {
00672   SiPixelRecHitQuality::QualWordType qualWord(0);
00673   
00674   SiPixelRecHitQuality::thePacking.setProbabilityXY ( probabilityXY() ,
00675                                                       qualWord );
00676   
00677   SiPixelRecHitQuality::thePacking.setProbabilityQ  ( probabilityQ_ , 
00678                                                       qualWord );
00679   
00680   SiPixelRecHitQuality::thePacking.setQBin          ( (int)qBin_, 
00681                                                       qualWord );
00682   
00683   SiPixelRecHitQuality::thePacking.setIsOnEdge      ( isOnEdge_,
00684                                                       qualWord );
00685   
00686   SiPixelRecHitQuality::thePacking.setHasBadPixels  ( hasBadPixels_,
00687                                                       qualWord );
00688   
00689   SiPixelRecHitQuality::thePacking.setSpansTwoROCs  ( spansTwoROCs_,
00690                                                       qualWord );
00691   
00692   SiPixelRecHitQuality::thePacking.setHasFilledProb ( hasFilledProb_,
00693                                                       qualWord );
00694   
00695   return qualWord;
00696 }