CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SiPixelLorentzAngle.cc
Go to the documentation of this file.
1 
2 #include <memory>
3 #include <string>
4 #include <iostream>
5 #include <fstream>
6 #include <TMath.h>
7 #include "CalibTracker/SiPixelLorentzAngle/interface/SiPixelLorentzAngle.h"
8 
9 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
10 #include "FWCore/Framework/interface/ESHandle.h"
11 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
12 #include "DataFormats/GeometryVector/interface/LocalVector.h"
13 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
14 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
15 #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h"
16 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
17 #include "DataFormats/TrackReco/interface/TrackExtra.h"
18 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
19 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
20 #include "Geometry/CommonTopologies/interface/StripTopology.h"
21 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
22 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
23 #include "DataFormats/TrackReco/interface/TrackFwd.h"
24 
25 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
26 #include "Geometry/TrackerGeometryBuilder/interface/GluedGeomDet.h"
27 #include "TrackingTools/TrackFitters/interface/TrajectoryStateCombiner.h"
28 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
29 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
30 int lower_bin_;
31 
32 using namespace std;
33 using namespace edm;
34 using namespace reco;
35 
36 SiPixelLorentzAngle::SiPixelLorentzAngle(edm::ParameterSet const& conf) :
37  conf_(conf), filename_(conf.getParameter<std::string>("fileName")), filenameFit_(conf.getParameter<std::string>("fileNameFit")), ptmin_(conf.getParameter<double>("ptMin")), simData_(conf.getParameter<bool>("simData")), normChi2Max_(conf.getParameter<double>("normChi2Max")), clustSizeYMin_(conf.getParameter<int>("clustSizeYMin")), residualMax_(conf.getParameter<double>("residualMax")), clustChargeMax_(conf.getParameter<double>("clustChargeMax")),hist_depth_(conf.getParameter<int>("binsDepth")), hist_drift_(conf.getParameter<int>("binsDrift"))
38 {
39  // anglefinder_=new TrackLocalAngle(conf);
40  hist_x_ = 50;
41  hist_y_ = 100;
42  min_x_ = -500.;
43  max_x_ = 500.;
44  min_y_ = -1500.;
45  max_y_ = 500.;
46  width_ = 0.0285;
47  min_depth_ = -100.;
48  max_depth_ = 400.;
49  min_drift_ = -1000.; //-200.;(conf.getParameter<double>("residualMax"))
50  max_drift_ = 1000.; //400.;
51 
52 }
53 
54 // Virtual destructor needed.
55 SiPixelLorentzAngle::~SiPixelLorentzAngle() { }
56 
57 void SiPixelLorentzAngle::beginJob()
58 {
59 
60  // cout << "started SiPixelLorentzAngle" << endl;
61  hFile_ = new TFile (filename_.c_str(), "RECREATE" );
62  int bufsize = 64000;
63  // create tree structure
64  SiPixelLorentzAngleTree_ = new TTree("SiPixelLorentzAngleTree_","SiPixel LorentzAngle tree", bufsize);
65  SiPixelLorentzAngleTree_->Branch("run", &run_, "run/I", bufsize);
66  SiPixelLorentzAngleTree_->Branch("event", &event_, "event/I", bufsize);
67  SiPixelLorentzAngleTree_->Branch("module", &module_, "module/I", bufsize);
68  SiPixelLorentzAngleTree_->Branch("ladder", &ladder_, "ladder/I", bufsize);
69  SiPixelLorentzAngleTree_->Branch("layer", &layer_, "layer/I", bufsize);
70  SiPixelLorentzAngleTree_->Branch("isflipped", &isflipped_, "isflipped/I", bufsize);
71  SiPixelLorentzAngleTree_->Branch("pt", &pt_, "pt/F", bufsize);
72  SiPixelLorentzAngleTree_->Branch("eta", &eta_, "eta/F", bufsize);
73  SiPixelLorentzAngleTree_->Branch("phi", &phi_, "phi/F", bufsize);
74  SiPixelLorentzAngleTree_->Branch("chi2", &chi2_, "chi2/D", bufsize);
75  SiPixelLorentzAngleTree_->Branch("ndof", &ndof_, "ndof/D", bufsize);
76  SiPixelLorentzAngleTree_->Branch("trackhit", &trackhit_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
77  SiPixelLorentzAngleTree_->Branch("simhit", &simhit_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
78  SiPixelLorentzAngleTree_->Branch("npix", &pixinfo_.npix, "npix/I", bufsize);
79  SiPixelLorentzAngleTree_->Branch("rowpix", pixinfo_.row, "row[npix]/F", bufsize);
80  SiPixelLorentzAngleTree_->Branch("colpix", pixinfo_.col, "col[npix]/F", bufsize);
81  SiPixelLorentzAngleTree_->Branch("adc", pixinfo_.adc, "adc[npix]/F", bufsize);
82  SiPixelLorentzAngleTree_->Branch("xpix", pixinfo_.x, "x[npix]/F", bufsize);
83  SiPixelLorentzAngleTree_->Branch("ypix", pixinfo_.y, "y[npix]/F", bufsize);
84  SiPixelLorentzAngleTree_->Branch("clust", &clust_, "x/F:y/F:charge/F:size_x/I:size_y/I:maxPixelCol/I:maxPixelRow:minPixelCol/I:minPixelRow/I", bufsize);
85  SiPixelLorentzAngleTree_->Branch("rechit", &rechit_, "x/F:y/F", bufsize);
86 
87  SiPixelLorentzAngleTreeForward_ = new TTree("SiPixelLorentzAngleTreeForward_","SiPixel LorentzAngle tree forward", bufsize);
88  SiPixelLorentzAngleTreeForward_->Branch("run", &run_, "run/I", bufsize);
89  SiPixelLorentzAngleTreeForward_->Branch("event", &event_, "event/I", bufsize);
90  SiPixelLorentzAngleTreeForward_->Branch("side", &sideF_, "side/I", bufsize);
91  SiPixelLorentzAngleTreeForward_->Branch("disk", &diskF_, "disk/I", bufsize);
92  SiPixelLorentzAngleTreeForward_->Branch("blade", &bladeF_, "blade/I", bufsize);
93  SiPixelLorentzAngleTreeForward_->Branch("panel", &panelF_, "panel/I", bufsize);
94  SiPixelLorentzAngleTreeForward_->Branch("module", &moduleF_, "module/I", bufsize);
95  SiPixelLorentzAngleTreeForward_->Branch("pt", &pt_, "pt/F", bufsize);
96  SiPixelLorentzAngleTreeForward_->Branch("eta", &eta_, "eta/F", bufsize);
97  SiPixelLorentzAngleTreeForward_->Branch("phi", &phi_, "phi/F", bufsize);
98  SiPixelLorentzAngleTreeForward_->Branch("chi2", &chi2_, "chi2/D", bufsize);
99  SiPixelLorentzAngleTreeForward_->Branch("ndof", &ndof_, "ndof/D", bufsize);
100  SiPixelLorentzAngleTreeForward_->Branch("trackhit", &trackhitF_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
101  SiPixelLorentzAngleTreeForward_->Branch("simhit", &simhitF_, "x/F:y/F:alpha/D:beta/D:gamma_/D", bufsize);
102  SiPixelLorentzAngleTreeForward_->Branch("npix", &pixinfoF_.npix, "npix/I", bufsize);
103  SiPixelLorentzAngleTreeForward_->Branch("rowpix", pixinfoF_.row, "row[npix]/F", bufsize);
104  SiPixelLorentzAngleTreeForward_->Branch("colpix", pixinfoF_.col, "col[npix]/F", bufsize);
105  SiPixelLorentzAngleTreeForward_->Branch("adc", pixinfoF_.adc, "adc[npix]/F", bufsize);
106  SiPixelLorentzAngleTreeForward_->Branch("xpix", pixinfoF_.x, "x[npix]/F", bufsize);
107  SiPixelLorentzAngleTreeForward_->Branch("ypix", pixinfoF_.y, "y[npix]/F", bufsize);
108  SiPixelLorentzAngleTreeForward_->Branch("clust", &clustF_, "x/F:y/F:charge/F:size_x/I:size_y/I:maxPixelCol/I:maxPixelRow:minPixelCol/I:minPixelRow/I", bufsize);
109  SiPixelLorentzAngleTreeForward_->Branch("rechit", &rechitF_, "x/F:y/F", bufsize);
110 
111 
112  //book histograms
113  char name[128];
114  for(int i_module = 1; i_module<=8; i_module++){
115  for(int i_layer = 1; i_layer<=3; i_layer++){
116  sprintf(name, "h_drift_depth_adc_layer%i_module%i", i_layer, i_module);
117  _h_drift_depth_adc_[i_module + (i_layer -1) * 8] = new TH2F(name,name,hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
118  sprintf(name, "h_drift_depth_adc2_layer%i_module%i", i_layer, i_module);
119  _h_drift_depth_adc2_[i_module + (i_layer -1) * 8] = new TH2F(name,name,hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
120  sprintf(name, "h_drift_depth_noadc_layer%i_module%i", i_layer, i_module);
121  _h_drift_depth_noadc_[i_module + (i_layer -1) * 8] = new TH2F(name,name,hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
122  sprintf(name, "h_drift_depth_layer%i_module%i", i_layer, i_module);
123  _h_drift_depth_[i_module + (i_layer -1) * 8] = new TH2F(name,name,hist_drift_ , min_drift_, max_drift_, hist_depth_, min_depth_, max_depth_);
124  sprintf(name, "h_mean_layer%i_module%i", i_layer, i_module);
125  _h_mean_[i_module + (i_layer -1) * 8] = new TH1F(name,name,hist_depth_, min_depth_, max_depth_);
126  }
127  }
128 
129  // just for some expaining plots
130  h_cluster_shape_adc_ = new TH2F("h_cluster_shape_adc","cluster shape with adc weight", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
131  h_cluster_shape_noadc_ = new TH2F("h_cluster_shape_noadc","cluster shape without adc weight", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
132  h_cluster_shape_ = new TH2F("h_cluster_shape","cluster shape", hist_x_, min_x_, max_x_, hist_y_, min_y_, max_y_);
133  h_cluster_shape_adc_rot_ = new TH2F("h_cluster_shape_adc_rot","cluster shape with adc weight", hist_x_, min_x_, max_x_, hist_y_, -max_y_, -min_y_);
134  h_cluster_shape_noadc_rot_ = new TH2F("h_cluster_shape_noadc_rot","cluster shape without adc weight", hist_x_, min_x_, max_x_, hist_y_, -max_y_, -min_y_);
135  h_cluster_shape_rot_ = new TH2F("h_cluster_shape_rot","cluster shape", hist_x_, min_x_, max_x_, hist_y_, -max_y_, -min_y_);
136  h_tracks_ = new TH1F("h_tracks","h_tracks",2,0.,2.);
137  event_counter_ = 0;
138  trackEventsCounter_ = 0;
139  // trackcounter_ = 0;
140  hitCounter_ = 0;
141  usedHitCounter_ = 0;
142  pixelTracksCounter_ = 0;
143 // edm::ESHandle<TrackerGeometry> estracker; //this block should not be in beginJob()
144 // c.get<TrackerDigiGeometryRecord>().get(estracker);
145 // tracker=&(* estracker);
146 }
147 
148 
149 // Functions that gets called by framework every event
150 void SiPixelLorentzAngle::analyze(const edm::Event& e, const edm::EventSetup& es)
151 {
152  event_counter_++;
153  // if(event_counter_ % 500 == 0) cout << "event number " << event_counter_ << endl;
154  cout << "event number " << event_counter_ << endl;
155 
156  edm::ESHandle<TrackerGeometry> estracker;
157  es.get<TrackerDigiGeometryRecord>().get(estracker);
158  tracker=&(* estracker);
159 
160  TrackerHitAssociator* associate;
161  if(simData_) associate = new TrackerHitAssociator(e); else associate = 0;
162  // restet values
163  module_=-1;
164  layer_=-1;
165  ladder_ = -1;
166  isflipped_ = -1;
167  pt_ = -999;
168  eta_ = 999;
169  phi_ = 999;
170  pixinfo_.npix = 0;
171 
172  run_ = e.id().run();
173  event_ = e.id().event();
174 
175  // get the association map between tracks and trajectories
176  edm::Handle<TrajTrackAssociationCollection> trajTrackCollectionHandle;
177  e.getByLabel(conf_.getParameter<std::string>("src"),trajTrackCollectionHandle);
178  if(trajTrackCollectionHandle->size() >0){
179  trackEventsCounter_++;
180  for(TrajTrackAssociationCollection::const_iterator it = trajTrackCollectionHandle->begin(); it!=trajTrackCollectionHandle->end();++it){
181  const Track& track = *it->val;
182  const Trajectory& traj = *it->key;
183 
184  // get the trajectory measurements
185  std::vector<TrajectoryMeasurement> tmColl = traj.measurements();
186  // TrajectoryStateOnSurface tsos = tsoscomb( itTraj->forwardPredictedState(), itTraj->backwardPredictedState() );
187  pt_ = track.pt();
188  eta_ = track.eta();
189  phi_ = track.phi();
190  chi2_ = traj.chiSquared();
191  ndof_ = traj.ndof();
192  if(pt_ < ptmin_) continue;
193  // iterate over trajectory measurements
194  std::vector<PSimHit> matched;
195  h_tracks_->Fill(0);
196  bool pixeltrack = false;
197  for(std::vector<TrajectoryMeasurement>::const_iterator itTraj = tmColl.begin(); itTraj != tmColl.end(); itTraj++) {
198  if(! itTraj->updatedState().isValid()) continue;
199  TransientTrackingRecHit::ConstRecHitPointer recHit = itTraj->recHit();
200  if(! recHit->isValid() || recHit->geographicalId().det() != DetId::Tracker ) continue;
201  unsigned int subDetID = (recHit->geographicalId().subdetId());
202  if( subDetID == PixelSubdetector::PixelBarrel || subDetID == PixelSubdetector::PixelEndcap){
203  if(!pixeltrack){
204  h_tracks_->Fill(1);
205  pixelTracksCounter_++;
206  }
207  pixeltrack = true;
208  }
209  if( subDetID == PixelSubdetector::PixelBarrel){
210 
211  hitCounter_++;
212 
213  DetId detIdObj = recHit->geographicalId();
214  const PixelGeomDetUnit * theGeomDet = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(detIdObj) );
215  if(!theGeomDet) continue;
216 
217  const PixelTopology * topol = &(theGeomDet->specificTopology());
218 
219  if(!topol) continue;
220  PXBDetId pxbdetIdObj(detIdObj);
221  layer_ = pxbdetIdObj.layer();
222  ladder_ = pxbdetIdObj.ladder();
223  module_ = pxbdetIdObj.module();
224  float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
225  float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
226  if ( tmp2<tmp1 ) isflipped_ = 1;
227  else isflipped_ = 0;
228  const SiPixelRecHit * recHitPix = dynamic_cast<const SiPixelRecHit *>((*recHit).hit());
229  if(!recHitPix) continue;
230  rechit_.x = recHitPix->localPosition().x();
231  rechit_.y = recHitPix->localPosition().y();
232  SiPixelRecHit::ClusterRef const& cluster = recHitPix->cluster();
233 
234  // fill entries in clust_
235  clust_.x = (cluster)->x();
236  clust_.y = (cluster)->y();
237  clust_.charge = (cluster->charge())/1000.;
238  clust_.size_x = cluster->sizeX();
239  clust_.size_y = cluster->sizeY();
240  clust_.maxPixelCol = cluster->maxPixelCol();
241  clust_.maxPixelRow = cluster->maxPixelRow();
242  clust_.minPixelCol = cluster->minPixelCol();
243  clust_.minPixelRow = cluster->minPixelRow();
244  // fill entries in pixinfo_:
245  fillPix(*cluster ,topol, pixinfo_);
246  // fill the trackhit info
247  TrajectoryStateOnSurface tsos=itTraj->updatedState();
248  if(!tsos.isValid()){
249  cout << "tsos not valid" << endl;
250  continue;
251  }
252  LocalVector trackdirection=tsos.localDirection();
253  LocalPoint trackposition=tsos.localPosition();
254 
255  if(trackdirection.z()==0) continue;
256  // the local position and direction
257  trackhit_.alpha = atan2(trackdirection.z(),trackdirection.x());
258  trackhit_.beta = atan2(trackdirection.z(),trackdirection.y());
259  trackhit_.gamma = atan2(trackdirection.x(),trackdirection.y());
260  trackhit_.x = trackposition.x();
261  trackhit_.y = trackposition.y();
262 
263 
264  // fill entries in simhit_:
265  if(simData_){
266  matched.clear();
267  matched = associate->associateHit((*recHitPix));
268  float dr_start=9999.;
269  for (std::vector<PSimHit>::iterator isim = matched.begin(); isim != matched.end(); ++isim){
270  DetId simdetIdObj((*isim).detUnitId());
271  if (simdetIdObj == detIdObj) {
272  float sim_x1 = (*isim).entryPoint().x(); // width (row index, in col direction)
273  float sim_y1 = (*isim).entryPoint().y(); // length (col index, in row direction)
274  float sim_x2 = (*isim).exitPoint().x();
275  float sim_y2 = (*isim).exitPoint().y();
276  float sim_xpos = 0.5*(sim_x1+sim_x2);
277  float sim_ypos = 0.5*(sim_y1+sim_y2);
278  float sim_px = (*isim).momentumAtEntry().x();
279  float sim_py = (*isim).momentumAtEntry().y();
280  float sim_pz = (*isim).momentumAtEntry().z();
281 
282  float dr = (sim_xpos-(recHitPix->localPosition().x()))*(sim_xpos-recHitPix->localPosition().x()) +
283  (sim_ypos-recHitPix->localPosition().y())*(sim_ypos-recHitPix->localPosition().y());
284  if(dr<dr_start) {
285  simhit_.x = sim_xpos;
286  simhit_.y = sim_ypos;
287  simhit_.alpha = atan2(sim_pz, sim_px);
288  simhit_.beta = atan2(sim_pz, sim_py);
289  simhit_.gamma = atan2(sim_px, sim_py);
290  dr_start = dr;
291  }
292  }
293  } // end of filling simhit_
294  }
295  // is one pixel in cluster a large pixel ? (hit will be excluded)
296  bool large_pix = false;
297  for (int j = 0; j < pixinfo_.npix; j++){
298  int colpos = static_cast<int>(pixinfo_.col[j]);
299  if (pixinfo_.row[j] == 0 || pixinfo_.row[j] == 79 || pixinfo_.row[j] == 80 || pixinfo_.row[j] == 159 || colpos % 52 == 0 || colpos % 52 == 51 ){
300  large_pix = true;
301  }
302  }
303 
304  double residual = TMath::Sqrt( (trackhit_.x - rechit_.x) * (trackhit_.x - rechit_.x) + (trackhit_.y - rechit_.y) * (trackhit_.y - rechit_.y) );
305 
306  SiPixelLorentzAngleTree_->Fill();
307  if( !large_pix && (chi2_/ndof_) < normChi2Max_ && cluster->sizeY() >= clustSizeYMin_ && residual < residualMax_ && (cluster->charge() < clustChargeMax_)){
308  usedHitCounter_++;
309  // iterate over pixels in hit
310  for (int j = 0; j < pixinfo_.npix; j++){
311  // use trackhits
312  float dx = (pixinfo_.x[j] - (trackhit_.x - width_/2. / TMath::Tan(trackhit_.alpha))) * 10000.;
313  float dy = (pixinfo_.y[j] - (trackhit_.y - width_/2. / TMath::Tan(trackhit_.beta))) * 10000.;
314  float depth = dy * tan(trackhit_.beta);
315  float drift = dx - dy * tan(trackhit_.gamma);
316  _h_drift_depth_adc_[module_ + (layer_ -1) * 8]->Fill(drift, depth, pixinfo_.adc[j]);
317  _h_drift_depth_adc2_[module_ + (layer_ -1) * 8]->Fill(drift, depth, pixinfo_.adc[j]*pixinfo_.adc[j]);
318  _h_drift_depth_noadc_[module_ + (layer_ -1) * 8]->Fill(drift, depth);
319  if( layer_ == 3 && module_==1 && isflipped_){
320  float dx_rot = dx * TMath::Cos(trackhit_.gamma) + dy * TMath::Sin(trackhit_.gamma);
321  float dy_rot = dy * TMath::Cos(trackhit_.gamma) - dx * TMath::Sin(trackhit_.gamma) ;
322  h_cluster_shape_adc_->Fill(dx, dy, pixinfo_.adc[j]);
323  h_cluster_shape_noadc_->Fill(dx, dy);
324  h_cluster_shape_adc_rot_->Fill(dx_rot, dy_rot, pixinfo_.adc[j]);
325  h_cluster_shape_noadc_rot_->Fill(dx_rot, dy_rot);
326  }
327  }
328  }
329  } else if (subDetID == PixelSubdetector::PixelEndcap) {
330  DetId detIdObj = recHit->geographicalId();
331  const PixelGeomDetUnit * theGeomDet = dynamic_cast<const PixelGeomDetUnit*> ( tracker->idToDet(detIdObj) );
332  if(!theGeomDet) continue;
333 
334  const PixelTopology * topol = &(theGeomDet->specificTopology());
335 
336  if(!topol) continue;
337  PXFDetId pxfdetIdObj(detIdObj);
338  sideF_ = pxfdetIdObj.side();
339  diskF_ = pxfdetIdObj.disk();
340  bladeF_ = pxfdetIdObj.blade();
341  panelF_ = pxfdetIdObj.panel();
342  moduleF_ = pxfdetIdObj.module();
343  // float tmp1 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp();
344  // float tmp2 = theGeomDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp();
345  // if ( tmp2<tmp1 ) isflipped_ = 1;
346  // else isflipped_ = 0;
347  const SiPixelRecHit * recHitPix = dynamic_cast<const SiPixelRecHit *>((*recHit).hit());
348  if(!recHitPix) continue;
349  rechitF_.x = recHitPix->localPosition().x();
350  rechitF_.y = recHitPix->localPosition().y();
351  SiPixelRecHit::ClusterRef const& cluster = recHitPix->cluster();
352 
353  // fill entries in clust_
354  clustF_.x = (cluster)->x();
355  clustF_.y = (cluster)->y();
356  clustF_.charge = (cluster->charge())/1000.;
357  clustF_.size_x = cluster->sizeX();
358  clustF_.size_y = cluster->sizeY();
359  clustF_.maxPixelCol = cluster->maxPixelCol();
360  clustF_.maxPixelRow = cluster->maxPixelRow();
361  clustF_.minPixelCol = cluster->minPixelCol();
362  clustF_.minPixelRow = cluster->minPixelRow();
363  // fill entries in pixinfo_:
364  fillPix(*cluster ,topol, pixinfoF_);
365  // fill the trackhit info
366  TrajectoryStateOnSurface tsos=itTraj->updatedState();
367  if(!tsos.isValid()){
368  cout << "tsos not valid" << endl;
369  continue;
370  }
371  LocalVector trackdirection=tsos.localDirection();
372  LocalPoint trackposition=tsos.localPosition();
373 
374  if(trackdirection.z()==0) continue;
375  // the local position and direction
376  trackhitF_.alpha = atan2(trackdirection.z(),trackdirection.x());
377  trackhitF_.beta = atan2(trackdirection.z(),trackdirection.y());
378  trackhitF_.gamma = atan2(trackdirection.x(),trackdirection.y());
379  trackhitF_.x = trackposition.x();
380  trackhitF_.y = trackposition.y();
381 
382 
383  // fill entries in simhit_:
384  if(simData_){
385  matched.clear();
386  matched = associate->associateHit((*recHitPix));
387  float dr_start=9999.;
388  for (std::vector<PSimHit>::iterator isim = matched.begin(); isim != matched.end(); ++isim){
389  DetId simdetIdObj((*isim).detUnitId());
390  if (simdetIdObj == detIdObj) {
391  float sim_x1 = (*isim).entryPoint().x(); // width (row index, in col direction)
392  float sim_y1 = (*isim).entryPoint().y(); // length (col index, in row direction)
393  float sim_x2 = (*isim).exitPoint().x();
394  float sim_y2 = (*isim).exitPoint().y();
395  float sim_xpos = 0.5*(sim_x1+sim_x2);
396  float sim_ypos = 0.5*(sim_y1+sim_y2);
397  float sim_px = (*isim).momentumAtEntry().x();
398  float sim_py = (*isim).momentumAtEntry().y();
399  float sim_pz = (*isim).momentumAtEntry().z();
400 
401  float dr = (sim_xpos-(recHitPix->localPosition().x()))*(sim_xpos-recHitPix->localPosition().x()) +
402  (sim_ypos-recHitPix->localPosition().y())*(sim_ypos-recHitPix->localPosition().y());
403  if(dr<dr_start) {
404  simhitF_.x = sim_xpos;
405  simhitF_.y = sim_ypos;
406  simhitF_.alpha = atan2(sim_pz, sim_px);
407  simhitF_.beta = atan2(sim_pz, sim_py);
408  simhitF_.gamma = atan2(sim_px, sim_py);
409  dr_start = dr;
410  }
411  }
412  } // end of filling simhit_
413  }
414  SiPixelLorentzAngleTreeForward_->Fill();
415  }
416  } //end iteration over trajectory measurements
417  } //end iteration over trajectories
418  }
419 
420 
421 
422 
423 
424 }
425 
426 void SiPixelLorentzAngle::endJob()
427 {
428  // produce histograms with the average adc counts
429  for(int i_ring = 1; i_ring<=24; i_ring++){
430  _h_drift_depth_[i_ring]->Divide(_h_drift_depth_adc_[i_ring], _h_drift_depth_noadc_[i_ring]);
431  }
432 
433  h_drift_depth_adc_slice_ = new TH1F("h_drift_depth_adc_slice", "slice of adc histogram", hist_drift_ , min_drift_, max_drift_);
434 
435  TF1 *f1 = new TF1("f1","[0] + [1]*x",50., 235.);
436  f1->SetParName(0,"p0");
437  f1->SetParName(1,"p1");
438  f1->SetParameter(0,0);
439  f1->SetParameter(1,0.4);
440  ofstream fLorentzFit( filenameFit_.c_str(), ios::trunc );
441  cout.precision( 4 );
442  fLorentzFit << "module" << "\t" << "layer" << "\t" << "offset" << "\t" << "error" << "\t" << "slope" << "\t" << "error" << "\t" "rel.err" << "\t" "pull" << "\t" << "chi2" << "\t" << "prob" << endl;
443  //loop over modlues and layers to fit the lorentz angle
444  for( int i_layer = 1; i_layer<=3; i_layer++){
445  for(int i_module = 1; i_module<=8; i_module++){
446  //loop over bins in depth (z-local-coordinate) (in order to fit slices)
447  for( int i = 1; i <= hist_depth_; i++){
448  findMean(i, (i_module + (i_layer - 1) * 8));
449  }// end loop over bins in depth
450  _h_mean_[i_module + (i_layer - 1) * 8]->Fit(f1,"ERQ");
451  double p0 = f1->GetParameter(0);
452  double e0 = f1->GetParError(0);
453  double p1 = f1->GetParameter(1);
454  double e1 = f1->GetParError(1);
455  double chi2 = f1->GetChisquare();
456  double prob = f1->GetProb();
457  fLorentzFit << setprecision( 4 ) << i_module << "\t" << i_layer << "\t" << p0 << "\t" << e0 << "\t" << p1 << setprecision( 3 ) << "\t" << e1 << "\t" << e1 / p1 *100. << "\t"<< (p1 - 0.424) / e1 << "\t"<< chi2 << "\t" << prob << endl;
458  }
459  } // end loop over modules and layers
460  fLorentzFit.close();
461  hFile_->cd();
462  for(int i_module = 1; i_module<=8; i_module++){
463  for(int i_layer = 1; i_layer<=3; i_layer++){
464  _h_drift_depth_adc_[i_module + (i_layer -1) * 8]->Write();
465  _h_drift_depth_adc2_[i_module + (i_layer -1) * 8]->Write();
466  _h_drift_depth_noadc_[i_module + (i_layer -1) * 8]->Write();
467  _h_drift_depth_[i_module + (i_layer -1) * 8]->Write();
468  _h_mean_[i_module + (i_layer -1) * 8]->Write();
469  }
470  }
471  h_cluster_shape_adc_->Write();
472  h_cluster_shape_noadc_->Write();
473  h_cluster_shape_adc_rot_->Write();
474  h_cluster_shape_noadc_rot_->Write();
475  h_tracks_->Write();
476 
477  hFile_->Write();
478  hFile_->Close();
479  cout << "events: " << event_counter_ << endl;
480  cout << "events with tracks: " << trackEventsCounter_ << endl;
481  cout << "events with pixeltracks: " << pixelTracksCounter_ << endl;
482  cout << "hits in the pixel: " << hitCounter_ << endl;
483  cout << "number of used Hits: " << usedHitCounter_ << endl;
484 }
485 
486 inline void SiPixelLorentzAngle::fillPix(const SiPixelCluster & LocPix, const PixelTopology * topol, Pixinfo& pixinfo)
487 
488 {
489  const std::vector<SiPixelCluster::Pixel>& pixvector = LocPix.pixels();
490  pixinfo.npix = 0;
491  for(std::vector<SiPixelCluster::Pixel>::const_iterator itPix = pixvector.begin(); itPix != pixvector.end(); itPix++){
492  // for(pixinfo.npix = 0; pixinfo.npix < static_cast<int>(pixvector.size()); ++pixinfo.npix) {
493  pixinfo.row[pixinfo.npix] = itPix->x;
494  pixinfo.col[pixinfo.npix] = itPix->y;
495  pixinfo.adc[pixinfo.npix] = itPix->adc;
496  LocalPoint lp = topol->localPosition(MeasurementPoint(itPix->x + 0.5, itPix->y+0.5));
497  pixinfo.x[pixinfo.npix] = lp.x();
498  pixinfo.y[pixinfo.npix]= lp.y();
499  pixinfo.npix++;
500  }
501 }
502 
503 void SiPixelLorentzAngle::findMean(int i, int i_ring)
504 {
505  double nentries = 0;
506 
507  h_drift_depth_adc_slice_->Reset("ICE");
508 
509  // determine sigma and sigma^2 of the adc counts and average adc counts
510  //loop over bins in drift width
511  for( int j = 1; j<= hist_drift_; j++){
512  if(_h_drift_depth_noadc_[i_ring]->GetBinContent(j, i) >= 1){
513  double adc_error2 = (_h_drift_depth_adc2_[i_ring]->GetBinContent(j,i) - _h_drift_depth_adc_[i_ring]->GetBinContent(j,i)*_h_drift_depth_adc_[i_ring]->GetBinContent(j, i) / _h_drift_depth_noadc_[i_ring]->GetBinContent(j, i)) / _h_drift_depth_noadc_[i_ring]->GetBinContent(j, i);
514  _h_drift_depth_adc_[i_ring]->SetBinError(j, i, sqrt(adc_error2));
515  double error2 = adc_error2 / (_h_drift_depth_noadc_[i_ring]->GetBinContent(j,i) - 1.);
516  _h_drift_depth_[i_ring]->SetBinError(j,i,sqrt(error2));
517  }
518  else{
519  _h_drift_depth_[i_ring]->SetBinError(j,i,0);
520  _h_drift_depth_adc_[i_ring]->SetBinError(j, i, 0);
521  }
522  h_drift_depth_adc_slice_->SetBinContent(j, _h_drift_depth_adc_[i_ring]->GetBinContent(j,i));
523  h_drift_depth_adc_slice_->SetBinError(j, _h_drift_depth_adc_[i_ring]->GetBinError(j,i));
524  nentries += _h_drift_depth_noadc_[i_ring]->GetBinContent(j,i);
525  } // end loop over bins in drift width
526 
527  double mean = h_drift_depth_adc_slice_->GetMean(1);
528  double error = 0;
529  if(nentries != 0){
530  error = h_drift_depth_adc_slice_->GetRMS(1) / sqrt(nentries);
531  }
532 
533  _h_mean_[i_ring]->SetBinContent(i, mean);
534  _h_mean_[i_ring]->SetBinError(i, error);
535 
536 }
edm::EventID::run
RunNumber_t run() const
Definition: EventID.h:42
Surface::toGlobal
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:78
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
edm::EventID::event
EventNumber_t event() const
Definition: EventID.h:44
i
int i
Definition: DBlmapReader.cc:9
PXFDetId::panel
unsigned int panel() const
panel id
Definition: PXFDetId.h:52
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:71
Topology::localPosition
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
Pixinfo::adc
float adc[maxpix]
Definition: SiPixelLorentzAngle.h:54
TrajectoryStateOnSurface::localDirection
LocalVector localDirection() const
Definition: TrajectoryStateOnSurface.h:176
edm::AssociationMap::const_iterator
const iterator
Definition: AssociationMap.h:47
SiPixelLorentzAngle::fillPix
void fillPix(const SiPixelCluster &LocPix, const PixelTopology *topol, Pixinfo &pixinfo)
Definition: SiPixelLorentzAngle.cc:486
SiPixelLorentzAngle::_h_drift_depth_adc_
std::map< int, TH2F * > _h_drift_depth_adc_
Definition: SiPixelLorentzAngle.h:169
shallow::drift
LocalVector drift(const StripGeomDetUnit *, const MagneticField &, const SiStripLorentzAngle &)
Definition: ShallowTools.cc:39
PV3DBase::y
T y() const
Definition: PV3DBase.h:62
reco::TrackBase::phi
double phi() const
azimuthal angle of momentum vector
Definition: TrackBase.h:139
PXBDetId::ladder
unsigned int ladder() const
ladder id
Definition: PXBDetId.h:39
Pixinfo::row
float row[maxpix]
Definition: SiPixelLorentzAngle.h:52
dt_dqm_sourceclient_common_cff.reco
tuple reco
Definition: dt_dqm_sourceclient_common_cff.py:105
PXBDetId::layer
unsigned int layer() const
layer id
Definition: PXBDetId.h:35
estimatePileup_makeJSON.trunc
list trunc
Definition: estimatePileup_makeJSON.py:214
SiPixelLorentzAngle::_h_drift_depth_noadc_
std::map< int, TH2F * > _h_drift_depth_noadc_
Definition: SiPixelLorentzAngle.h:171
PXFDetId::blade
unsigned int blade() const
blade id
Definition: PXFDetId.h:48
SiPixelLorentzAngle::_h_drift_depth_
std::map< int, TH2F * > _h_drift_depth_
Definition: SiPixelLorentzAngle.h:172
SiPixelLorentzAngle::tracker
const TrackerGeometry * tracker
Definition: SiPixelLorentzAngle.h:194
Hit::beta
double beta
Definition: SiPixelLorentzAngle.h:63
Trajectory::measurements
DataContainer const & measurements() const
Definition: Trajectory.h:203
MeasurementPoint
Measurement2DPoint MeasurementPoint
Measurement points are two-dimensional by default.
Definition: MeasurementPoint.h:12
reco::TrackBase::eta
double eta() const
pseudorapidity of momentum vector
Definition: TrackBase.h:141
Fit
Definition: Fit.h:34
HcalObjRepresent::Fill
void Fill(HcalDetId &id, double val, std::vector< TH2F > &depth)
Definition: HcalObjRepresent.h:611
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
lower_bin_
int lower_bin_
Definition: SiPixelLorentzAngle.cc:30
reco::TrackBase::pt
double pt() const
track transverse momentum
Definition: TrackBase.h:131
PV3DBase::z
T z() const
Definition: PV3DBase.h:63
SiPixelLorentzAngle::conf_
edm::ParameterSet conf_
Definition: SiPixelLorentzAngle.h:144
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
PXBDetId::module
unsigned int module() const
det id
Definition: PXBDetId.h:43
j
int j
Definition: DBlmapReader.cc:9
PXFDetId::module
unsigned int module() const
det id
Definition: PXFDetId.h:56
TrackerGeometry::idToDet
virtual const GeomDet * idToDet(DetId) const
Definition: TrackerGeometry.cc:137
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
PXFDetId::disk
unsigned int disk() const
disk id
Definition: PXFDetId.h:43
SiPixelLorentzAngle::analyze
virtual void analyze(const edm::Event &e, const edm::EventSetup &c)
Definition: SiPixelLorentzAngle.cc:150
dbtoconf.conf
tuple conf
Definition: dbtoconf.py:185
Local3DPoint
Point3DBase< float, LocalTag > Local3DPoint
Definition: LocalPoint.h:9
DetId
Definition: DetId.h:20
Trajectory::ndof
int ndof(bool bon=true) const
Definition: Trajectory.cc:74
BaseTrackerRecHit::localPosition
virtual LocalPoint localPosition() const
Definition: BaseTrackerRecHit.cc:28
Hit::alpha
double alpha
Definition: SiPixelLorentzAngle.h:62
Pixinfo::x
float x[maxpix]
Definition: SiPixelLorentzAngle.h:55
SiPixelRecHit::cluster
ClusterRef cluster() const
Definition: SiPixelRecHit.h:41
edm::EventSetup::get
const T & get() const
Definition: EventSetup.h:55
PixelGeomDetUnit::specificTopology
virtual const PixelTopology & specificTopology() const
Returns a reference to the pixel proxy topology.
Definition: PixelGeomDetUnit.cc:19
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:56
SiPixelCluster
Pixel cluster – collection of neighboring pixels above threshold.
Definition: SiPixelCluster.h:23
p1
double p1[4]
Definition: TauolaWrapper.h:89
GeomDet::surface
const BoundPlane & surface() const
The nominal surface of the GeomDet.
Definition: GeomDet.h:35
TrackerHitAssociator::associateHit
std::vector< PSimHit > associateHit(const TrackingRecHit &thit)
Definition: TrackerHitAssociator.cc:120
PXFDetId::side
unsigned int side() const
positive or negative id
Definition: PXFDetId.h:38
Pixinfo::col
float col[maxpix]
Definition: SiPixelLorentzAngle.h:53
SiPixelLorentzAngle::_h_drift_depth_adc2_
std::map< int, TH2F * > _h_drift_depth_adc2_
Definition: SiPixelLorentzAngle.h:170
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
Hit::gamma
double gamma
Definition: SiPixelLorentzAngle.h:64
vdt::x
x
Definition: VDTMath.h:216
PV3DBase::x
T x() const
Definition: PV3DBase.h:61
SiPixelLorentzAngle::findMean
void findMean(int i, int i_ring)
Definition: SiPixelLorentzAngle.cc:503
SiPixelLorentzAngle::_h_mean_
std::map< int, TH1F * > _h_mean_
Definition: SiPixelLorentzAngle.h:174
SiPixelCluster::pixels
const std::vector< Pixel > pixels() const
Definition: SiPixelCluster.h:142
Trajectory::chiSquared
double chiSquared() const
Definition: Trajectory.h:242
SiPixelRecHit
Our base class.
Definition: SiPixelRecHit.h:22
Pixinfo::y
float y[maxpix]
Definition: SiPixelLorentzAngle.h:56
Generated for CMSSW Reference Manual by doxygen 1.8.5