CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
DeDxDiscriminatorProducer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: DeDxDiscriminatorProducer
4 // Class: DeDxDiscriminatorProducer
5 //
13 //
14 // Original Author: andrea
15 // Created: Thu May 31 14:09:02 CEST 2007
16 // Code Updates: loic Quertenmont (querten)
17 // Created: Thu May 10 14:09:02 CEST 2008
18 // $Id: DeDxDiscriminatorProducer.cc,v 1.3 2013/02/27 13:28:30 muzaffar Exp $
19 //
20 //
21 
22 
23 // system include files
24 #include <memory>
27 
31 
35 
41 
42 
43 #include "TFile.h"
44 
45 using namespace reco;
46 using namespace std;
47 using namespace edm;
48 
50 {
51 
52  produces<ValueMap<DeDxData> >();
53 
54  m_tracksTag = iConfig.getParameter<edm::InputTag>("tracks");
55  m_trajTrackAssociationTag = iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation");
56 
57  usePixel = iConfig.getParameter<bool>("UsePixel");
58  useStrip = iConfig.getParameter<bool>("UseStrip");
59  if(!usePixel && !useStrip)
60  edm::LogWarning("DeDxHitsProducer") << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
61 
62  Formula = iConfig.getUntrackedParameter<unsigned>("Formula" , 0);
63  Reccord = iConfig.getUntrackedParameter<string> ("Reccord" , "SiStripDeDxMip_3D_Rcd");
64  ProbabilityMode = iConfig.getUntrackedParameter<string> ("ProbabilityMode" , "Accumulation");
65 
66 
67  MinTrackMomentum = iConfig.getUntrackedParameter<double> ("minTrackMomentum" , 0.0);
68  MaxTrackMomentum = iConfig.getUntrackedParameter<double> ("maxTrackMomentum" , 99999.0);
69  MinTrackEta = iConfig.getUntrackedParameter<double> ("minTrackEta" , -5.0);
70  MaxTrackEta = iConfig.getUntrackedParameter<double> ("maxTrackEta" , 5.0);
71  MaxNrStrips = iConfig.getUntrackedParameter<unsigned>("maxNrStrips" , 255);
72  MinTrackHits = iConfig.getUntrackedParameter<unsigned>("MinTrackHits" , 3);
73 
74  shapetest = iConfig.getParameter<bool>("ShapeTest");
75  useCalibration = iConfig.getParameter<bool>("UseCalibration");
76  m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
77 
78  Prob_ChargePath = NULL;
79 }
80 
81 
83 
84 // ------------ method called once each job just before starting event loop ------------
86 {
88  if( strcmp(Reccord.c_str(),"SiStripDeDxMip_3D_Rcd")==0){
89  iSetup.get<SiStripDeDxMip_3D_Rcd>() .get(DeDxMapHandle_);
90  }else if(strcmp(Reccord.c_str(),"SiStripDeDxPion_3D_Rcd")==0){
91  iSetup.get<SiStripDeDxPion_3D_Rcd>().get(DeDxMapHandle_);
92  }else if(strcmp(Reccord.c_str(),"SiStripDeDxKaon_3D_Rcd")==0){
93  iSetup.get<SiStripDeDxKaon_3D_Rcd>().get(DeDxMapHandle_);
94  }else if(strcmp(Reccord.c_str(),"SiStripDeDxProton_3D_Rcd")==0){
95  iSetup.get<SiStripDeDxProton_3D_Rcd>().get(DeDxMapHandle_);
96  }else if(strcmp(Reccord.c_str(),"SiStripDeDxElectron_3D_Rcd")==0){
97  iSetup.get<SiStripDeDxElectron_3D_Rcd>().get(DeDxMapHandle_);
98  }else{
99 // printf("The reccord %s is not known by the DeDxDiscriminatorProducer\n", Reccord.c_str());
100 // printf("Program will exit now\n");
101  exit(0);
102  }
103  DeDxMap_ = *DeDxMapHandle_.product();
104 
105  double xmin = DeDxMap_.rangeX().min;
106  double xmax = DeDxMap_.rangeX().max;
107  double ymin = DeDxMap_.rangeY().min;
108  double ymax = DeDxMap_.rangeY().max;
109  double zmin = DeDxMap_.rangeZ().min;
110  double zmax = DeDxMap_.rangeZ().max;
111 
112  if(Prob_ChargePath)delete Prob_ChargePath;
113  Prob_ChargePath = new TH3D ("Prob_ChargePath" , "Prob_ChargePath" , DeDxMap_.numberOfBinsX(), xmin, xmax, DeDxMap_.numberOfBinsY() , ymin, ymax, DeDxMap_.numberOfBinsZ(), zmin, zmax);
114 
115 
116 
117  if(strcmp(ProbabilityMode.c_str(),"Accumulation")==0){
118 // printf("LOOOP ON P\n");
119  for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
120 // printf("LOOOP ON PATH\n");
121  for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
122 // printf("LOOOP ON CHARGE\n");
123 
124  double Ni = 0;
125  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);}
126 
127  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
128  double tmp = 0;
129  for(int l=0;l<=k;l++){ tmp+=DeDxMap_.binContent(i,j,l);}
130 
131  if(Ni>0){
132  Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
133 // printf("P=%6.2f Path=%6.2f Charge%8.2f --> Prob=%8.3f\n",Prob_ChargePath->GetXaxis()->GetBinCenter(i), Prob_ChargePath->GetYaxis()->GetBinCenter(j), Prob_ChargePath->GetZaxis()->GetBinCenter(k), tmp/Ni);
134  }else{
135  Prob_ChargePath->SetBinContent (i, j, k, 0);
136  }
137  }
138  }
139  }
140  }else if(strcmp(ProbabilityMode.c_str(),"Integral")==0){
141 // printf("LOOOP ON P\n");
142  for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
143 // printf("LOOOP ON PATH\n");
144  for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
145 // printf("LOOOP ON CHARGE\n");
146 
147  double Ni = 0;
148  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);}
149 
150  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
151  double tmp = DeDxMap_.binContent(i,j,k);
152 
153  if(Ni>0){
154  Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
155 // printf("P=%6.2f Path=%6.2f Charge%8.2f --> Prob=%8.3f\n",Prob_ChargePath->GetXaxis()->GetBinCenter(i), Prob_ChargePath->GetYaxis()->GetBinCenter(j), Prob_ChargePath->GetZaxis()->GetBinCenter(k), tmp/Ni);
156  }else{
157  Prob_ChargePath->SetBinContent (i, j, k, 0);
158  }
159  }
160  }
161  }
162  }else{
163 // printf("The ProbabilityMode: %s is unknown\n",ProbabilityMode.c_str());
164 // printf("The program will stop now\n");
165  exit(0);
166  }
167 
168 
169 
170 /*
171  for(int i=0;i<Prob_ChargePath->GetXaxis()->GetNbins();i++){
172  for(int j=0;j<Prob_ChargePath->GetYaxis()->GetNbins();j++){
173  double tmp = DeDxMap_.binContent(i,j);
174  Prob_ChargePath->SetBinContent (i, j, tmp);
175  printf("%i %i --> %f\n",i,j,tmp);
176  }
177  }
178 */
179 
180 
181 
182  if(MODsColl.size()==0){
184  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
185  m_tracker = tkGeom.product();
186 
187  vector<GeomDet*> Det = tkGeom->dets();
188  for(unsigned int i=0;i<Det.size();i++){
189  DetId Detid = Det[i]->geographicalId();
190  int SubDet = Detid.subdetId();
191 
192  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
193  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
194 
195  StripGeomDetUnit* DetUnit = dynamic_cast<StripGeomDetUnit*> (Det[i]);
196  if(!DetUnit)continue;
197 
198  const StripTopology& Topo = DetUnit->specificTopology();
199  unsigned int NAPV = Topo.nstrips()/128;
200 
201  double Eta = DetUnit->position().basicVector().eta();
202  double R = DetUnit->position().basicVector().transverse();
203  double Thick = DetUnit->surface().bounds().thickness();
204 
205  stModInfo* MOD = new stModInfo;
206  MOD->DetId = Detid.rawId();
207  MOD->SubDet = SubDet;
208  MOD->Eta = Eta;
209  MOD->R = R;
210  MOD->Thickness = Thick;
211  MOD->NAPV = NAPV;
212  MODsColl[MOD->DetId] = MOD;
213  }
214  }
215 
216  MakeCalibrationMap();
217  }
218 
219 }
220 
222 {
223  MODsColl.clear();
224 }
225 
226 
227 
229 {
230  auto_ptr<ValueMap<DeDxData> > trackDeDxDiscrimAssociation(new ValueMap<DeDxData> );
231  ValueMap<DeDxData>::Filler filler(*trackDeDxDiscrimAssociation);
232 
233  Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
234  iEvent.getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
235  const TrajTrackAssociationCollection TrajToTrackMap = *trajTrackAssociationHandle.product();
236 
237  edm::Handle<reco::TrackCollection> trackCollectionHandle;
238  iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
239 
241  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
242  m_tracker = tkGeom.product();
243 
244  std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.size() );
245 
246  unsigned track_index = 0;
247  for(TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap.begin(); it!=TrajToTrackMap.end(); ++it, track_index++) {
248  dEdxDiscrims[track_index] = DeDxData(-1, -2, 0 );
249 
250  const Track track = *it->val;
251  const Trajectory traj = *it->key;
252 
253  if(track.eta() <MinTrackEta || track.eta()>MaxTrackEta ){continue;}
254  if(track.p() <MinTrackMomentum || track.p() >MaxTrackMomentum){continue;}
255  if(track.found()<MinTrackHits ){continue;}
256 
257  std::vector<double> vect_probs;
258  vector<TrajectoryMeasurement> measurements = traj.measurements();
259 
260  unsigned int NClusterSaturating = 0;
261  for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
262 
263  TrajectoryStateOnSurface trajState = measurement_it->updatedState();
264  if( !trajState.isValid() ) continue;
265 
266  const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
267  const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
268  const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
269  const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
270 
271  double Prob;
272  if(sistripsimplehit){
273  Prob = GetProbability(DeDxTools::GetCluster(sistripsimplehit), trajState,sistripsimplehit->geographicalId());
274  if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripsimplehit)->amplitudes()))) Prob=-1.0;
275  if(Prob>=0) vect_probs.push_back(Prob);
276 
277  if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripsimplehit),sistripsimplehit->geographicalId())>0)NClusterSaturating++;
278 
279  }else if(sistripmatchedhit){
280  Prob = GetProbability(DeDxTools::GetCluster(sistripmatchedhit->monoHit()), trajState, sistripmatchedhit->monoId());
281  if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripmatchedhit->monoHit())->amplitudes()))) Prob=-1.0;
282  if(Prob>=0) vect_probs.push_back(Prob);
283 
284  Prob = GetProbability(DeDxTools::GetCluster(sistripmatchedhit->stereoHit()), trajState,sistripmatchedhit->stereoId());
285  if(Prob>=0) vect_probs.push_back(Prob);
286 
287  if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripmatchedhit->monoHit()),sistripmatchedhit->monoId()) >0)NClusterSaturating++;
288  if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripmatchedhit->stereoHit()),sistripmatchedhit->stereoId())>0)NClusterSaturating++;
289  }else if(sistripsimple1dhit){
290  Prob = GetProbability(DeDxTools::GetCluster(sistripsimple1dhit), trajState, sistripsimple1dhit->geographicalId());
291  if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripsimple1dhit)->amplitudes()))) Prob=-1.0;
292  if(Prob>=0) vect_probs.push_back(Prob);
293  if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripsimple1dhit),sistripsimple1dhit->geographicalId())>0)NClusterSaturating++;
294  }else{
295  }
296  }
297 
298  double estimator = ComputeDiscriminator(vect_probs);
299  int size = vect_probs.size();
300  float Error = -1;
301 
302  //WARNING: Since the dEdX Error is not properly computed for the moment
303  //It was decided to store the number of saturating cluster in that dataformat
304  Error = NClusterSaturating;
305  dEdxDiscrims[track_index] = DeDxData(estimator, Error, size );
306 
307 // printf("%i --> %g\n",size,estimator);
308  }
309 
310  filler.insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
311  filler.fill();
312  iEvent.put(trackDeDxDiscrimAssociation);
313 }
314 
315 
317  const uint32_t & detId){
318  const vector<uint8_t>& ampls = cluster->amplitudes();
319 
320  stModInfo* MOD = NULL;
321  if(useCalibration)MOD = MODsColl[detId];
322 
323  int SaturatingStrip = 0;
324  for(unsigned int i=0;i<ampls.size();i++){
325  int StripCharge = ampls[i];
326  if(MOD){StripCharge = (int)(StripCharge / MOD->Gain);}
327  if(StripCharge>=254)SaturatingStrip++;
328  }
329  return SaturatingStrip;
330 }
331 
333  const uint32_t & detId)
334 {
335  // Get All needed variables
336  LocalVector trackDirection = trajState.localDirection();
337  double cosine = trackDirection.z()/trackDirection.mag();
338  const vector<uint8_t>& ampls = cluster->amplitudes();
339  // uint32_t detId = cluster->geographicalId();
340  // int firstStrip = cluster->firstStrip();
341  stModInfo* MOD = MODsColl[detId];
342 
343 
344  // Sanity Checks
345  if( ampls.size()>MaxNrStrips) {return -1;}
346 // if( DeDxDiscriminatorTools::IsSaturatingStrip (ampls)) {return -1;}
347 // if( DeDxDiscriminatorTools::IsSpanningOver2APV (firstStrip, ampls.size())) {return -1;}
348 // if(!DeDxDiscriminatorTools::IsFarFromBorder (trajState, m_tracker->idToDetUnit(DetId(detId)) )) {return -1;}
349 
350 
351  // Find Probability for this given Charge and Path
352  double charge = 0;
353  if(useCalibration){
354  for(unsigned int i=0;i<ampls.size();i++){
355  int CalibratedCharge = ampls[i];
356  CalibratedCharge = (int)(CalibratedCharge / MOD->Gain);
357  if(CalibratedCharge>=1024){
358  CalibratedCharge = 255;
359  }else if(CalibratedCharge>254){
360  CalibratedCharge = 254;
361  }
362  charge+=CalibratedCharge;
363  }
364  }else{
365  charge = DeDxDiscriminatorTools::charge(ampls);
366  }
367  double path = DeDxDiscriminatorTools::path(cosine,MOD->Thickness);
368 
369  int BinX = Prob_ChargePath->GetXaxis()->FindBin(trajState.localMomentum().mag());
370  int BinY = Prob_ChargePath->GetYaxis()->FindBin(path);
371  int BinZ = Prob_ChargePath->GetZaxis()->FindBin(charge/path);
372  return Prob_ChargePath->GetBinContent(BinX,BinY,BinZ);
373 }
374 
375 
376 double DeDxDiscriminatorProducer::ComputeDiscriminator(std::vector<double>& vect_probs)
377 {
378  double estimator = -1;
379  int size = vect_probs.size();
380  if(size<=0)return estimator;
381 
382  if(Formula==0){
383  double P = 1;
384  for(int i=0;i<size;i++){
385  if(vect_probs[i]<=0.0001){P *= pow(0.0001 , 1.0/size);}
386  else {P *= pow(vect_probs[i], 1.0/size);}
387  }
388  estimator = P;
389  }else if(Formula==1){
390  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
391  for(int i=0;i<size;i++){if(vect_probs[i]<=0.0001)vect_probs[i] = 0.0001; }
392 
393  double SumJet = 0.;
394  for(int i=0;i<size;i++){ SumJet+= log(vect_probs[i]); }
395 
396  double Loginvlog=log(-SumJet);
397  double Prob =1.;
398  double lfact=1.;
399 
400  for(int l=1; l!=size; l++){
401  lfact*=l;
402  Prob+=exp(l*Loginvlog-log(1.*lfact));
403  }
404 
405  double LogProb=log(Prob);
406  double ProbJet=std::min(exp(std::max(LogProb+SumJet,-30.)),1.);
407  estimator = -log10(ProbJet)/4.;
408  estimator = 1-estimator;
409  }else if(Formula==2){
410  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
411  double P = 1.0/(12*size);
412  for(int i=1;i<=size;i++){
413  P += pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
414  }
415  P *= (3.0/size);
416  estimator = P;
417  }else{
418  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
419  double P = 1.0/(12*size);
420  for(int i=1;i<=size;i++){
421  P += vect_probs[i-1] * pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
422  }
423  P *= (3.0/size);
424  estimator = P;
425  }
426 
427  return estimator;
428 }
429 
430 
432  if(!useCalibration)return;
433 
434 
435  TChain* t1 = new TChain("SiStripCalib/APVGain");
436  t1->Add(m_calibrationPath.c_str());
437 
438  unsigned int tree_DetId;
439  unsigned char tree_APVId;
440  double tree_Gain;
441 
442  t1->SetBranchAddress("DetId" ,&tree_DetId );
443  t1->SetBranchAddress("APVId" ,&tree_APVId );
444  t1->SetBranchAddress("Gain" ,&tree_Gain );
445 
446  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
447  t1->GetEntry(ientry);
448  stModInfo* MOD = MODsColl[tree_DetId];
449  MOD->Gain = tree_Gain;
450  }
451 
452  delete t1;
453 
454 }
455 
456 //define this as a plug-in
double path(double cosine, double thickness)
double p() const
momentum vector magnitude
Definition: TrackBase.h:129
T getParameter(std::string const &) const
virtual int nstrips() const =0
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
const SiStripCluster * GetCluster(const TrackerSingleRecHit *hit)
Definition: DeDxTools.h:27
virtual void produce(edm::Event &, const edm::EventSetup &) override
const_iterator end() const
last iterator over the map (read only)
LocalVector localDirection() const
double GetProbability(const SiStripCluster *cluster, TrajectoryStateOnSurface trajState, const uint32_t &)
fixed size error matrix
Definition: Error.h:29
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:53
DeDxDiscriminatorProducer(const edm::ParameterSet &)
#define P
#define NULL
Definition: scimark2.h:8
const Bounds & bounds() const
Definition: Surface.h:128
#define min(a, b)
Definition: mlp_lapack.h:161
virtual void beginRun(edm::Run const &run, const edm::EventSetup &) override
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
double charge(const std::vector< uint8_t > &Ampls)
uint32_t geographicalId() const
uint32_t rawId() const
get the raw id
Definition: DetId.h:45
LocalVector localMomentum() const
virtual float thickness() const =0
DataContainer const & measurements() const
Definition: Trajectory.h:215
int iEvent
Definition: GenABIO.cc:243
T mag() const
Definition: PV3DBase.h:67
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:41
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:141
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
T z() const
Definition: PV3DBase.h:64
int j
Definition: DBlmapReader.cc:9
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:39
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
int k[5][pyjets_maxn]
DEFINE_FWK_MODULE(CosmicTrackingParticleSelector)
Definition: DetId.h:20
bool shapeSelection(const std::vector< uint8_t > &ampls)
Definition: DeDxTools.cc:134
size_type size() const
map size
const T & get() const
Definition: EventSetup.h:55
T const * product() const
Definition: ESHandle.h:62
int ClusterSaturatingStrip(const SiStripCluster *cluster, const uint32_t &)
T const * product() const
Definition: Handle.h:74
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
unsigned short found() const
Number of valid hits on track.
Definition: Track.h:100
T transverse() const
Another name for perp()
double ComputeDiscriminator(std::vector< double > &vect_probs)
const_iterator begin() const
first iterator over the map (read only)
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:56
tuple size
Write out results.
const std::vector< uint8_t > & amplitudes() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
Definition: Run.h:36