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