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  RawHits mono;
55 
56  mono.trajectoryMeasurement = &(*it);
57 
58  mono.angleCosine = cosine;
59  const std::vector<uint8_t> & amplitudes = projectedHit->cluster()->amplitudes();
60  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
61  mono.NSaturating =0;
62  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
63 
64  mono.detId= projectedHit->originalId();
65  hits.push_back(mono);
66 
67  }else if(const SiStripRecHit2D* singleHit=dynamic_cast<const SiStripRecHit2D*>(recHit)){
68  if(!useStrip) continue;
69 
70  RawHits mono;
71 
72  mono.trajectoryMeasurement = &(*it);
73 
74  mono.angleCosine = cosine;
75  const std::vector<uint8_t> & amplitudes = singleHit->cluster()->amplitudes();
76  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
77  mono.NSaturating =0;
78  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
79 
80  mono.detId= singleHit->geographicalId();
81  hits.push_back(mono);
82 
83  }else if(const SiStripRecHit1D* single1DHit=dynamic_cast<const SiStripRecHit1D*>(recHit)){
84  if(!useStrip) continue;
85 
86  RawHits mono;
87 
88  mono.trajectoryMeasurement = &(*it);
89 
90  mono.angleCosine = cosine;
91  const std::vector<uint8_t> & amplitudes = single1DHit->cluster()->amplitudes();
92  mono.charge = accumulate(amplitudes.begin(), amplitudes.end(), 0);
93  mono.NSaturating =0;
94  for(unsigned int i=0;i<amplitudes.size();i++){if(amplitudes[i]>=254)mono.NSaturating++;}
95 
96  mono.detId= single1DHit->geographicalId();
97  hits.push_back(mono);
98 
99 
100  }else if(const SiPixelRecHit* pixelHit=dynamic_cast<const SiPixelRecHit*>(recHit)){
101  if(!usePixel) continue;
102 
103  RawHits pixel;
104 
105  pixel.trajectoryMeasurement = &(*it);
106 
107  pixel.angleCosine = cosine;
108  pixel.charge = pixelHit->cluster()->charge();
109  pixel.NSaturating=-1;
110  pixel.detId= pixelHit->geographicalId();
111  hits.push_back(pixel);
112  }
113 
114  }
115  // return hits;
116  }
117 
118 
119 
120 
121 double genericAverage(const reco::DeDxHitCollection &hits, float expo )
122 {
123  double result=0;
124  size_t n = hits.size();
125  for(size_t i = 0; i< n; i ++)
126  {
127  result+=pow(hits[i].charge(),expo);
128  }
129  return (n>0)?pow(result/n,1./expo):0.;
130 }
131 
132 
133 bool shapeSelection(const std::vector<uint8_t> & ampls)
134 {
135  // ---------------- COMPTAGE DU NOMBRE DE MAXIMA --------------------------
136  //----------------------------------------------------------------------------
137 // printf("ShapeTest \n");
138  Int_t NofMax=0; Int_t recur255=1; Int_t recur254=1;
139  bool MaxOnStart=false;bool MaxInMiddle=false, MaxOnEnd =false;
140  Int_t MaxPos=0;
141  // Début avec max
142  if(ampls.size()!=1 && ((ampls[0]>ampls[1])
143  || (ampls.size()>2 && ampls[0]==ampls[1] && ampls[1]>ampls[2] && ampls[0]!=254 && ampls[0]!=255)
144  || (ampls.size()==2 && ampls[0]==ampls[1] && ampls[0]!=254 && ampls[0]!=255)) ){
145  NofMax=NofMax+1; MaxOnStart=true; }
146 
147  // Maximum entouré
148  if(ampls.size()>2){
149  for (unsigned int i =1; i < ampls.size()-1; i++) {
150  if( (ampls[i]>ampls[i-1] && ampls[i]>ampls[i+1])
151  || (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) ){
152  NofMax=NofMax+1; MaxInMiddle=true; MaxPos=i;
153  }
154  if(ampls[i]==255 && ampls[i]==ampls[i-1]) {
155  recur255=recur255+1;
156  MaxPos=i-(recur255/2);
157  if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
158  }
159  if(ampls[i]==254 && ampls[i]==ampls[i-1]) {
160  recur254=recur254+1;
161  MaxPos=i-(recur254/2);
162  if(ampls[i]>ampls[i+1]){NofMax=NofMax+1;MaxInMiddle=true;}
163  }
164  }
165  }
166  // Fin avec un max
167  if(ampls.size()>1){
168  if(ampls[ampls.size()-1]>ampls[ampls.size()-2]
169  || (ampls.size()>2 && ampls[ampls.size()-1]==ampls[ampls.size()-2] && ampls[ampls.size()-2]>ampls[ampls.size()-3] )
170  || ampls[ampls.size()-1]==255){
171  NofMax=NofMax+1; MaxOnEnd=true; }
172  }
173  // Si une seule strip touchée
174  if(ampls.size()==1){ NofMax=1;}
175 
176 
177 
178 
179  // --- SELECTION EN FONCTION DE LA FORME POUR LES MAXIMA UNIQUES ---------
180  //------------------------------------------------------------------------
181  /*
182  ____
183  | |____
184  ____| | |
185  | | | |____
186  ____| | | | |
187  | | | | | |____
188  __|____|____|____|____|____|____|__
189  C_Mnn C_Mn C_M C_D C_Dn C_Dnn
190  */
191 // bool shapetest=true;
192  bool shapecdtn=false;
193 
194 // Float_t C_M; Float_t C_D; Float_t C_Mn; Float_t C_Dn; Float_t C_Mnn; Float_t C_Dnn;
195  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;
196  Int_t CDPos;
197  Float_t coeff1=1.7; Float_t coeff2=2.0;
198  Float_t coeffn=0.10; Float_t coeffnn=0.02; Float_t noise=4.0;
199 
200  if(NofMax==1){
201 
202  if(MaxOnStart==true){
203  C_M=(Float_t)ampls[0]; C_D=(Float_t)ampls[1];
204  if(ampls.size()<3) shapecdtn=true ;
205  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;}
206  else if(ampls.size()>3){ C_Dn=(Float_t)ampls[2]; C_Dnn=(Float_t)ampls[3] ;
207  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
208  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){
209  shapecdtn=true;}
210  }
211  }
212 
213  if(MaxOnEnd==true){
214  C_M=(Float_t)ampls[ampls.size()-1]; C_D=(Float_t)ampls[ampls.size()-2];
215  if(ampls.size()<3) shapecdtn=true ;
216  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;}
217  else if(ampls.size()>3){C_Dn=(Float_t)ampls[ampls.size()-3] ; C_Dnn=(Float_t)ampls[ampls.size()-4] ;
218  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
219  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise){
220  shapecdtn=true;}
221  }
222  }
223 
224  if(MaxInMiddle==true){
225  C_M=(Float_t)ampls[MaxPos];
226  int LeftOfMaxPos=MaxPos-1;if(LeftOfMaxPos<=0)LeftOfMaxPos=0;
227  int RightOfMaxPos=MaxPos+1;if(RightOfMaxPos>=(int)ampls.size())RightOfMaxPos=ampls.size()-1;
228  //int after = RightOfMaxPos; int before = LeftOfMaxPos; if (after>=(int)ampls.size() || before<0) std::cout<<"invalid read MaxPos:"<<MaxPos <<"size:"<<ampls.size() <<std::endl;
229  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;}
230  if(C_Mn<coeff1*coeffn*C_M+coeff2*coeffnn*C_D+2*noise || C_M==255){
231  if(ampls.size()==3) shapecdtn=true ;
232  else if(ampls.size()>3){
233  if(CDPos>MaxPos){
234  if(ampls.size()-CDPos-1==0){
235  C_Dn=0; C_Dnn=0;
236  }
237  if(ampls.size()-CDPos-1==1){
238  C_Dn=(Float_t)ampls[CDPos+1];
239  C_Dnn=0;
240  }
241  if(ampls.size()-CDPos-1>1){
242  C_Dn=(Float_t)ampls[CDPos+1];
243  C_Dnn=(Float_t)ampls[CDPos+2];
244  }
245  if(MaxPos>=2){
246  C_Mnn=(Float_t)ampls[MaxPos-2];
247  }
248  else if(MaxPos<2) C_Mnn=0;
249  }
250  if(CDPos<MaxPos){
251  if(CDPos==0){
252  C_Dn=0; C_Dnn=0;
253  }
254  if(CDPos==1){
255  C_Dn=(Float_t)ampls[0];
256  C_Dnn=0;
257  }
258  if(CDPos>1){
259  C_Dn=(Float_t)ampls[CDPos-1];
260  C_Dnn=(Float_t)ampls[CDPos-2];
261  }
262  if(ampls.size()-LeftOfMaxPos>1 && MaxPos+2<(int)(ampls.size())-1){
263  C_Mnn=(Float_t)ampls[MaxPos+2];
264  }else C_Mnn=0;
265  }
266  if((C_Dn<=coeff1*coeffn*C_D+coeff2*coeffnn*C_M+2*noise || C_D==255)
267  && C_Mnn<=coeff1*coeffn*C_Mn+coeff2*coeffnn*C_M+2*noise
268  && C_Dnn<=coeff1*coeffn*C_Dn+coeff2*coeffnn*C_D+2*noise) {
269  shapecdtn=true;
270  }
271 
272  }
273  }
274  }
275  }
276  if(ampls.size()==1){shapecdtn=true;}
277 
278  return shapecdtn;
279 }
280 
281 
282 
283 
284 
285 
286 }
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:122
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:107
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:133
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:121
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:23
Generated for CMSSW Reference Manual by doxygen 1.8.5