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