Reference Manual

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