CMS 3D CMS Logo

SiStripLAProfileBooker.cc
Go to the documentation of this file.
1 /* VI Janurary 2012
2  * This file need to be migrated to the new interface of matched hit as the mono/stero are not there anymore
3  * what is returned are hits w/o localpoistion, just cluster and id
4  */
5 #include <string>
6 #include <iostream>
7 #include <fstream>
8 #include <functional>
9 
10 
14 
20 
25 
35 
40 
43 
44 #include <TProfile.h>
45 #include <TStyle.h>
46 #include <TF1.h>
47 
48 #include<list>
49 
50 class DetIdLess
51  : public std::binary_function< const SiStripRecHit2D*,const SiStripRecHit2D*,bool> {
52 public:
53 
54  DetIdLess() {}
55 
56  bool operator()( const SiStripRecHit2D* a, const SiStripRecHit2D* b) const {
57  return *(a->cluster())<*(b->cluster());
58  }
59 };
60 
61 
62  //Constructor
63 
65  conf_(conf)
66 {
67 }
68 
69  //BeginRun
70 
72 
73  //Retrieve tracker topology from geometry
75  c.get<TrackerTopologyRcd>().get(tTopoHandle);
76  const TrackerTopology* const tTopo = tTopoHandle.product();
77 
78  //get magnetic field and geometry from ES
80  c.get<IdealMagneticFieldRecord>().get(esmagfield);
81  const MagneticField * magfield=&(*esmagfield);
82 
84  c.get<TrackerDigiGeometryRecord>().get(estracker);
85  tracker=&(*estracker);
86 
87  std::vector<uint32_t> activeDets;
88  edm::ESHandle<SiStripDetCabling> tkmechstruct=nullptr;
89  if (conf_.getParameter<bool>("UseStripCablingDB")){
90  c.get<SiStripDetCablingRcd>().get(tkmechstruct);
91  activeDets.clear();
92  tkmechstruct->addActiveDetectorsRawIds(activeDets);
93  }
94  else {
95  const TrackerGeometry::DetIdContainer& Id = estracker->detIds();
96  TrackerGeometry::DetIdContainer::const_iterator Iditer;
97  activeDets.clear();
98  for(Iditer=Id.begin();Iditer!=Id.end();Iditer++){
99  activeDets.push_back(Iditer->rawId());
100  }
101  }
102 
103  edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
104  //Get Ids;
105  double ModuleRangeMin=conf_.getParameter<double>("ModuleXMin");
106  double ModuleRangeMax=conf_.getParameter<double>("ModuleXMax");
107  double TIBRangeMin=conf_.getParameter<double>("TIBXMin");
108  double TIBRangeMax=conf_.getParameter<double>("TIBXMax");
109  double TOBRangeMin=conf_.getParameter<double>("TOBXMin");
110  double TOBRangeMax=conf_.getParameter<double>("TOBXMax");
111  int TIB_bin=conf_.getParameter<int>("TIB_bin");
112  int TOB_bin=conf_.getParameter<int>("TOB_bin");
113  int SUM_bin=conf_.getParameter<int>("SUM_bin");
114 
115  hFile = new TFile (conf_.getUntrackedParameter<std::string>("treeName").c_str(), "RECREATE" );
116 
117  Hit_Tree = hFile->mkdir("Hit_Tree");
118  Track_Tree = hFile->mkdir("Track_Tree");
119  Event_Tree = hFile->mkdir("Event_Tree");
120 
121  HitsTree = new TTree("HitsTree", "HitsTree");
122 
123  HitsTree->Branch("RunNumber", &RunNumber, "RunNumber/I");
124  HitsTree->Branch("EventNumber", &EventNumber, "EventNumber/I");
125  HitsTree->Branch("TanTrackAngle", &TanTrackAngle, "TanTrackAngle/F");
126  HitsTree->Branch("TanTrackAngleParallel", &TanTrackAngleParallel, "TanTrackAngleParallel/F");
127  HitsTree->Branch("ClSize", &ClSize, "ClSize/I");
128  HitsTree->Branch("HitCharge", &HitCharge, "HitCharge/I");
129  HitsTree->Branch("Hit_Std_Dev", &hit_std_dev, "hit_std_dev/F");
130  HitsTree->Branch("Type", &Type, "Type/I");
131  HitsTree->Branch("Layer", &Layer, "Layer/I");
132  HitsTree->Branch("Wheel", &Wheel, "Wheel/I");
133  HitsTree->Branch("bw_fw", &bw_fw, "bw_fw/I");
134  HitsTree->Branch("Ext_Int", &Ext_Int, "Ext_Int/I");
135  HitsTree->Branch("MonoStereo", &MonoStereo, "MonoStereo/I");
136  HitsTree->Branch("MagField", &MagField, "MagField/F");
137  HitsTree->Branch("SignCorrection", &SignCorrection, "SignCorrection/F");
138  HitsTree->Branch("XGlobal", &XGlobal, "XGlobal/F");
139  HitsTree->Branch("YGlobal", &YGlobal, "YGlobal/F");
140  HitsTree->Branch("ZGlobal", &ZGlobal, "ZGlobal/F");
141  HitsTree->Branch("ParticleCharge", &ParticleCharge, "ParticleCharge/I");
142  HitsTree->Branch("Momentum", &Momentum, "Momentum/F");
143  HitsTree->Branch("pt", &pt, "pt/F");
144  HitsTree->Branch("chi2norm", &chi2norm, "chi2norm/F");
145  HitsTree->Branch("EtaTrack", &EtaTrack, "EtaTrack/F");
146  HitsTree->Branch("PhiTrack", &PhiTrack, "PhiTrack/F");
147  HitsTree->Branch("TrajSize", &trajsize, "trajsize/I");
148  HitsTree->Branch("HitNr", &HitNr, "HitNr/I");
149  HitsTree->Branch("HitPerTrack", &HitPerTrack, "HitPerTrack/I");
150  HitsTree->Branch("id_detector", &id_detector, "id_detector/I");
151  HitsTree->Branch("thick_detector", &thick_detector, "thick_detector/F");
152  HitsTree->Branch("pitch_detector", &pitch_detector, "pitch_detector/F");
153  HitsTree->Branch("Amplitudes", Amplitudes, "Amplitudes[ClSize]/I");
154 
155  HitsTree->SetDirectory(Hit_Tree);
156 
157  TrackTree = new TTree("TrackTree", "TrackTree");
158 
159  TrackTree->Branch("TrackCounter", &TrackCounter, "TrackCounter/I");
160 
161  TrackTree->SetDirectory(Track_Tree);
162 
163  EventTree = new TTree("EventTree", "EventTree");
164 
165  EventTree->Branch("EventCounter", &EventCounter, "EventCounter/I");
166 
167  EventTree->SetDirectory(Event_Tree);
168 
169 
170  // use SistripHistoId for producing histogram id (and title)
171  SiStripHistoId hidmanager;
172 
173  // create SiStripFolderOrganizer
174  SiStripFolderOrganizer folder_organizer;
175 
177 
178  //get all detids
179 
180  for(std::vector<uint32_t>::const_iterator Id = activeDets.begin(); Id!=activeDets.end(); Id++){
181 
182  // for(Iditer=Id.begin();Iditer!=Id.end();Iditer++){ //loop on detids
183  DetId Iditero=DetId(*Id);
184  DetId *Iditer=&Iditero;
185  if((Iditer->subdetId() == int(StripSubdetector::TIB)) || (Iditer->subdetId() == int(StripSubdetector::TOB))){ //include only barrel
186 
187  int module_bin = 0;
188  if(Iditer->subdetId() == int(StripSubdetector::TIB)){
189  module_bin = TIB_bin;
190  }else{
191  module_bin = TOB_bin;
192  }
193 
194  // create a TProfile for each module
195  StripSubdetector subid(*Iditer);
196  std::string hid;
197  //Mono single sided detectors
198  LocalPoint p;
199  auto stripdet = tracker->idToDet(subid);
200  if(!stripdet->isLeaf())continue;
201  const StripTopology& topol=(const StripTopology&)stripdet->topology();
202  float thickness=stripdet->specificSurface().bounds().thickness();
203 
204  folder_organizer.setDetectorFolder(Iditer->rawId(), tTopo);
205  hid = hidmanager.createHistoId(TkTag.label(),"det",Iditer->rawId());
206  MonitorElement * profile=dbe_->bookProfile(hid,hid,module_bin,ModuleRangeMin,ModuleRangeMax,20,0,5,"");
207  detparameters *param=new detparameters;
208  histos[Iditer->rawId()] = profile;
209  detmap[Iditer->rawId()] = param;
210  param->thickness = thickness*10000;
211  param->pitch = topol.localPitch(p)*10000;
212 
213  const GlobalPoint globalp = (stripdet->surface()).toGlobal(p);
214  GlobalVector globalmagdir = magfield->inTesla(globalp);
215  param->magfield=(stripdet->surface()).toLocal(globalmagdir);
216 
217  profile->setAxisTitle("tan(#theta_{t})",1);
218  profile->setAxisTitle("Cluster size",2);
219 
220  // create a summary histo if it does not exist already
222  unsigned int layerid;
223  getlayer(subid,tTopo,name,layerid);
224  name+=TkTag.label();
225  if(summaryhisto.find(layerid)==(summaryhisto.end())){
226  folder_organizer.setSiStripFolder();
227  MonitorElement * summaryprofile=nullptr;
229  summaryprofile=dbe_->bookProfile(name,name,SUM_bin,TIBRangeMin,TIBRangeMax,20,0,5,"");
230  else if (subid.subdetId()==int (StripSubdetector::TOB)||subid.subdetId()==int (StripSubdetector::TEC))
231  summaryprofile=dbe_->bookProfile(name,name,SUM_bin,TOBRangeMin,TOBRangeMax,20,0,5,"");
232  if(summaryprofile){
233  detparameters *summaryparam=new detparameters;
234  summaryhisto[layerid] = summaryprofile;
235  summarydetmap[layerid] = summaryparam;
236  summaryparam->thickness = thickness*10000;
237  summaryparam->pitch = topol.localPitch(p)*10000;
238  summaryprofile->setAxisTitle("tan(#theta_{t})",1);
239  summaryprofile->setAxisTitle("Cluster size",2);
240  }
241  }
242 
243  }
244  }
245 
246  trackcollsize = 0;
247  trajsize = 0;
248  RunNumber = 0;
249  EventNumber = 0;
250  hitcounter = 0;
251  hitcounter_2ndloop = 0;
252  worse_double_hit = 0;
253  better_double_hit = 0;
254  eventcounter = 0;
255 
256  EventCounter = 1;
257  TrackCounter = 1;
258 
259 }
260 
262  detparmap::iterator detpariter;
263  for( detpariter=detmap.begin(); detpariter!=detmap.end();++detpariter)delete detpariter->second;
264  for( detpariter=summarydetmap.begin(); detpariter!=summarydetmap.end();++detpariter)delete detpariter->second;
265  delete hFile;
266 }
267 
268 // Analyzer: Functions that gets called by framework every event
269 
271 {
272  //Retrieve tracker topology from geometry
273  edm::ESHandle<TrackerTopology> tTopoHandle;
274  es.get<TrackerTopologyRcd>().get(tTopoHandle);
275  const TrackerTopology* const tTopo = tTopoHandle.product();
276 
277  RunNumber = e.id().run();
278  EventNumber = e.id().event();
279 
280  eventcounter++;
281 
282  EventTree->Fill();
283 
284  //Analysis of Trajectory-RecHits
285 
286  edm::InputTag TkTag = conf_.getParameter<edm::InputTag>("Tracks");
287 
289  e.getByLabel(TkTag,trackCollection);
290 
292  e.getByLabel(TkTag,TrajectoryCollection);
293 
295  e.getByLabel(TkTag, TrajTrackMap);
296 
297  const reco::TrackCollection *tracks=trackCollection.product();
298 
299  // FIXME this has to be changed to use pointers to clusters...
300  std::map<const SiStripRecHit2D*,std::pair<float,float>,DetIdLess> hitangleassociation;
301  std::list<SiStripRecHit2D> cache; // ugly, inefficient, effective in making the above working
302 
303  trackcollsize=tracks->size();
304  trajsize=TrajectoryCollection->size();
305 
306  edm::LogInfo("SiStripLAProfileBooker::analyze") <<" Number of tracks in event = "<<trackcollsize<<"\n";
307  edm::LogInfo("SiStripLAProfileBooker::analyze") <<" Number of trajectories in event = "<<trajsize<<"\n";
308 
310 
311  for(TrajTrackIter = TrajTrackMap->begin(); TrajTrackIter!= TrajTrackMap->end(); TrajTrackIter++){ //loop on trajectories
312 
313  if(TrajTrackIter->key->foundHits()>=5){
314 
315  TrackTree->Fill();
316 
317  ParticleCharge = -99;
318  Momentum = -99;
319  pt = -99;
320  chi2norm = -99;
321  HitPerTrack = -99;
322  EtaTrack = -99;
323  PhiTrack = -99;
324 
325  ParticleCharge = TrajTrackIter->val->charge();
326  pt = TrajTrackIter->val->pt();
327  Momentum = TrajTrackIter->val->p();
328  chi2norm = TrajTrackIter->val->normalizedChi2();
329  EtaTrack = TrajTrackIter->val->eta();
330  PhiTrack = TrajTrackIter->val->phi();
331  HitPerTrack = TrajTrackIter->key->foundHits();
332 
333  std::vector<TrajectoryMeasurement> TMeas=TrajTrackIter->key->measurements();
334  std::vector<TrajectoryMeasurement>::iterator itm;
335 
336  for (itm=TMeas.begin();itm!=TMeas.end();itm++){ //loop on hits
337 
338  int i;
339  for(i=0;i<100;i++){Amplitudes[i]=0;}
340 
341  TanTrackAngle = -99;
343  ClSize = -99;
344  HitCharge = 0;
345  Type = -99;
346  Layer = -99;
347  Wheel = -99;
348  bw_fw = -99;
349  Ext_Int = -99;
350  MonoStereo = -99;
351  MagField = -99;
352  SignCorrection = -99;
353  XGlobal = -99;
354  YGlobal = -99;
355  ZGlobal = -99;
356  barycenter = -99;
357  hit_std_dev = -99;
358  sumx = 0;
359  id_detector=-1;
360  thick_detector=-1;
361  pitch_detector=-1;
362  HitNr = 1;
363 
364  SiStripRecHit2D lhit;
365  TrajectoryStateOnSurface tsos=itm->updatedState();
366  const TransientTrackingRecHit::ConstRecHitPointer thit=itm->recHit();
367  if((thit->geographicalId().subdetId() == int(StripSubdetector::TIB)) || thit->geographicalId().subdetId()== int(StripSubdetector::TOB)){ //include only barrel
368  const SiStripMatchedRecHit2D* matchedhit=dynamic_cast<const SiStripMatchedRecHit2D*>((*thit).hit());
369  const ProjectedSiStripRecHit2D* phit=dynamic_cast<const ProjectedSiStripRecHit2D*>((*thit).hit());
370  const SiStripRecHit2D* hit=dynamic_cast<const SiStripRecHit2D*>((*thit).hit());
371  if(phit) {lhit = phit->originalHit(); hit = &lhit;}
372 
373  LocalVector trackdirection=tsos.localDirection();
374 
375  if(matchedhit){//if matched hit...
376 
377  GluedGeomDet * gdet=(GluedGeomDet *)tracker->idToDet(matchedhit->geographicalId());
378 
379  GlobalVector gtrkdir=gdet->toGlobal(trackdirection);
380 
381  // THIS THE POINTER TO THE MONO HIT OF A MATCHED HIT
382 
383  // top be migrated to the more direct interface of matchedhit
384  cache.push_back(matchedhit->monoHit());
385  const SiStripRecHit2D * monohit = &cache.back();
386  const SiStripRecHit2D::ClusterRef & monocluster=monohit->cluster();
387  const GeomDetUnit * monodet=gdet->monoDet();
388  // this does not exists anymore! either project the matched or use CPE
389  const LocalPoint monoposition = monohit->localPosition();
390 
392  id_detector = detid.rawId();
394  const StripTopology& mtopol=(StripTopology&)monodet->topology();
395  pitch_detector = mtopol.localPitch(monoposition);
396  const GlobalPoint monogposition = (monodet->surface()).toGlobal(monoposition);
397  ClSize = (monocluster->amplitudes()).size();
398 
399  const auto & amplitudes = monocluster->amplitudes();
400 
401  barycenter = monocluster->barycenter()- 0.5;
402  uint16_t FirstStrip = monocluster->firstStrip();
403  auto begin=amplitudes.begin();
404  nstrip=0;
405  for(auto idigi=begin; idigi!=amplitudes.end(); idigi++){
406  Amplitudes[nstrip]=*idigi;
407  sumx+=pow(((FirstStrip+idigi-begin)-barycenter),2)*(*idigi);
408  HitCharge+=*idigi;
409  }
411 
412 
413  XGlobal = monogposition.x();
414  YGlobal = monogposition.y();
415  ZGlobal = monogposition.z();
416 
417  Type = detid.subdetId();
418  MonoStereo=detid.stereo();
419 
420  if(detid.subdetId() == int (StripSubdetector::TIB)){
421 
422  Layer = tTopo->tibLayer(detid);
423  bw_fw = tTopo->tibStringInfo(detid)[0];
424  Ext_Int = tTopo->tibStringInfo(detid)[1];
425  }
426  if(detid.subdetId() == int (StripSubdetector::TOB)){
427 
428  Layer = tTopo->tobLayer(detid);
429  bw_fw = tTopo->tobRodInfo(detid)[0];
430  }
431  if(detid.subdetId() == int (StripSubdetector::TID)){
432 
433  Wheel = tTopo->tidWheel(detid);
434  bw_fw = tTopo->tidModuleInfo(detid)[0];
435  }
436  if(detid.subdetId() == int (StripSubdetector::TEC)){
437 
438  Wheel = tTopo->tecWheel(detid);
439  bw_fw = tTopo->tecPetalInfo(detid)[0];
440  }
441 
442 
443  LocalVector monotkdir=monodet->toLocal(gtrkdir);
444 
445  if(monotkdir.z()!=0){
446 
447  // THE LOCAL ANGLE (MONO)
448  float tanangle = monotkdir.x()/monotkdir.z();
449  TanTrackAngleParallel = monotkdir.y()/monotkdir.z();
450  TanTrackAngle = tanangle;
451  detparmap::iterator TheDet=detmap.find(detid.rawId());
452  LocalVector localmagdir;
453  if(TheDet!=detmap.end())localmagdir=TheDet->second->magfield;
454  MagField = localmagdir.mag();
455  if(MagField != 0.){
456  LocalVector monoylocal(0,1,0);
457  float signcorrection = (localmagdir * monoylocal)/(MagField);
458  if(signcorrection!=0)SignCorrection=1/signcorrection;
459  }
460 
461  std::map<const SiStripRecHit2D *,std::pair<float,float>,DetIdLess>::iterator alreadystored=hitangleassociation.find(monohit);
462 
463 
464  if(alreadystored != hitangleassociation.end()){//decide which hit take
465  if(itm->estimate() > alreadystored->second.first){
466  worse_double_hit++;}
467  if(itm->estimate() < alreadystored->second.first){
469  hitangleassociation.insert(std::make_pair(monohit, std::make_pair(itm->estimate(),tanangle)));
470 
471  }}
472  else{
473  hitangleassociation.insert(make_pair(monohit, std::make_pair(itm->estimate(),tanangle)));
474  HitsTree->Fill();
475  hitcounter++;}
476 
477  // THIS THE POINTER TO THE STEREO HIT OF A MATCHED HIT
478 
479  // top be migrated to the more direct interface of matchedhit
480  cache.push_back(matchedhit->stereoHit());
481  const SiStripRecHit2D * stereohit = &cache.back();
482  const SiStripRecHit2D::ClusterRef & stereocluster=stereohit->cluster();
483  const GeomDetUnit * stereodet=gdet->stereoDet();
484  // this does not exists anymore! either project the matched or use CPE
485  const LocalPoint stereoposition = stereohit->localPosition();
487  id_detector = detid.rawId();
489  const StripTopology& stopol=(StripTopology&)stereodet->topology();
490  pitch_detector = stopol.localPitch(stereoposition);
491  const GlobalPoint stereogposition = (stereodet->surface()).toGlobal(stereoposition);
492 
493  ClSize = (stereocluster->amplitudes()).size();
494 
495  const auto & amplitudes = stereocluster->amplitudes();
496 
497  barycenter = stereocluster->barycenter()- 0.5;
498  uint16_t FirstStrip = stereocluster->firstStrip();
499  auto begin=amplitudes.begin();
500  nstrip=0;
501  for(auto idigi=begin; idigi!=amplitudes.end(); idigi++){
502  Amplitudes[nstrip]=*idigi;
503  sumx+=pow(((FirstStrip+idigi-begin)-barycenter),2)*(*idigi);
504  HitCharge+=*idigi;
505  }
507 
508  XGlobal = stereogposition.x();
509  YGlobal = stereogposition.y();
510  ZGlobal = stereogposition.z();
511 
512  Type = detid.subdetId();
513  MonoStereo=detid.stereo();
514 
515  if(detid.subdetId() == int (StripSubdetector::TIB)){
516 
517  Layer = tTopo->tibLayer(detid);
518  bw_fw = tTopo->tibStringInfo(detid)[0];
519  Ext_Int = tTopo->tibStringInfo(detid)[1];
520  }
521  if(detid.subdetId() == int (StripSubdetector::TOB)){
522 
523  Layer = tTopo->tobLayer(detid);
524  bw_fw = tTopo->tobRodInfo(detid)[0];
525  }
526  if(detid.subdetId() == int (StripSubdetector::TID)){
527 
528  Wheel = tTopo->tidWheel(detid);
529  bw_fw = tTopo->tidModuleInfo(detid)[0];
530  }
531  if(detid.subdetId() == int (StripSubdetector::TEC)){
532 
533  Wheel = tTopo->tecWheel(detid);
534  bw_fw = tTopo->tecPetalInfo(detid)[0];
535  }
536 
537 
538  LocalVector stereotkdir=stereodet->toLocal(gtrkdir);
539 
540  if(stereotkdir.z()!=0){
541 
542  // THE LOCAL ANGLE (STEREO)
543  float tanangle = stereotkdir.x()/stereotkdir.z();
544  TanTrackAngleParallel = stereotkdir.y()/stereotkdir.z();
545  TanTrackAngle = tanangle;
546  detparmap::iterator TheDet=detmap.find(detid.rawId());
547  LocalVector localmagdir;
548  if(TheDet!=detmap.end())localmagdir=TheDet->second->magfield;
549  MagField = localmagdir.mag();
550  LocalVector stereoylocal(0,1,0);
551  if(MagField != 0.){
552  float signcorrection = (localmagdir * stereoylocal)/(MagField);
553  if(signcorrection!=0)SignCorrection=1/signcorrection;}
554 
555  std::map<const SiStripRecHit2D *,std::pair<float,float>,DetIdLess>::iterator alreadystored=hitangleassociation.find(stereohit);
556 
557  if(alreadystored != hitangleassociation.end()){//decide which hit take
558  if(itm->estimate() > alreadystored->second.first){
559  worse_double_hit++;}
560  if(itm->estimate() < alreadystored->second.first){
562  hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(),tanangle)));
563 
564  }}
565  else{
566  hitangleassociation.insert(std::make_pair(stereohit, std::make_pair(itm->estimate(),tanangle)));
567  HitsTree->Fill();
568  hitcounter++;}
569 
570  }
571  }
572  }
573  else if(hit){
574 
575 
576  // hit= POINTER TO THE RECHIT
577 
578  const SiStripRecHit2D::ClusterRef & cluster=hit->cluster();
579 
581  const LocalPoint position = hit->localPosition();
583  id_detector = detid.rawId();
585  const StripTopology& topol=(StripTopology&)gdet->topology();
586  pitch_detector = topol.localPitch(position);
587  const GlobalPoint gposition = (gdet->surface()).toGlobal(position);
588 
589  ClSize = (cluster->amplitudes()).size();
590 
591  const auto & amplitudes = cluster->amplitudes();
592 
593  barycenter = cluster->barycenter()- 0.5;
594  uint16_t FirstStrip = cluster->firstStrip();
595  nstrip=0;
596  auto begin =amplitudes.begin();
597  for(auto idigi=amplitudes.begin(); idigi!=amplitudes.end(); idigi++){
598  Amplitudes[nstrip]=*idigi;
599  sumx+=pow(((FirstStrip+idigi-begin)-barycenter),2)*(*idigi);
600  HitCharge+=*idigi;
601  }
603 
604  XGlobal = gposition.x();
605  YGlobal = gposition.y();
606  ZGlobal = gposition.z();
607 
608  Type = detid.subdetId();
609  MonoStereo=detid.stereo();
610 
611  if(detid.subdetId() == int (StripSubdetector::TIB)){
612 
613  Layer = tTopo->tibLayer(detid);
614  bw_fw = tTopo->tibStringInfo(detid)[0];
615  Ext_Int = tTopo->tibStringInfo(detid)[1];
616  }
617  if(detid.subdetId() == int (StripSubdetector::TOB)){
618 
619  Layer = tTopo->tobLayer(detid);
620  bw_fw = tTopo->tobRodInfo(detid)[0];
621  }
622  if(detid.subdetId() == int (StripSubdetector::TID)){
623 
624  Wheel = tTopo->tidWheel(detid);
625  bw_fw = tTopo->tidModuleInfo(detid)[0];
626  }
627  if(detid.subdetId() == int (StripSubdetector::TEC)){
628 
629  Wheel = tTopo->tecWheel(detid);
630  bw_fw = tTopo->tecPetalInfo(detid)[0];
631  }
632 
633  if(trackdirection.z()!=0){
634 
635  // THE LOCAL ANGLE
636  float tanangle = trackdirection.x()/trackdirection.z();
637  TanTrackAngleParallel = trackdirection.y()/trackdirection.z();
638  TanTrackAngle = tanangle;
639  detparmap::iterator TheDet=detmap.find(detid.rawId());
640  LocalVector localmagdir;
641  if(TheDet!=detmap.end())localmagdir=TheDet->second->magfield;
642  MagField = localmagdir.mag();
643  if(MagField != 0.){
644  LocalVector ylocal(0,1,0);
645  float signcorrection = (localmagdir * ylocal)/(MagField);
646  if(signcorrection!=0)SignCorrection=1/signcorrection;}
647 
648  std::map<const SiStripRecHit2D *,std::pair<float,float>, DetIdLess>::iterator alreadystored=hitangleassociation.find(hit);
649 
650  if(alreadystored != hitangleassociation.end()){//decide which hit take
651  if(itm->estimate() > alreadystored->second.first){
652  worse_double_hit++;}
653  if(itm->estimate() < alreadystored->second.first){
655  hitangleassociation.insert(std::make_pair(hit, std::make_pair(itm->estimate(),tanangle)));
656 
657  }}
658  else{
659  hitangleassociation.insert(std::make_pair(hit,std::make_pair(itm->estimate(), tanangle) ) );
660  HitsTree->Fill();
661  hitcounter++;}
662 
663 
664  }
665  }
666  }
667  }
668  }
669  }
670  std::map<const SiStripRecHit2D *,std::pair<float,float>,DetIdLess>::iterator hitsiter;
671 
672 
673  for(hitsiter=hitangleassociation.begin();hitsiter!=hitangleassociation.end();hitsiter++){
674 
676 
677  const SiStripRecHit2D* hit=hitsiter->first;
678  const SiStripRecHit2D::ClusterRef & cluster=hit->cluster();
679 
680  size=(cluster->amplitudes()).size();
681 
683 
684  float tangent = hitsiter->second.second;
685 
686  //Sign and XZ plane projection correction applied in TrackLocalAngle (TIB|TOB layers)
687 
688  detparmap::iterator thedet=detmap.find(detid.rawId());
689  LocalVector localmagdir;
690  if(thedet!=detmap.end())localmagdir=thedet->second->magfield;
691  float localmagfield = localmagdir.mag();
692 
693  if(localmagfield != 0.){
694 
695  LocalVector ylocal(0,1,0);
696 
697  float normprojection = (localmagdir * ylocal)/(localmagfield);
698 
699  if(normprojection == 0.)LogDebug("SiStripLAProfileBooker::analyze")<<"Error: YBprojection = 0";
700 
701  else{
702  float signprojcorrection = 1/normprojection;
703  tangent*=signprojcorrection;
704  }
705  }
706 
707  //Filling histograms
708 
709  histomap::iterator thehisto=histos.find(detid.rawId());
710 
711  if(thehisto==histos.end())edm::LogError("SiStripLAProfileBooker::analyze")<<"Error: the profile associated to"<<detid.rawId()<<"does not exist! ";
712  else thehisto->second->Fill(tangent,size);
713 
714  //Summary histograms
716  unsigned int layerid;
717  getlayer(detid,tTopo,name,layerid);
718  histomap::iterator thesummaryhisto=summaryhisto.find(layerid);
719  if(thesummaryhisto==summaryhisto.end())edm::LogError("SiStripLAProfileBooker::analyze")<<"Error: the profile associated to subdet "<<name<<"does not exist! ";
720  else thesummaryhisto->second->Fill(tangent,size);
721 
722  }
723 
724 
725 }
726 
727 
728 //Makename function
729 
730 void SiStripLAProfileBooker::getlayer(const DetId & detid, const TrackerTopology* tTopo, std::string &name,unsigned int &layerid){
731  int layer=0;
732  std::stringstream layernum;
733 
734  if(detid.subdetId() == int (StripSubdetector::TIB)){
735 
736  name+="TIB_Layer_";
737  layer = tTopo->tibLayer(detid);
738  }
739 
740  else if(detid.subdetId() == int (StripSubdetector::TID)){
741 
742  name+="TID_Ring_";
743  layer = tTopo->tidRing(detid);
744  }
745 
746  else if(detid.subdetId() == int (StripSubdetector::TOB)){
747 
748  name+="TOB_Layer_";
749  layer = tTopo->tobLayer(detid);
750 
751  }
752 
753  else if(detid.subdetId() == int (StripSubdetector::TEC)){
754 
755  name+="TEC_Ring_";
756  layer = tTopo->tecRing(detid);
757  }
758  layernum<<layer;
759  name+=layernum.str();
760  layerid=detid.subdetId()*10+layer;
761 
762 }
763 
765 
766  std::string outputFile_ =conf_.getUntrackedParameter<std::string>("fileName", "LorentzAngle.root");
767  dbe_->save(outputFile_);
768 
769  hFile->Write();
770  hFile->Close();
771 }
#define LogDebug(id)
RunNumber_t run() const
Definition: EventID.h:39
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:41
T getUntrackedParameter(std::string const &, T const &) const
void addActiveDetectorsRawIds(std::vector< uint32_t > &) const
unsigned int tibLayer(const DetId &id) const
unsigned int tidRing(const DetId &id) const
const GeomDetUnit * monoDet() const
Definition: GluedGeomDet.h:20
const_iterator end() const
last iterator over the map (read only)
virtual const Topology & topology() const
Definition: GeomDet.cc:81
std::vector< unsigned int > tidModuleInfo(const DetId &id) const
LocalVector localDirection() const
GlobalPoint toGlobal(const Local2DPoint &lp) const
Conversion to the global R.F. from the R.F. of the GeomDet.
Definition: GeomDet.h:54
unsigned int tecRing(const DetId &id) const
ring id
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
T y() const
Definition: PV3DBase.h:63
LocalPoint toLocal(const GlobalPoint &gp) const
Conversion to the R.F. of the GeomDet.
Definition: GeomDet.h:69
const Bounds & bounds() const
Definition: Surface.h:120
unsigned int tidWheel(const DetId &id) const
LocalVector toLocal(const reco::Track::Vector &v, const Surface &s)
std::vector< unsigned int > tibStringInfo(const DetId &id) const
const Plane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:42
void setDetectorFolder(uint32_t rawdetid, const TrackerTopology *tTopo)
uint32_t rawId() const
get the raw id
Definition: DetId.h:43
std::vector< unsigned int > tecPetalInfo(const DetId &id) const
T mag() const
Definition: PV3DBase.h:67
std::vector< unsigned int > tobRodInfo(const DetId &id) const
std::shared_ptr< TrackingRecHit const > ConstRecHitPointer
void getlayer(const DetId &detid, const TrackerTopology *tTopo, std::string &name, unsigned int &layerid)
T sqrt(T t)
Definition: SSEVec.h:18
T z() const
Definition: PV3DBase.h:64
void analyze(const edm::Event &e, const edm::EventSetup &c) override
const DetIdContainer & detIds() const override
Returm a vector of all GeomDet DetIds (including those of GeomDetUnits)
void save(const std::string &filename, const std::string &path="", const std::string &pattern="", const std::string &rewrite="", const uint32_t run=0, const uint32_t lumi=0, SaveReferenceTag ref=SaveWithReference, int minStatus=dqm::qstatus::STATUS_OK, const std::string &fileupdate="RECREATE")
Definition: DQMStore.cc:2702
ClusterRef cluster() const
MonitorElement * bookProfile(const char *name, const char *title, int nchX, double lowX, double highX, int nchY, double lowY, double highY, const char *option="s")
Definition: DQMStore.cc:1372
int subdetId() const
get the contents of the subdetector field (not cast into any detector&#39;s numbering enum) ...
Definition: DetId.h:37
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:475
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
def cache(function)
SiStripRecHit2D originalHit() const
Definition: DetId.h:18
unsigned int stereo() const
stereo
SiStripRecHit2D stereoHit() const
virtual float thickness() const =0
std::string createHistoId(std::string description, std::string id_type, uint32_t component_id)
const T & get() const
Definition: EventSetup.h:58
bool operator()(const SiStripRecHit2D *a, const SiStripRecHit2D *b) const
double b
Definition: hdecay.h:120
std::vector< DetId > DetIdContainer
virtual float localPitch(const LocalPoint &) const =0
std::string const & label() const
Definition: InputTag.h:36
std::vector< Trajectory > TrajectoryCollection
SiStripRecHit2D monoHit() const
void beginRun(edm::Run const &, const edm::EventSetup &c) override
edm::EventID id() const
Definition: EventBase.h:60
#define begin
Definition: vmac.h:32
double a
Definition: hdecay.h:121
static int position[264][3]
Definition: ReadPGInfo.cc:509
const TrackerGeomDet * idToDet(DetId) const override
LocalPoint localPosition() const final
const_iterator begin() const
first iterator over the map (read only)
DetId geographicalId() const
SiStripLAProfileBooker(const edm::ParameterSet &conf)
T x() const
Definition: PV3DBase.h:62
unsigned int tecWheel(const DetId &id) const
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
T const * product() const
Definition: ESHandle.h:86
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
const GeomDetUnit * stereoDet() const
Definition: GluedGeomDet.h:21
const Plane & specificSurface() const
Same as surface(), kept for backward compatibility.
Definition: GeomDet.h:45
unsigned int tobLayer(const DetId &id) const
Definition: Run.h:43
const TrackerGeometry * tracker