CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DeDxTools.cc
Go to the documentation of this file.
1 #include "RecoTracker/DeDx/interface/DeDxTools.h"
2 #include <vector>
3 
4 #include <numeric>
5 namespace DeDxTools {
6 using namespace std;
7 using namespace reco;
8 
9 
10  void trajectoryRawHits(const edm::Ref<std::vector<Trajectory> >& trajectory, vector<RawHits>& hits, bool usePixel, bool useStrip)
11  {
12 
13  // vector<RawHits> hits;
14 
15  const vector<TrajectoryMeasurement> & measurements = trajectory->measurements();
16  for(vector<TrajectoryMeasurement>::const_iterator it = measurements.begin(); it!=measurements.end(); it++){
17 
18  //FIXME: check that "updated" State is the best one (wrt state in the middle)
19  TrajectoryStateOnSurface trajState=it->updatedState();
20  if( !trajState.isValid()) continue;
21 
22  const TrackingRecHit * recHit=(*it->recHit()).hit();
23 
24  LocalVector trackDirection = trajState.localDirection();
25  double cosine = trackDirection.z()/trackDirection.mag();
26 
27  if(const SiStripMatchedRecHit2D* matchedHit=dynamic_cast<const SiStripMatchedRecHit2D*>(recHit)){
28  if(!useStrip) continue;
29 
30  RawHits mono,stereo;
31  mono.trajectoryMeasurement = &(*it);
32  stereo.trajectoryMeasurement = &(*it);
33  mono.angleCosine = cosine;
34  stereo.angleCosine = cosine;
35  const std::vector<uint8_t> & amplitudes = matchedHit->monoCluster().amplitudes();
36  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
37  mono.NSaturating =0;
38  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
39 
40  const std::vector<uint8_t> & amplitudesSt = matchedHit->stereoCluster().amplitudes();
41  stereo.charge = accumulate(amplitudesSt.begin(), amplitudesSt.end(), 0);
42  stereo.NSaturating =0;
43  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)stereo.NSaturating++;}
44 
45  mono.detId= matchedHit->monoId();
46  stereo.detId= matchedHit->stereoId();
47 
48  hits.push_back(mono);
49  hits.push_back(stereo);
50 
51  }else if(const ProjectedSiStripRecHit2D* projectedHit=dynamic_cast<const ProjectedSiStripRecHit2D*>(recHit)) {
52  if(!useStrip) continue;
53 
54  const SiStripRecHit2D* singleHit=&(projectedHit->originalHit());
55  RawHits mono;
56 
57  mono.trajectoryMeasurement = &(*it);
58 
59  mono.angleCosine = cosine;
60  const std::vector<uint8_t> & amplitudes = singleHit->cluster()->amplitudes();
61  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
62  mono.NSaturating =0;
63  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
64 
65  mono.detId= singleHit->geographicalId();
66  hits.push_back(mono);
67 
68  }else if(const SiStripRecHit2D* singleHit=dynamic_cast<const SiStripRecHit2D*>(recHit)){
69  if(!useStrip) continue;
70 
71  RawHits mono;
72 
73  mono.trajectoryMeasurement = &(*it);
74 
75  mono.angleCosine = cosine;
76  const std::vector<uint8_t> & amplitudes = singleHit->cluster()->amplitudes();
77  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
78  mono.NSaturating =0;
79  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
80 
81  mono.detId= singleHit->geographicalId();
82  hits.push_back(mono);
83 
84  }else if(const SiStripRecHit1D* single1DHit=dynamic_cast<const SiStripRecHit1D*>(recHit)){
85  if(!useStrip) continue;
86 
87  RawHits mono;
88 
89  mono.trajectoryMeasurement = &(*it);
90 
91  mono.angleCosine = cosine;
92  const std::vector<uint8_t> & amplitudes = single1DHit->cluster()->amplitudes();
93  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
94  mono.NSaturating =0;
95  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
96 
97  mono.detId= single1DHit->geographicalId();
98  hits.push_back(mono);
99 
100 
101  }else if(const SiPixelRecHit* pixelHit=dynamic_cast<const SiPixelRecHit*>(recHit)){
102  if(!usePixel) continue;
103 
104  RawHits pixel;
105 
106  pixel.trajectoryMeasurement = &(*it);
107 
108  pixel.angleCosine = cosine;
109  pixel.charge = pixelHit->cluster()->charge();
110  pixel.NSaturating=-1;
111  pixel.detId= pixelHit->geographicalId();
112  hits.push_back(pixel);
113  }
114 
115  }
116  // return hits;
117  }
118 
119 
120 
121 
122 double genericAverage(const reco::DeDxHitCollection &hits, float expo )
123 {
124  double result=0;
125  size_t n = hits.size();
126  for(size_t i = 0; i< n; i ++)
127  {
128  result+=pow(hits[i].charge(),expo);
129  }
130  return (n>0)?pow(result/n,1./expo):0.;
131 }
132 
133 
134 bool shapeSelection(const std::vector<uint8_t> & ampls)
135 {
136  // ---------------- COMPTAGE DU NOMBRE DE MAXIMA --------------------------
137  //----------------------------------------------------------------------------
138 // printf("ShapeTest \n");
139  Int_t NofMax=0; Int_t recur255=1; Int_t recur254=1;
140  bool MaxOnStart=false;bool MaxInMiddle=false, MaxOnEnd =false;
141  Int_t MaxPos=0;
142  // Début avec max
143  if(ampls.size()!=1 && ((ampls[0]>ampls[1])
144  || (ampls.size()>2 && ampls[0]==ampls[1] && ampls[1]>ampls[2] && ampls[0]!=254 && ampls[0]!=255)
145  || (ampls.size()==2 && ampls[0]==ampls[1] && ampls[0]!=254 && ampls[0]!=255)) ){
146  NofMax=NofMax+1; MaxOnStart=true; }
147 
148  // Maximum entouré
149  if(ampls.size()>2){
150  for (unsigned int i =1; i < ampls.size()-1; i++) {
151  if( (ampls[i]>ampls[i-1] && ampls[i]>ampls[i+1])
152  || (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) ){
153  NofMax=NofMax+1; MaxInMiddle=true; MaxPos=i;
154  }
155  if(ampls[i]==255 && ampls[i]==ampls[i-1]) {
156  recur255=recur255+1;
157  MaxPos=i-(recur255/2);
158  if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
159  }
160  if(ampls[i]==254 && ampls[i]==ampls[i-1]) {
161  recur254=recur254+1;
162  MaxPos=i-(recur254/2);
163  if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
164  }
165  }
166  }
167  // Fin avec un max
168  if(ampls.size()>1){
169  if(ampls[ampls.size()-1]>ampls[ampls.size()-2]
170  || (ampls.size()>2 && ampls[ampls.size()-1]==ampls[ampls.size()-2] && ampls[ampls.size()-2]>ampls[ampls.size()-3] )
171  || ampls[ampls.size()-1]==255){
172  NofMax=NofMax+1; MaxOnEnd=true; }
173  }
174  // Si une seule strip touchée
175  if(ampls.size()==1){ NofMax=1;}
176 
177 
178 
179 
180  // --- SELECTION EN FONCTION DE LA FORME POUR LES MAXIMA UNIQUES ---------
181  //------------------------------------------------------------------------
182  /*
183  ____
184  | |____
185  ____| | |
186  | | | |____
187  ____| | | | |
188  | | | | | |____
189  __|____|____|____|____|____|____|__
190  C_Mnn C_Mn C_M C_D C_Dn C_Dnn
191  */
192 // bool shapetest=true;
193  bool shapecdtn=false;
194 
195 // Float_t C_M; Float_t C_D; Float_t C_Mn; Float_t C_Dn; Float_t C_Mnn; Float_t C_Dnn;
196  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;
197  Int_t CDPos;
198  Float_t coeff1=1.7; Float_t coeff2=2.0;
199  Float_t coeffn=0.10; Float_t coeffnn=0.02; Float_t noise=4.0;
200 
201  if(NofMax==1){
202 
203  if(MaxOnStart==true){
204  C_M=(Float_t)ampls[0]; C_D=(Float_t)ampls[1];
205  if(ampls.size()<3) shapecdtn=true ;
206  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;}
207  else if(ampls.size()>3){ C_Dn=(Float_t)ampls[2]; C_Dnn=(Float_t)ampls[3] ;
208  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
209  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){
210  shapecdtn=true;}
211  }
212  }
213 
214  if(MaxOnEnd==true){
215  C_M=(Float_t)ampls[ampls.size()-1]; C_D=(Float_t)ampls[ampls.size()-2];
216  if(ampls.size()<3) shapecdtn=true ;
217  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;}
218  else if(ampls.size()>3){C_Dn=(Float_t)ampls[ampls.size()-3] ; C_Dnn=(Float_t)ampls[ampls.size()-4] ;
219  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
220  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){
221  shapecdtn=true;}
222  }
223  }
224 
225  if(MaxInMiddle==true){
226  C_M=(Float_t)ampls[MaxPos];
227  int LeftOfMaxPos=MaxPos-1;if(LeftOfMaxPos<=0)LeftOfMaxPos=0;
228  int RightOfMaxPos=MaxPos+1;if(RightOfMaxPos>=(int)ampls.size())RightOfMaxPos=ampls.size()-1;
229  //int after = RightOfMaxPos; int before = LeftOfMaxPos; if (after>=(int)ampls.size() || before<0) std::cout<<"invalid read MaxPos:"<<MaxPos <<"size:"<<ampls.size() <<std::endl;
230  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;}
231  if(C_Mn<coeff1*coeffn*C_M+coeff2*coeffnn*C_D+2*noise || C_M==255){
232  if(ampls.size()==3) shapecdtn=true ;
233  else if(ampls.size()>3){
234  if(CDPos>MaxPos){
235  if(ampls.size()-CDPos-1==0){
236  C_Dn=0; C_Dnn=0;
237  }
238  if(ampls.size()-CDPos-1==1){
239  C_Dn=(Float_t)ampls[CDPos+1];
240  C_Dnn=0;
241  }
242  if(ampls.size()-CDPos-1>1){
243  C_Dn=(Float_t)ampls[CDPos+1];
244  C_Dnn=(Float_t)ampls[CDPos+2];
245  }
246  if(MaxPos>=2){
247  C_Mnn=(Float_t)ampls[MaxPos-2];
248  }
249  else if(MaxPos<2) C_Mnn=0;
250  }
251  if(CDPos<MaxPos){
252  if(CDPos==0){
253  C_Dn=0; C_Dnn=0;
254  }
255  if(CDPos==1){
256  C_Dn=(Float_t)ampls[0];
257  C_Dnn=0;
258  }
259  if(CDPos>1){
260  C_Dn=(Float_t)ampls[CDPos-1];
261  C_Dnn=(Float_t)ampls[CDPos-2];
262  }
263  if(ampls.size()-LeftOfMaxPos>1 && MaxPos+2<(int)(ampls.size())-1){
264  C_Mnn=(Float_t)ampls[MaxPos+2];
265  }else C_Mnn=0;
266  }
267  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
268  && C_Mnn<=coeff1*coeffn*C_Mn+coeff2*coeffnn*C_M+2*noise
269  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise) {
270  shapecdtn=true;
271  }
272 
273  }
274  }
275  }
276  }
277  if(ampls.size()==1){shapecdtn=true;}
278 
279  return shapecdtn;
280 }
281 
282 
283 
284 
285 
286 
287 }
i
int i
Definition: DBlmapReader.cc:9
DeDxTools::trajectoryRawHits
void trajectoryRawHits(const edm::Ref< std::vector< Trajectory > > &trajectory, std::vector< RawHits > &hits, bool usePixel, bool useStrip)
Definition: DeDxTools.cc:10
TrajectoryStateOnSurface::localDirection
LocalVector localDirection() const
Definition: TrajectoryStateOnSurface.h:177
reco::DeDxHitCollection
std::vector< DeDxHit > DeDxHitCollection
Definition: DeDxHit.h:49
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
DeDxTools::RawHits::angleCosine
double angleCosine
Definition: DeDxTools.h:21
PV3DBase::z
T z() const
Definition: PV3DBase.h:64
query.result
tuple result
Definition: query.py:137
DeDxTools::shapeSelection
bool shapeSelection(const std::vector< uint8_t > &ampls)
Definition: DeDxTools.cc:134
DeDxTools::RawHits::trajectoryMeasurement
const TrajectoryMeasurement * trajectoryMeasurement
Definition: DeDxTools.h:23
math::LocalVector
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< float >, ROOT::Math::LocalCoordinateSystemTag > LocalVector
vector in local coordinate system
Definition: Vector3D.h:25
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
DeDxTools::genericAverage
double genericAverage(const reco::DeDxHitCollection &, float expo=1.)
Definition: DeDxTools.cc:122
SiPixelRecHit
Pixel Reconstructed Hit.
Generated for CMSSW Reference Manual by doxygen 1.8.5