Reference Manual

Go to the documentation of this file.
00001 #include "RecoTracker/DeDx/interface/DeDxTools.h"
00002 #include <vector>
00003 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00004 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00005 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
00006 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
00007 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
00008 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
00009 #include "DataFormats/TrackReco/interface/DeDxHit.h"
00010 
00011 #include <numeric>
00012 namespace DeDxTools {
00013 using namespace std;
00014 using namespace reco;
00015 
00016                    
00017   void trajectoryRawHits(const edm::Ref<std::vector<Trajectory> >& trajectory, vector<RawHits>& hits, bool usePixel, bool useStrip)
00018   {
00019 
00020     //    vector<RawHits> hits;
00021 
00022     const vector<TrajectoryMeasurement> & measurements = trajectory->measurements();
00023     for(vector<TrajectoryMeasurement>::const_iterator it = measurements.begin(); it!=measurements.end(); it++){
00024 
00025       //FIXME: check that "updated" State is the best one (wrt state in the middle) 
00026       TrajectoryStateOnSurface trajState=it->updatedState();
00027       if( !trajState.isValid()) continue;
00028      
00029       const TrackingRecHit * recHit=(*it->recHit()).hit();
00030 
00031        LocalVector trackDirection = trajState.localDirection();
00032        double cosine = trackDirection.z()/trackDirection.mag();
00033               
00034        if(const SiStripMatchedRecHit2D* matchedHit=dynamic_cast<const SiStripMatchedRecHit2D*>(recHit)){
00035            if(!useStrip) continue;
00036 
00037            RawHits mono,stereo; 
00038            mono.trajectoryMeasurement = &(*it);
00039            stereo.trajectoryMeasurement = &(*it);
00040            mono.angleCosine = cosine; 
00041            stereo.angleCosine = cosine;
00042            const std::vector<uint8_t> &  amplitudes = matchedHit->monoHit()->cluster()->amplitudes(); 
00043            mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
00044            mono.NSaturating =0;
00045            for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
00046        
00047            const std::vector<uint8_t> & amplitudesSt = matchedHit->stereoHit()->cluster()->amplitudes();
00048            stereo.charge = accumulate(amplitudesSt.begin(), amplitudesSt.end(), 0);
00049            stereo.NSaturating =0;
00050            for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)stereo.NSaturating++;}
00051    
00052            mono.detId= matchedHit->monoHit()->geographicalId();
00053            stereo.detId= matchedHit->stereoHit()->geographicalId();
00054 
00055            hits.push_back(mono);
00056            hits.push_back(stereo);
00057 
00058         }else if(const ProjectedSiStripRecHit2D* projectedHit=dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit)) {
00059            if(!useStrip) continue;
00060 
00061            const SiStripRecHit2D* singleHit=&(projectedHit->originalHit());
00062            RawHits mono;
00063 
00064            mono.trajectoryMeasurement = &(*it);
00065 
00066            mono.angleCosine = cosine; 
00067            const std::vector<uint8_t> & amplitudes = singleHit->cluster()->amplitudes();
00068            mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
00069            mono.NSaturating =0;
00070            for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
00071 
00072            mono.detId= singleHit->geographicalId();
00073            hits.push_back(mono);
00074       
00075         }else if(const SiStripRecHit2D* singleHit=dynamic_cast<const SiStripRecHit2D*>(recHit)){
00076            if(!useStrip) continue;
00077 
00078            RawHits mono;
00079                
00080            mono.trajectoryMeasurement = &(*it);
00081 
00082            mono.angleCosine = cosine; 
00083            const std::vector<uint8_t> & amplitudes = singleHit->cluster()->amplitudes();
00084            mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
00085            mono.NSaturating =0;
00086            for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
00087 
00088            mono.detId= singleHit->geographicalId();
00089            hits.push_back(mono);
00090 
00091         }else if(const SiStripRecHit1D* single1DHit=dynamic_cast<const SiStripRecHit1D*>(recHit)){
00092            if(!useStrip) continue;
00093 
00094            RawHits mono;
00095 
00096            mono.trajectoryMeasurement = &(*it);
00097 
00098            mono.angleCosine = cosine;
00099            const std::vector<uint8_t> & amplitudes = single1DHit->cluster()->amplitudes();
00100            mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
00101            mono.NSaturating =0;
00102            for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
00103 
00104            mono.detId= single1DHit->geographicalId();
00105            hits.push_back(mono);
00106 
00107       
00108         }else if(const SiPixelRecHit* pixelHit=dynamic_cast<const SiPixelRecHit*>(recHit)){
00109            if(!usePixel) continue;
00110 
00111            RawHits pixel;
00112 
00113            pixel.trajectoryMeasurement = &(*it);
00114 
00115            pixel.angleCosine = cosine; 
00116            pixel.charge = pixelHit->cluster()->charge();
00117            pixel.NSaturating=-1;
00118            pixel.detId= pixelHit->geographicalId();
00119            hits.push_back(pixel);
00120        }
00121        
00122     }
00123     // return hits;
00124   }
00125 
00126 
00127 
00128 
00129 double genericAverage(const reco::DeDxHitCollection &hits, float expo )
00130 {
00131  double result=0;
00132  size_t n = hits.size();
00133  for(size_t i = 0; i< n; i ++)
00134  {
00135     result+=pow(hits[i].charge(),expo); 
00136  }
00137  return (n>0)?pow(result/n,1./expo):0.;
00138 }
00139 
00140 
00141 bool shapeSelection(const std::vector<uint8_t> & ampls)
00142 {
00143   // ----------------  COMPTAGE DU NOMBRE DE MAXIMA   --------------------------
00144   //----------------------------------------------------------------------------
00145 //      printf("ShapeTest \n");
00146          Int_t NofMax=0; Int_t recur255=1; Int_t recur254=1;
00147          bool MaxOnStart=false;bool MaxInMiddle=false, MaxOnEnd =false;
00148          Int_t MaxPos=0;
00149         // Début avec max
00150          if(ampls.size()!=1 && ((ampls[0]>ampls[1])
00151             || (ampls.size()>2 && ampls[0]==ampls[1] && ampls[1]>ampls[2] && ampls[0]!=254 && ampls[0]!=255) 
00152             || (ampls.size()==2 && ampls[0]==ampls[1] && ampls[0]!=254 && ampls[0]!=255)) ){
00153           NofMax=NofMax+1;  MaxOnStart=true;  }
00154 
00155         // Maximum entouré
00156          if(ampls.size()>2){
00157           for (unsigned int i =1; i < ampls.size()-1; i++) {
00158                 if( (ampls[i]>ampls[i-1] && ampls[i]>ampls[i+1]) 
00159                     || (ampls.size()>3 && i>0 && i<ampls.size()-2 && ampls[i]==ampls[i+1] && ampls[i]>ampls[i-1] && ampls[i]>ampls[i+2] && ampls[i]!=254 && ampls[i]!=255) ){ 
00160                  NofMax=NofMax+1; MaxInMiddle=true;  MaxPos=i; 
00161                 }
00162                 if(ampls[i]==255 && ampls[i]==ampls[i-1]) {
00163                         recur255=recur255+1;
00164                         MaxPos=i-(recur255/2);
00165                         if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
00166                 }
00167                 if(ampls[i]==254 && ampls[i]==ampls[i-1]) {
00168                         recur254=recur254+1;
00169                         MaxPos=i-(recur254/2);
00170                         if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
00171                 }
00172             }
00173          }
00174         // Fin avec un max
00175          if(ampls.size()>1){
00176           if(ampls[ampls.size()-1]>ampls[ampls.size()-2]
00177              || (ampls.size()>2 && ampls[ampls.size()-1]==ampls[ampls.size()-2] && ampls[ampls.size()-2]>ampls[ampls.size()-3] ) 
00178              ||  ampls[ampls.size()-1]==255){
00179            NofMax=NofMax+1;  MaxOnEnd=true;   }
00180          }
00181         // Si une seule strip touchée
00182         if(ampls.size()==1){    NofMax=1;}
00183 
00184 
00185 
00186 
00187   // ---  SELECTION EN FONCTION DE LA FORME POUR LES MAXIMA UNIQUES ---------
00188   //------------------------------------------------------------------------
00189   /*
00190                ____
00191               |    |____
00192           ____|    |    |
00193          |    |    |    |____
00194      ____|    |    |    |    |
00195     |    |    |    |    |    |____
00196   __|____|____|____|____|____|____|__
00197     C_Mnn C_Mn C_M  C_D  C_Dn C_Dnn
00198   */
00199 //   bool shapetest=true;
00200    bool shapecdtn=false;
00201 
00202 //      Float_t C_M;    Float_t C_D;    Float_t C_Mn;   Float_t C_Dn;   Float_t C_Mnn;  Float_t C_Dnn;
00203         Float_t C_M=0.0;        Float_t C_D=0.0;        Float_t C_Mn=10000;     Float_t C_Dn=10000;     Float_t C_Mnn=10000;    Float_t C_Dnn=10000;
00204         Int_t CDPos;
00205         Float_t coeff1=1.7;     Float_t coeff2=2.0;
00206         Float_t coeffn=0.10;    Float_t coeffnn=0.02; Float_t noise=4.0;
00207 
00208         if(NofMax==1){
00209 
00210                 if(MaxOnStart==true){
00211                         C_M=(Float_t)ampls[0]; C_D=(Float_t)ampls[1];
00212                                 if(ampls.size()<3) shapecdtn=true ;
00213                                 else if(ampls.size()==3){C_Dn=(Float_t)ampls[2] ; if(C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255) shapecdtn=true;}
00214                                 else if(ampls.size()>3){ C_Dn=(Float_t)ampls[2];  C_Dnn=(Float_t)ampls[3] ;
00215                                                         if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
00216                                                            && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){
00217                                                          shapecdtn=true;}
00218                                 }
00219                 }
00220 
00221                 if(MaxOnEnd==true){
00222                         C_M=(Float_t)ampls[ampls.size()-1]; C_D=(Float_t)ampls[ampls.size()-2];
00223                                 if(ampls.size()<3) shapecdtn=true ;
00224                                 else if(ampls.size()==3){C_Dn=(Float_t)ampls[0] ; if(C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255) shapecdtn=true;}
00225                                 else if(ampls.size()>3){C_Dn=(Float_t)ampls[ampls.size()-3] ; C_Dnn=(Float_t)ampls[ampls.size()-4] ; 
00226                                                         if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
00227                                                            && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){ 
00228                                                          shapecdtn=true;}
00229                                 }
00230                 }
00231 
00232                 if(MaxInMiddle==true){
00233                         C_M=(Float_t)ampls[MaxPos];
00234                         int LeftOfMaxPos=MaxPos-1;if(LeftOfMaxPos<=0)LeftOfMaxPos=0;
00235                         int RightOfMaxPos=MaxPos+1;if(RightOfMaxPos>=(int)ampls.size())RightOfMaxPos=ampls.size()-1;
00236                         //int after = RightOfMaxPos; int before = LeftOfMaxPos; if (after>=(int)ampls.size() ||  before<0)  std::cout<<"invalid read MaxPos:"<<MaxPos <<"size:"<<ampls.size() <<std::endl; 
00237                         if(ampls[LeftOfMaxPos]<ampls[RightOfMaxPos]){ C_D=(Float_t)ampls[RightOfMaxPos]; C_Mn=(Float_t)ampls[LeftOfMaxPos];CDPos=RightOfMaxPos;} else{ C_D=(Float_t)ampls[LeftOfMaxPos]; C_Mn=(Float_t)ampls[RightOfMaxPos];CDPos=LeftOfMaxPos;}
00238                         if(C_Mn<coeff1*coeffn*C_M+coeff2*coeffnn*C_D+2*noise || C_M==255){ 
00239                                 if(ampls.size()==3) shapecdtn=true ;
00240                                 else if(ampls.size()>3){
00241                                         if(CDPos>MaxPos){
00242                                                 if(ampls.size()-CDPos-1==0){
00243                                                         C_Dn=0; C_Dnn=0;
00244                                                 }
00245                                                 if(ampls.size()-CDPos-1==1){
00246                                                         C_Dn=(Float_t)ampls[CDPos+1];
00247                                                         C_Dnn=0;
00248                                                 }
00249                                                 if(ampls.size()-CDPos-1>1){
00250                                                         C_Dn=(Float_t)ampls[CDPos+1];
00251                                                         C_Dnn=(Float_t)ampls[CDPos+2];
00252                                                 }
00253                                                 if(MaxPos>=2){
00254                                                         C_Mnn=(Float_t)ampls[MaxPos-2];
00255                                                 }
00256                                                 else if(MaxPos<2) C_Mnn=0;
00257                                         }
00258                                         if(CDPos<MaxPos){
00259                                                 if(CDPos==0){
00260                                                         C_Dn=0; C_Dnn=0;
00261                                                 }
00262                                                 if(CDPos==1){
00263                                                         C_Dn=(Float_t)ampls[0];
00264                                                         C_Dnn=0;
00265                                                 }
00266                                                 if(CDPos>1){
00267                                                         C_Dn=(Float_t)ampls[CDPos-1];
00268                                                         C_Dnn=(Float_t)ampls[CDPos-2];
00269                                                 }
00270                                                 if(ampls.size()-LeftOfMaxPos>1 && MaxPos+2<(int)(ampls.size())-1){
00271                                                         C_Mnn=(Float_t)ampls[MaxPos+2];
00272                                                 }else C_Mnn=0;                                                  
00273                                         }
00274                                         if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
00275                                            && C_Mnn<=coeff1*coeffn*C_Mn+coeff2*coeffnn*C_M+2*noise
00276                                            && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise) {
00277                                                 shapecdtn=true;
00278                                         }
00279 
00280                                 }
00281                         }                       
00282                 }
00283         }
00284         if(ampls.size()==1){shapecdtn=true;}
00285 
00286    return shapecdtn;
00287 } 
00288 
00289 
00290 
00291 
00292 
00293 
00294 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/RecoTracker/DeDx/src/DeDxTools.cc
