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.17 2010/05/25 14:40:08 querten 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  useCalibration = iConfig.getParameter<bool>("UseCalibration");
75  m_calibrationPath = iConfig.getParameter<string>("calibrationPath");
76 
77  Prob_ChargePath = NULL;
78 }
79 
80 
82 
83 // ------------ method called once each job just before starting event loop ------------
85 {
87  if( strcmp(Reccord.c_str(),"SiStripDeDxMip_3D_Rcd")==0){
88  iSetup.get<SiStripDeDxMip_3D_Rcd>() .get(DeDxMapHandle_);
89  }else if(strcmp(Reccord.c_str(),"SiStripDeDxPion_3D_Rcd")==0){
90  iSetup.get<SiStripDeDxPion_3D_Rcd>().get(DeDxMapHandle_);
91  }else if(strcmp(Reccord.c_str(),"SiStripDeDxKaon_3D_Rcd")==0){
92  iSetup.get<SiStripDeDxKaon_3D_Rcd>().get(DeDxMapHandle_);
93  }else if(strcmp(Reccord.c_str(),"SiStripDeDxProton_3D_Rcd")==0){
94  iSetup.get<SiStripDeDxProton_3D_Rcd>().get(DeDxMapHandle_);
95  }else if(strcmp(Reccord.c_str(),"SiStripDeDxElectron_3D_Rcd")==0){
96  iSetup.get<SiStripDeDxElectron_3D_Rcd>().get(DeDxMapHandle_);
97  }else{
98 // printf("The reccord %s is not known by the DeDxDiscriminatorProducer\n", Reccord.c_str());
99 // printf("Program will exit now\n");
100  exit(0);
101  }
102  DeDxMap_ = *DeDxMapHandle_.product();
103 
104  double xmin = DeDxMap_.rangeX().min;
105  double xmax = DeDxMap_.rangeX().max;
106  double ymin = DeDxMap_.rangeY().min;
107  double ymax = DeDxMap_.rangeY().max;
108  double zmin = DeDxMap_.rangeZ().min;
109  double zmax = DeDxMap_.rangeZ().max;
110 
111  if(Prob_ChargePath)delete Prob_ChargePath;
112  Prob_ChargePath = new TH3D ("Prob_ChargePath" , "Prob_ChargePath" , DeDxMap_.numberOfBinsX(), xmin, xmax, DeDxMap_.numberOfBinsY() , ymin, ymax, DeDxMap_.numberOfBinsZ(), zmin, zmax);
113 
114 
115 
116  if(strcmp(ProbabilityMode.c_str(),"Accumulation")==0){
117 // printf("LOOOP ON P\n");
118  for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
119 // printf("LOOOP ON PATH\n");
120  for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
121 // printf("LOOOP ON CHARGE\n");
122 
123  double Ni = 0;
124  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);}
125 
126  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
127  double tmp = 0;
128  for(int l=0;l<=k;l++){ tmp+=DeDxMap_.binContent(i,j,l);}
129 
130  if(Ni>0){
131  Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
132 // 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);
133  }else{
134  Prob_ChargePath->SetBinContent (i, j, k, 0);
135  }
136  }
137  }
138  }
139  }else if(strcmp(ProbabilityMode.c_str(),"Integral")==0){
140 // printf("LOOOP ON P\n");
141  for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
142 // printf("LOOOP ON PATH\n");
143  for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
144 // printf("LOOOP ON CHARGE\n");
145 
146  double Ni = 0;
147  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);}
148 
149  for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
150  double tmp = DeDxMap_.binContent(i,j,k);
151 
152  if(Ni>0){
153  Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
154 // 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);
155  }else{
156  Prob_ChargePath->SetBinContent (i, j, k, 0);
157  }
158  }
159  }
160  }
161  }else{
162 // printf("The ProbabilityMode: %s is unknown\n",ProbabilityMode.c_str());
163 // printf("The program will stop now\n");
164  exit(0);
165  }
166 
167 
168 
169 /*
170  for(int i=0;i<Prob_ChargePath->GetXaxis()->GetNbins();i++){
171  for(int j=0;j<Prob_ChargePath->GetYaxis()->GetNbins();j++){
172  double tmp = DeDxMap_.binContent(i,j);
173  Prob_ChargePath->SetBinContent (i, j, tmp);
174  printf("%i %i --> %f\n",i,j,tmp);
175  }
176  }
177 */
178 
179 
180 
181  if(MODsColl.size()==0){
183  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
184  m_tracker = tkGeom.product();
185 
186  vector<GeomDet*> Det = tkGeom->dets();
187  for(unsigned int i=0;i<Det.size();i++){
188  DetId Detid = Det[i]->geographicalId();
189  int SubDet = Detid.subdetId();
190 
191  if( SubDet == StripSubdetector::TIB || SubDet == StripSubdetector::TID ||
192  SubDet == StripSubdetector::TOB || SubDet == StripSubdetector::TEC ){
193 
194  StripGeomDetUnit* DetUnit = dynamic_cast<StripGeomDetUnit*> (Det[i]);
195  if(!DetUnit)continue;
196 
197  const StripTopology& Topo = DetUnit->specificTopology();
198  unsigned int NAPV = Topo.nstrips()/128;
199 
200  double Eta = DetUnit->position().basicVector().eta();
201  double R = DetUnit->position().basicVector().transverse();
202  double Thick = DetUnit->surface().bounds().thickness();
203 
204  stModInfo* MOD = new stModInfo;
205  MOD->DetId = Detid.rawId();
206  MOD->SubDet = SubDet;
207  MOD->Eta = Eta;
208  MOD->R = R;
209  MOD->Thickness = Thick;
210  MOD->NAPV = NAPV;
211  MODsColl[MOD->DetId] = MOD;
212  }
213  }
214 
215  MakeCalibrationMap();
216  }
217 
218 }
219 
221 {
222  MODsColl.clear();
223 
224 /*
225  TFile* file = new TFile("MipsMap.root", "RECREATE");
226  Prob_ChargePath->Write();
227  file->Write();
228  file->Close();
229 */
230 }
231 
232 
233 
235 {
236  auto_ptr<ValueMap<DeDxData> > trackDeDxDiscrimAssociation(new ValueMap<DeDxData> );
237  ValueMap<DeDxData>::Filler filler(*trackDeDxDiscrimAssociation);
238 
239  Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
240  iEvent.getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
241  const TrajTrackAssociationCollection TrajToTrackMap = *trajTrackAssociationHandle.product();
242 
243  edm::Handle<reco::TrackCollection> trackCollectionHandle;
244  iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
245 
247  iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
248  m_tracker = tkGeom.product();
249 
250  std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.size() );
251 
252  unsigned track_index = 0;
253  for(TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap.begin(); it!=TrajToTrackMap.end(); ++it, track_index++) {
254  dEdxDiscrims[track_index] = DeDxData(-1, -2, 0 );
255 
256  const Track track = *it->val;
257  const Trajectory traj = *it->key;
258 
259  if(track.eta() <MinTrackEta || track.eta()>MaxTrackEta ){continue;}
260  if(track.p() <MinTrackMomentum || track.p() >MaxTrackMomentum){continue;}
261  if(track.found()<MinTrackHits ){continue;}
262 
263  std::vector<double> vect_probs;
264  vector<TrajectoryMeasurement> measurements = traj.measurements();
265 
266  unsigned int NClusterSaturating = 0;
267  for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
268 
269  TrajectoryStateOnSurface trajState = measurement_it->updatedState();
270  if( !trajState.isValid() ) continue;
271 
272  const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
273  const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
274  const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
275  const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
276 
277  double Prob;
278  if(sistripsimplehit)
279  {
280  Prob = GetProbability((sistripsimplehit->cluster()).get(), trajState); if(Prob>=0) vect_probs.push_back(Prob);
281  if(ClusterSaturatingStrip((sistripsimplehit->cluster()).get())>0)NClusterSaturating++;
282  }else if(sistripmatchedhit){
283  Prob = GetProbability((sistripmatchedhit->monoHit()->cluster()).get(), trajState); if(Prob>=0) vect_probs.push_back(Prob);
284  Prob = GetProbability((sistripmatchedhit->stereoHit()->cluster()).get(), trajState); if(Prob>=0) vect_probs.push_back(Prob);
285  if(ClusterSaturatingStrip((sistripmatchedhit->monoHit()->cluster()).get()) >0)NClusterSaturating++;
286  if(ClusterSaturatingStrip((sistripmatchedhit->stereoHit()->cluster()).get())>0)NClusterSaturating++;
287  }else if(sistripsimple1dhit){
288  Prob = GetProbability((sistripsimple1dhit->cluster()).get(), trajState); if(Prob>=0) vect_probs.push_back(Prob);
289  if(ClusterSaturatingStrip((sistripsimple1dhit->cluster()).get())>0)NClusterSaturating++;
290  }else{
291  }
292  }
293 
294  double estimator = ComputeDiscriminator(vect_probs);
295  int size = vect_probs.size();
296  float Error = -1;
297 
298  //WARNING: Since the dEdX Error is not properly computed for the moment
299  //It was decided to store the number of saturating cluster in that dataformat
300  Error = NClusterSaturating;
301  dEdxDiscrims[track_index] = DeDxData(estimator, Error, size );
302 
303 // printf("%i --> %g\n",size,estimator);
304  }
305 
306  filler.insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
307  filler.fill();
308  iEvent.put(trackDeDxDiscrimAssociation);
309 }
310 
311 
313  const vector<uint8_t>& ampls = cluster->amplitudes();
314 
315  stModInfo* MOD = NULL;
316  if(useCalibration)MOD = MODsColl[cluster->geographicalId()];
317 
318  int SaturatingStrip = 0;
319  for(unsigned int i=0;i<ampls.size();i++){
320  int StripCharge = ampls[i];
321  if(MOD){StripCharge = (int)(StripCharge / MOD->Gain);}
322  if(StripCharge>=254)SaturatingStrip++;
323  }
324  return SaturatingStrip;
325 }
326 
328 {
329  // Get All needed variables
330  LocalVector trackDirection = trajState.localDirection();
331  double cosine = trackDirection.z()/trackDirection.mag();
332  const vector<uint8_t>& ampls = cluster->amplitudes();
333  uint32_t detId = cluster->geographicalId();
334 // int firstStrip = cluster->firstStrip();
335  stModInfo* MOD = MODsColl[detId];
336 
337 
338  // Sanity Checks
339  if( ampls.size()>MaxNrStrips) {return -1;}
340 // if( DeDxDiscriminatorTools::IsSaturatingStrip (ampls)) {return -1;}
341 // if( DeDxDiscriminatorTools::IsSpanningOver2APV (firstStrip, ampls.size())) {return -1;}
342 // if(!DeDxDiscriminatorTools::IsFarFromBorder (trajState, m_tracker->idToDetUnit(DetId(detId)) )) {return -1;}
343 
344 
345  // Find Probability for this given Charge and Path
346  double charge = 0;
347  if(useCalibration){
348  for(unsigned int i=0;i<ampls.size();i++){
349  int CalibratedCharge = ampls[i];
350  CalibratedCharge = (int)(CalibratedCharge / MOD->Gain);
351  if(CalibratedCharge>=1024){
352  CalibratedCharge = 255;
353  }else if(CalibratedCharge>254){
354  CalibratedCharge = 254;
355  }
356  charge+=CalibratedCharge;
357  }
358  }else{
359  charge = DeDxDiscriminatorTools::charge(ampls);
360  }
361  double path = DeDxDiscriminatorTools::path(cosine,MOD->Thickness);
362 
363  int BinX = Prob_ChargePath->GetXaxis()->FindBin(trajState.localMomentum().mag());
364  int BinY = Prob_ChargePath->GetYaxis()->FindBin(path);
365  int BinZ = Prob_ChargePath->GetZaxis()->FindBin(charge/path);
366 
367  return Prob_ChargePath->GetBinContent(BinX,BinY,BinZ);
368 }
369 
370 
371 double DeDxDiscriminatorProducer::ComputeDiscriminator(std::vector<double>& vect_probs)
372 {
373  double estimator = -1;
374  int size = vect_probs.size();
375  if(size<=0)return estimator;
376 
377  if(Formula==0){
378  double P = 1;
379  for(int i=0;i<size;i++){
380  if(vect_probs[i]<=0.0001){P *= pow(0.0001 , 1.0/size);}
381  else {P *= pow(vect_probs[i], 1.0/size);}
382  }
383  estimator = P;
384  }else if(Formula==1){
385  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
386  for(int i=0;i<size;i++){if(vect_probs[i]<=0.0001)vect_probs[i] = 0.0001; }
387 
388  double SumJet = 0.;
389  for(int i=0;i<size;i++){ SumJet+= log(vect_probs[i]); }
390 
391  double Loginvlog=log(-SumJet);
392  double Prob =1.;
393  double lfact=1.;
394 
395  for(int l=1; l!=size; l++){
396  lfact*=l;
397  Prob+=exp(l*Loginvlog-log(1.*lfact));
398  }
399 
400  double LogProb=log(Prob);
401  double ProbJet=std::min(exp(std::max(LogProb+SumJet,-30.)),1.);
402  estimator = -log10(ProbJet)/4.;
403  estimator = 1-estimator;
404  }else if(Formula==2){
405  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
406  double P = 1.0/(12*size);
407  for(int i=1;i<=size;i++){
408  P += pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
409  }
410  P *= (3.0/size);
411  estimator = P;
412  }else{
413  std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
414  double P = 1.0/(12*size);
415  for(int i=1;i<=size;i++){
416  P += vect_probs[i-1] * pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
417  }
418  P *= (3.0/size);
419  estimator = P;
420  }
421 
422  return estimator;
423 }
424 
425 
427  if(!useCalibration)return;
428 
429 
430  TChain* t1 = new TChain("SiStripCalib/APVGain");
431  t1->Add(m_calibrationPath.c_str());
432 
433  unsigned int tree_DetId;
434  unsigned char tree_APVId;
435  double tree_Gain;
436 
437  t1->SetBranchAddress("DetId" ,&tree_DetId );
438  t1->SetBranchAddress("APVId" ,&tree_APVId );
439  t1->SetBranchAddress("Gain" ,&tree_Gain );
440 
441  for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
442  t1->GetEntry(ientry);
443  stModInfo* MOD = MODsColl[tree_DetId];
444  MOD->Gain = tree_Gain;
445  }
446 
447  delete t1;
448 
449 }
450 
451 //define this as a plug-in
double path(double cosine, double thickness)
double p() const
momentum vector magnitude
Definition: TrackBase.h:128
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
double GetProbability(const SiStripCluster *cluster, TrajectoryStateOnSurface trajState)
const_iterator end() const
last iterator over the map (read only)
const SiStripRecHit2D * stereoHit() const
LocalVector localDirection() const
int ClusterSaturatingStrip(const SiStripCluster *cluster)
fixed size error matrix
Definition: Error.h:29
void insert(const H &h, I begin, I end)
Definition: ValueMap.h:52
DeDxDiscriminatorProducer(const edm::ParameterSet &)
#define P
#define NULL
Definition: scimark2.h:8
#define min(a, b)
Definition: mlp_lapack.h:161
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
virtual const StripTopology & specificTopology() const
Returns a reference to the strip proxy topology.
double charge(const std::vector< uint8_t > &Ampls)
int path() const
Definition: HLTadd.h:3
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:169
int iEvent
Definition: GenABIO.cc:243
T mag() const
Definition: PV3DBase.h:61
const Surface::PositionType & position() const
The position (origin of the R.F.)
Definition: GeomDet.h:43
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:140
const T & max(const T &a, const T &b)
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:84
virtual void beginRun(edm::Run &run, const edm::EventSetup &)
T z() const
Definition: PV3DBase.h:58
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:359
int k[5][pyjets_maxn]
DEFINE_FWK_MODULE(CosmicTrackingParticleSelector)
Definition: DetId.h:20
const Bounds & bounds() const
Definition: BoundSurface.h:89
Log< T >::type log(const T &t)
Definition: Log.h:22
size_type size() const
map size
virtual void produce(edm::Event &, const edm::EventSetup &)
const T & get() const
Definition: EventSetup.h:55
ClusterRef const & cluster() const
T const * product() const
Definition: ESHandle.h:62
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)
ClusterRef const & cluster() const
const_iterator begin() const
first iterator over the map (read only)
const SiStripRecHit2D * monoHit() const
const BasicVectorType & basicVector() const
Definition: PV3DBase.h:54
virtual const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:37
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:31