CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
SiStripElectronAnalyzer.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: RecoEgamma/Examples
4 // Class : SiStripElectronAnalyzer
5 //
6 // Implementation:
7 // <Notes on implementation>
8 //
9 // Original Author: Jim Pivarski
10 // Created: Fri May 26 16:49:38 EDT 2006
11 //
12 
13 #include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h"
14 #include "DataFormats/Common/interface/Handle.h"
15 #include "DataFormats/DetId/interface/DetId.h"
16 #include "DataFormats/EgammaCandidates/interface/Electron.h"
17 #include "DataFormats/EgammaCandidates/interface/ElectronFwd.h"
18 #include "DataFormats/EgammaCandidates/interface/SiStripElectron.h"
19 #include "DataFormats/EgammaCandidates/interface/SiStripElectronFwd.h"
20 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
21 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
22 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
23 #include "DataFormats/Math/interface/Point3D.h"
24 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
25 #include "DataFormats/TrackReco/interface/Track.h"
26 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
27 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
28 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
29 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
30 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
31 #include "FWCore/Framework/interface/ESHandle.h"
32 #include "FWCore/Framework/interface/Event.h"
33 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
34 #include "FWCore/MessageLogger/interface/MessageLogger.h"
35 #include "FWCore/ParameterSet/interface/ParameterSet.h"
36 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
37 #include "Geometry/CommonDetUnit/interface/GeomDetEnumerators.h"
38 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
39 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
40 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
41 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
42 #include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"
43 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
44 
45 #include "TBranch.h"
46 #include "TFile.h"
47 #include "TH1F.h"
48 #include "TNtuple.h"
49 #include "TTree.h"
50 
51 #include <memory>
52 #include <map>
53 #include <cmath>
54 
55 #define myMaxHits 1000
56 
57 class SiStripElectronAnalyzer : public edm::one::EDAnalyzer<> {
58 public:
59  explicit SiStripElectronAnalyzer(const edm::ParameterSet&);
60  ~SiStripElectronAnalyzer() override;
61 
62  void analyze(const edm::Event&, const edm::EventSetup&) override;
63  void beginJob() override;
64  virtual void initNtuple(void);
65  void endJob(void) override;
66 
67 private:
68  double unwrapPhi(double phi) const {
69  while (phi > M_PI) {
70  phi -= 2. * M_PI;
71  }
72  while (phi < -M_PI) {
73  phi += 2. * M_PI;
74  }
75  return phi;
76  }
77 
78  // ----------member data ---------------------------
79  std::string fileName_;
80 
81  TFile* file_;
82  TH1F* numCand_;
83  TH1F* numElectrons_;
84  TH1F* numSuperClusters_;
85  TH1F* energySuperClusters_;
86  TH1F* sizeSuperClusters_;
87  TH1F* emaxSuperClusters_;
88  TH1F* phiWidthSuperClusters_;
89 
90  TH1F* energySuperClustersEl_;
91  TH1F* sizeSuperClustersEl_;
92  TH1F* emaxSuperClustersEl_;
93  TH1F* phiWidthSuperClustersEl_;
94 
95  TH1F* ptDiff;
96  TH1F* pDiff;
97  TH1F* pElectronFailed;
98  TH1F* ptElectronFailed;
99  TH1F* pElectronPassed;
100  TH1F* ptElectronPassed;
101  TH1F* sizeSuperClustersPassed;
102  TH1F* sizeSuperClustersFailed;
103  TH1F* energySuperClustersPassed;
104  TH1F* energySuperClustersFailed;
105  TH1F* eOverPFailed;
106  TH1F* eOverPPassed;
107 
108  TH1F* numSiStereoHits_;
109  TH1F* numSiMonoHits_;
110  TH1F* numSiMatchedHits_;
111 
112  TTree* myTree_;
113 
114  int NShowers_;
115  float EShower_[myMaxHits];
116  float XShower_[myMaxHits];
117  float YShower_[myMaxHits];
118  float ZShower_[myMaxHits];
119 
120  int NStereoHits_;
121  float StereoHitX_[myMaxHits];
122  float StereoHitY_[myMaxHits];
123  float StereoHitZ_[myMaxHits];
124 
125  float StereoHitR_[myMaxHits];
126  float StereoHitPhi_[myMaxHits];
127  float StereoHitTheta_[myMaxHits];
128 
129  // errors in local coords
130  float StereoHitSigX_[myMaxHits];
131  float StereoHitSigY_[myMaxHits];
132  float StereoHitCorr_[myMaxHits];
133 
134  float StereoHitSignal_[myMaxHits];
135  float StereoHitNoise_[myMaxHits];
136  int StereoHitWidth_[myMaxHits];
137 
138  int StereoDetector_[myMaxHits];
139  int StereoLayer_[myMaxHits];
140 
141  // mono corresponds to "rphi" only hits
142  int NMonoHits_;
143  float MonoHitX_[myMaxHits];
144  float MonoHitY_[myMaxHits];
145  float MonoHitZ_[myMaxHits];
146 
147  float MonoHitR_[myMaxHits];
148  float MonoHitPhi_[myMaxHits];
149  float MonoHitTheta_[myMaxHits];
150 
151  // errors in local coords
152  float MonoHitSigX_[myMaxHits];
153  float MonoHitSigY_[myMaxHits];
154  float MonoHitCorr_[myMaxHits];
155 
156  float MonoHitSignal_[myMaxHits];
157  float MonoHitNoise_[myMaxHits];
158  int MonoHitWidth_[myMaxHits];
159 
160  int MonoDetector_[myMaxHits];
161  int MonoLayer_[myMaxHits];
162 
163  // matched hits
164  int NMatchedHits_;
165  float MatchedHitX_[myMaxHits];
166  float MatchedHitY_[myMaxHits];
167  float MatchedHitZ_[myMaxHits];
168 
169  float MatchedHitR_[myMaxHits];
170  float MatchedHitPhi_[myMaxHits];
171  float MatchedHitTheta_[myMaxHits];
172 
173  // errors in local coords
174  float MatchedHitSigX_[myMaxHits];
175  float MatchedHitSigY_[myMaxHits];
176  float MatchedHitCorr_[myMaxHits];
177 
178  float MatchedHitSignal_[myMaxHits];
179  float MatchedHitNoise_[myMaxHits];
180  int MatchedHitWidth_[myMaxHits];
181 
182  int MatchedDetector_[myMaxHits];
183  int MatchedLayer_[myMaxHits];
184 
185  std::string mctruthProducer_;
186  std::string mctruthCollection_;
187  std::string superClusterProducer_;
188  std::string superClusterCollection_;
189  std::string basicClusterProducer_;
190  std::string basicClusterCollection_;
191  std::string eBRecHitProducer_;
192  std::string eBRecHitCollection_;
193  std::string siElectronProducer_;
194  std::string siElectronCollection_;
195  std::string electronProducer_;
196  std::string electronCollection_;
197  std::string siHitProducer_;
198  std::string siRphiHitCollection_;
199  std::string siStereoHitCollection_;
200  std::string siMatchedHitCollection_;
201 };
202 
203 #include "FWCore/Framework/interface/MakerMacros.h"
204 DEFINE_FWK_MODULE(SiStripElectronAnalyzer);
205 
206 SiStripElectronAnalyzer::SiStripElectronAnalyzer(const edm::ParameterSet& iConfig) {
207  //now do what ever initialization is needed
208  fileName_ = iConfig.getParameter<std::string>("fileName");
209 
210  file_ = new TFile(fileName_.c_str(), "RECREATE");
211  numCand_ = new TH1F("numCandidates", "Number of candidates found", 10, -0.5, 9.5);
212  numElectrons_ = new TH1F("numElectrons", "Number of Electrons found", 10, -0.5, 9.5);
213  numSuperClusters_ = new TH1F("numSuperClusters", "Number of Ecal SuperClusters", 50, 0, 50);
214 
215  energySuperClusters_ = new TH1F("energySuperClusters", "Super Cluster Energy - all ", 200, 0, 2000.);
216  energySuperClustersEl_ = new TH1F("energySuperClustersEl", "Super Cluster Energy - Electron Cands ", 200, 0., 2000.);
217 
218  sizeSuperClusters_ = new TH1F("sizeSuperClusters", "Super Cluster Size - all ", 20, 0, 19);
219  sizeSuperClustersEl_ = new TH1F("sizeSuperClustersEl", "Super Cluster Size - Electron Cands ", 20, 0, 19);
220 
221  emaxSuperClusters_ = new TH1F("emaxSuperClusters", "Super Cluster Emax - all ", 200, 0, 2000.);
222  emaxSuperClustersEl_ = new TH1F("emaxSuperClustersEl", "Super Cluster Emax - Electron Cands ", 200, 0, 2000.);
223 
224  phiWidthSuperClusters_ = new TH1F("phiWidthSuperClusters", "Super Cluster Width - all ", 20, 0., 0.05);
225  phiWidthSuperClustersEl_ = new TH1F("phiWidthSuperClustersEl", "Super Cluster Width - Electron Cands ", 20, 0., 0.05);
226 
227  ptDiff = new TH1F("ptDiff", " ptDiff ", 20, -10., 10.);
228  pDiff = new TH1F("pDiff", " pDiff ", 100, -50., 50.);
229 
230  pElectronFailed = new TH1F("pElectronFailed", " pElectronFailed ", 55, 0., 110.);
231  ptElectronFailed = new TH1F("ptElectronFailed", " ptElectronFailed ", 55, 0., 110.);
232 
233  pElectronPassed = new TH1F("pElectronPassed", " pElectronPassed ", 55, 0., 110.);
234  ptElectronPassed = new TH1F("ptElectronPassed", " ptElectronPassed ", 55, 0., 110.);
235 
236  sizeSuperClustersFailed = new TH1F("sizeSuperClustersFailed", "Super Cluster Size - Failed ", 20, 0, 19);
237  sizeSuperClustersPassed = new TH1F("sizeSuperClustersPassed", "Super Cluster Size - Passed ", 20, 0, 19);
238 
239  energySuperClustersPassed = new TH1F("energySuperClustersPassed", "Super Cluster Energy - Passed ", 125, 0, 250.);
240  energySuperClustersFailed = new TH1F("energySuperClustersFailed", "Super Cluster Energy - Failed ", 125, 0, 250.);
241 
242  eOverPFailed = new TH1F("eOverPFailed", " E over P - Failed ", 50, 0, 10.);
243  eOverPPassed = new TH1F("eOverPPassed", " E over P - Passed ", 50, 0, 10.);
244 
245  numSiStereoHits_ = new TH1F("numSiStereoHits", "Number of Si StereoHits", 100, 0, 1000);
246  numSiMonoHits_ = new TH1F("numSiMonoHits", "Number of Si MonoHits", 100, 0, 1000);
247  numSiMatchedHits_ = new TH1F("numSiMatchedHits", "Number of Si MatchedHits", 100, 0, 1000);
248 
250 
251  mctruthProducer_ = iConfig.getParameter<std::string>("mctruthProducer");
252  mctruthCollection_ = iConfig.getParameter<std::string>("mctruthCollection");
253 
254  superClusterProducer_ = iConfig.getParameter<std::string>("superClusterProducer");
255  superClusterCollection_ = iConfig.getParameter<std::string>("superClusterCollection");
256 
257  eBRecHitProducer_ = iConfig.getParameter<std::string>("recHitProducer");
258  eBRecHitCollection_ = iConfig.getParameter<std::string>("recHitCollection");
259 
260  siElectronProducer_ = iConfig.getParameter<std::string>("siElectronProducer");
261  siElectronCollection_ = iConfig.getParameter<std::string>("siElectronCollection");
262 
263  electronProducer_ = iConfig.getParameter<std::string>("electronProducer");
264  electronCollection_ = iConfig.getParameter<std::string>("electronCollection");
265 
266  siHitProducer_ = iConfig.getParameter<std::string>("siHitProducer");
267  siRphiHitCollection_ = iConfig.getParameter<std::string>("siRphiHitCollection");
268  siStereoHitCollection_ = iConfig.getParameter<std::string>("siStereoHitCollection");
269  siMatchedHitCollection_ = iConfig.getParameter<std::string>("siMatchedHitCollection");
270 }
271 
272 // SiStripElectronAnalyzer::SiStripElectronAnalyzer(const SiStripElectronAnalyzer& rhs)
273 // {
274 // // do actual copying here;
275 // }
276 
277 SiStripElectronAnalyzer::~SiStripElectronAnalyzer() {
278  // do anything here that needs to be done at desctruction time
279  // (e.g. close files, deallocate resources etc.)
280 
281  file_->Write();
282  file_->Close();
283 }
284 
285 //
286 // assignment operators
287 //
288 // const SiStripElectronAnalyzer& SiStripElectronAnalyzer::operator=(const SiStripElectronAnalyzer& rhs)
289 // {
290 // //An exception safe implementation is
291 // SiStripElectronAnalyzer temp(rhs);
292 // swap(rhs);
293 //
294 // return *this;
295 // }
296 
297 //
298 // member functions
299 //
300 // init for TTree
301 void SiStripElectronAnalyzer::beginJob() {
302  myTree_ = new TTree("myTree", "my first Tree example");
303 
304  myTree_->Branch("NShowers", &NShowers_, "NShowers/I");
305 
306  // first specify the ECAL clusters
307  // need to explicitly include array length.
308  myTree_->Branch("EShower", &EShower_, "EShower[1000]/F");
309  myTree_->Branch("XShower", &XShower_, "XShower[1000]/F");
310  myTree_->Branch("YShower", &YShower_, "YShower[1000]/F");
311  myTree_->Branch("ZShower", &ZShower_, "ZShower[1000]/F");
312 
313  // second specify the Si Stereo Hits
314  myTree_->Branch("NStereoHits", &NStereoHits_, "NStereoHits/I");
315  myTree_->Branch("StereoHitX", &StereoHitX_, "StereoHitX[1000]/F");
316  myTree_->Branch("StereoHitY", &StereoHitY_, "StereoHitY[1000]/F");
317  myTree_->Branch("StereoHitZ", &StereoHitZ_, "StereoHitZ[1000]/F");
318 
319  myTree_->Branch("StereoHitR", &StereoHitR_, "StereoHitR[1000]/F");
320  myTree_->Branch("StereoHitPhi", &StereoHitPhi_, "StereoHitPhi[1000]/F");
321  myTree_->Branch("StereoHitTheta", &StereoHitTheta_, "StereoHitTheta[1000]/F");
322 
323  myTree_->Branch("StereoHitSigX", &StereoHitSigX_, "StereoHitSigX[1000]/F");
324  myTree_->Branch("StereoHitSigY", &StereoHitSigY_, "StereoHitSigY[1000]/F");
325  myTree_->Branch("StereoHitCorr", &StereoHitCorr_, "StereoHitCorr[1000]/F");
326 
327  myTree_->Branch("StereoHitSignal", &StereoHitSignal_, "StereoHitSignal[1000]/F");
328  myTree_->Branch("StereoHitNoise", &StereoHitNoise_, "StereoHitNoise[1000]/F");
329  myTree_->Branch("StereoHitWidth", &StereoHitWidth_, "StereoHitWidth[1000]/I");
330 
331  myTree_->Branch("StereoDetector", &StereoDetector_, "StereoDetector[1000]/I");
332  myTree_->Branch("StereoLayer", &StereoLayer_, "StereoLayer[1000]/I");
333 
334  // specify the Si mono (rphi) hits
335  myTree_->Branch("NMonoHits", &NMonoHits_, "NMonoHits/I");
336  myTree_->Branch("MonoHitX", &MonoHitX_, "MonoHitX[1000]/F");
337  myTree_->Branch("MonoHitY", &MonoHitY_, "MonoHitY[1000]/F");
338  myTree_->Branch("MonoHitZ", &MonoHitZ_, "MonoHitZ[1000]/F");
339 
340  myTree_->Branch("MonoHitR", &MonoHitR_, "MonoHitR[1000]/F");
341  myTree_->Branch("MonoHitPhi", &MonoHitPhi_, "MonoHitPhi[1000]/F");
342  myTree_->Branch("MonoHitTheta", &MonoHitTheta_, "MonoHitTheta[1000]/F");
343 
344  myTree_->Branch("MonoHitSigX", &MonoHitSigX_, "MonoHitSigX[1000]/F");
345  myTree_->Branch("MonoHitSigY", &MonoHitSigY_, "MonoHitSigY[1000]/F");
346  myTree_->Branch("MonoHitCorr", &MonoHitCorr_, "MonoHitCorr[1000]/F");
347 
348  myTree_->Branch("MonoHitSignal", &MonoHitSignal_, "MonoHitSignal[1000]/F");
349  myTree_->Branch("MonoHitNoise", &MonoHitNoise_, "MonoHitNoise[1000]/F");
350  myTree_->Branch("MonoHitWidth", &MonoHitWidth_, "MonoHitWidth[1000]/I");
351 
352  myTree_->Branch("MonoDetector", &MonoDetector_, "MonoDetector[1000]/I");
353  myTree_->Branch("MonoLayer", &MonoLayer_, "MonoLayer[1000]/I");
354 
355  // specify the Si matched (rphi) hits
356  myTree_->Branch("NMatchedHits", &NMatchedHits_, "NMatchedHits/I");
357  myTree_->Branch("MatchedHitX", &MatchedHitX_, "MatchedHitX[1000]/F");
358  myTree_->Branch("MatchedHitY", &MatchedHitY_, "MatchedHitY[1000]/F");
359  myTree_->Branch("MatchedHitZ", &MatchedHitZ_, "MatchedHitZ[1000]/F");
360 
361  myTree_->Branch("MatchedHitR", &MatchedHitR_, "MatchedHitR[1000]/F");
362  myTree_->Branch("MatchedHitPhi", &MatchedHitPhi_, "MatchedHitPhi[1000]/F");
363  myTree_->Branch("MatchedHitTheta", &MatchedHitTheta_, "MatchedHitTheta[1000]/F");
364 
365  myTree_->Branch("MatchedHitSigX", &MatchedHitSigX_, "MatchedHitSigX[1000]/F");
366  myTree_->Branch("MatchedHitSigY", &MatchedHitSigY_, "MatchedHitSigY[1000]/F");
367  myTree_->Branch("MatchedHitCorr", &MatchedHitCorr_, "MatchedHitCorr[1000]/F");
368 
369  myTree_->Branch("MatchedHitSignal", &MatchedHitSignal_, "MatchedHitSignal[1000]/F");
370  myTree_->Branch("MatchedHitNoise", &MatchedHitNoise_, "MatchedHitNoise[1000]/F");
371  myTree_->Branch("MatchedHitWidth", &MatchedHitWidth_, "MatchedHitWidth[1000]/I");
372 
373  myTree_->Branch("MatchedDetector", &MatchedDetector_, "MatchedDetector[1000]/I");
374  myTree_->Branch("MatchedLayer", &MatchedLayer_, "MatchedLayer[1000]/I");
375 }
376 
377 void SiStripElectronAnalyzer::initNtuple() {
378  LogDebug("") << " In initNtuple ";
379 
380  NShowers_ = -999;
381  for (int init = 0; init < myMaxHits; ++init) {
382  EShower_[init] = -999.;
383  XShower_[init] = -999.;
384  YShower_[init] = -999.;
385  ZShower_[init] = -999.;
386  }
387  NStereoHits_ = -999;
388 
389  for (int init = 0; init < myMaxHits; ++init) {
390  StereoHitX_[init] = -999.;
391  StereoHitY_[init] = -999.;
392  StereoHitZ_[init] = -999.;
393  StereoHitR_[init] = -999.;
394  StereoHitPhi_[init] = -999.;
395  StereoHitTheta_[init] = -999.;
396 
397  StereoHitSignal_[init] = -999.;
398  StereoHitNoise_[init] = -999.;
399  StereoHitWidth_[init] = -999;
400  ;
401  }
402 
403  NMonoHits_ = -999;
404  for (int init = 0; init < myMaxHits; ++init) {
405  MonoHitX_[init] = -999.;
406  MonoHitY_[init] = -999.;
407  MonoHitZ_[init] = -999.;
408  MonoHitR_[init] = -999.;
409  MonoHitPhi_[init] = -999.;
410  MonoHitTheta_[init] = -999.;
411 
412  MonoHitSignal_[init] = -999.;
413  MonoHitNoise_[init] = -999.;
414  MonoHitWidth_[init] = -999;
415  ;
416  }
417 
418  NMatchedHits_ = -999;
419  for (int init = 0; init < myMaxHits; ++init) {
420  MatchedHitX_[init] = -999.;
421  MatchedHitY_[init] = -999.;
422  MatchedHitZ_[init] = -999.;
423  MatchedHitR_[init] = -999.;
424  MatchedHitPhi_[init] = -999.;
425  MatchedHitTheta_[init] = -999.;
426 
427  MatchedHitSignal_[init] = -999.;
428  MatchedHitNoise_[init] = -999.;
429  MatchedHitWidth_[init] = -999;
430  ;
431  }
432 }
433 
434 // ------------ method called to produce the data ------------
435 void SiStripElectronAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
436  //Retrieve tracker topology from geometry
437  edm::ESHandle<TrackerTopology> tTopo;
438  iSetup.get<TrackerTopologyRcd>().get(tTopo);
439 
440  using namespace std; // so you can say "cout" and "endl"
441 
442  initNtuple();
443 
444  // http://cmsdoc.cern.ch/swdev/lxr/CMSSW/source/clhep/CLHEP/HepMC/GenParticle.h
445  // http://cmsdoc.cern.ch/swdev/lxr/CMSSW/source/clhep/CLHEP/HepMC/GenVertex.h
446  // removed by JED - causes trouble in release post 0_9_0
447  // edm::Handle<edm::HepMCProduct> mctruthHandle;
448  // iEvent.getByLabel(mctruthProducer_, mctruthCollection_, mctruthHandle);
449  // HepMC::GenEvent mctruth = mctruthHandle->getHepMCData();
450 
451  // for (HepMC::GenEvent::particle_const_iterator partIter = mctruth.particles_begin();
452  // partIter != mctruth.particles_end();
453  // ++partIter) {
454  // // for (HepMC::GenEvent::vertex_const_iterator vertIter = mctruth.vertices_begin();
455  // // vertIter != mctruth.vertices_end();
456  // // ++vertIter) {
457  // CLHEP::HepLorentzVector creation = (*partIter)->CreationVertex();
458  // CLHEP::HepLorentzVector momentum = (*partIter)->Momentum();
459  // HepPDT::ParticleID id = (*partIter)->particleID(); // electrons and positrons are 11 and -11
460  // edm::LogInfo("") << "MC particle id " << id.pid() << ", creationVertex " << creation << " cm, initialMomentum " << momentum << " GeV/c" << endl;
461  // }
462 
463  // load the rechits for the Ecal
464  edm::Handle<EcalRecHitCollection> pRecHits;
465  iEvent.getByLabel(eBRecHitProducer_, eBRecHitCollection_, pRecHits);
466  // Create a pointer to the RecHits - unused for now
467  // const EcalRecHitCollection *hitCollection = pRecHits.product();
468 
469  // http://cmsdoc.cern.ch/swdev/lxr/CMSSW/source/self/DataFormats/EgammaReco/interface/SuperCluster.h
470  edm::Handle<reco::SuperClusterCollection> clusterHandle;
471  iEvent.getByLabel(superClusterProducer_, superClusterCollection_, clusterHandle);
472 
475 
476  LogDebug("") << " Start loop over " << clusterHandle->end() - clusterHandle->begin() << " superClusters ";
477 
478  for (reco::SuperClusterCollection::const_iterator clusterIter = clusterHandle->begin();
479  clusterIter != clusterHandle->end();
480  ++clusterIter) {
481  double energy = clusterIter->energy();
482  math::XYZPoint position = clusterIter->position();
483  std::ostringstream str;
484 
485  str << " SuperCluster " << energy << " GeV, position " << position << " cm"
486  << "\n";
487 
488  energySuperClusters_->Fill(energy);
489  sizeSuperClusters_->Fill(clusterIter->clustersSize());
490  // this only makes sense for hybrid superclusters
491 
492  // try to point to the constituent clusters for this SuperCluster
493 
494  str << "About to loop over basicClusters"
495  << "\n";
496 
497  double emaxSuperCluster = 0.;
498  double phibar = 0.;
499  double phi2bar = 0.;
500  double eTotSuperCluster = 0.;
501 
502  for (reco::CaloCluster_iterator basicClusterIter = clusterIter->clustersBegin();
503  basicClusterIter != clusterIter->clustersEnd();
504  ++basicClusterIter) {
505  //std::vector<DetId> theIds= (*basicClusterIter)->getHitsByDetId();
506 
507  str << " basicCluster Energy " << (*basicClusterIter)->energy() << " Position " << (*basicClusterIter)->position()
508  << " \n"
509  << " Position phi " << (*basicClusterIter)->position().phi() << " recHits "
510  << (*basicClusterIter)->size() << " \n";
511 
512  double eCluster = (*basicClusterIter)->energy();
513  if (eCluster > emaxSuperCluster) {
514  emaxSuperCluster = eCluster;
515  }
516  eTotSuperCluster += eCluster;
517  double phiCluster = (*basicClusterIter)->position().phi();
518  phibar += eCluster * phiCluster;
519  phi2bar += eCluster * phiCluster * phiCluster;
520 
521  } // end of basicClusterIter loop
522 
523  phibar = phibar / eTotSuperCluster;
524  phi2bar = phi2bar / eTotSuperCluster;
525  double phiWidth = phi2bar - phibar * phibar;
526  if (phiWidth > 0.) {
527  phiWidth = std::pow(phiWidth, 0.5);
528  } else {
529  phiWidth = 0.;
530  }
531  str << " SuperCluster stats "
532  << "\n";
533  str << "phibar " << phibar << " phi2bar " << phi2bar << " eTotSuperCluster " << eTotSuperCluster << " phiWidth "
534  << phiWidth << std::endl;
535 
536  phiWidthSuperClusters_->Fill(phiWidth);
537 
538  emaxSuperClusters_->Fill(emaxSuperCluster);
539 
540  str << " Done with this SuperCluster " << std::endl;
541 
542  LogDebug("") << str.str();
543 
544  } // end of loop over superClusters
545 
546  LogDebug("") << " End loop over superClusters ";
547 
550 
552  //
553  // loop over all EcalRecHits and print out their x,y,z,E
554  // edm::LogInfo("") << " Dumping all recHits in this event " << endl ;
555  // for(EcalRecHitCollection::const_iterator _blah = hitCollection->begin();
556  // _blah != hitCollection->end() ; ++_blah ) {
557  // edm::LogInfo("") << "Ecal RecHit Energy: " << _blah->energy() << endl ;
558  // // " Position " << _blah.position() << endl ;
559  // }
560  //
561  // edm::LogInfo("") << "Dump finished " << endl ;
562  //
564 
565  // DataFormats/EgammaCandidates/src/SiStripElectron.cc
566  edm::Handle<reco::SiStripElectronCollection> siStripElectronHandle;
567  iEvent.getByLabel(siElectronProducer_, siElectronCollection_, siStripElectronHandle);
568 
570 
571  LogDebug("") << " Dumping Algo's guess of SiStripElectron Candidate Info ";
572  int numberOfElectrons = 0;
573  // need to check if fit succeeded
574  LogDebug("") << " Number of SiStripElectrons " << siStripElectronHandle->size();
575 
576  for (reco::SiStripElectronCollection::const_iterator electronIter = siStripElectronHandle->begin();
577  electronIter != siStripElectronHandle->end();
578  ++electronIter) {
579  LogDebug("") << "about to get stuff from electroncandidate " << numberOfElectrons << "\n"
580  << "supercluster energy = " << electronIter->superCluster()->energy() << "\n"
581  << "fit results are phi(r) = " << electronIter->phiAtOrigin() << " + " << electronIter->phiVsRSlope()
582  << "*r"
583  << "\n"
584  << " chi2 " << electronIter->chi2() << " ndof " << electronIter->ndof() << "\n"
585  << " Pt " << electronIter->pt() << "\n"
586  << "P, Px, Py, Pz " << electronIter->p() << " " << electronIter->px() << " " << electronIter->py()
587  << " " << electronIter->pz() << "\n"
588  << "you get the idea...";
589 
590  // make plots for supercluster that an electron has been associ w/. here
591  energySuperClustersEl_->Fill(electronIter->superCluster()->energy());
592  sizeSuperClustersEl_->Fill(electronIter->superCluster()->clustersSize());
593 
594  // loop over basicClusters to get energy
595  double emaxSuperCluster = 0.;
596  double phibar = 0.;
597  double phi2bar = 0.;
598  double eTotSuperCluster = 0.;
599 
600  for (reco::CaloCluster_iterator basicClusterIter = electronIter->superCluster()->clustersBegin();
601  basicClusterIter != electronIter->superCluster()->clustersEnd();
602  ++basicClusterIter) {
603  //std::vector<DetId> theIds= (*basicClusterIter)->getHitsByDetId();
604 
605  double eCluster = (*basicClusterIter)->energy();
606  if (eCluster > emaxSuperCluster) {
607  emaxSuperCluster = eCluster;
608  }
609  eTotSuperCluster += eCluster;
610  double phiCluster = (*basicClusterIter)->position().phi();
611  phibar += eCluster * phiCluster;
612  phi2bar += eCluster * phiCluster * phiCluster;
613  }
614 
615  phibar = phibar / eTotSuperCluster;
616  phi2bar = phi2bar / eTotSuperCluster;
617  double phiWidth = phi2bar - phibar * phibar;
618  if (phiWidth > 0.) {
619  phiWidth = std::pow(phiWidth, 0.5);
620  } else {
621  phiWidth = 0.;
622  }
623 
624  phiWidthSuperClustersEl_->Fill(phiWidth);
625 
626  emaxSuperClustersEl_->Fill(emaxSuperCluster);
627 
628  numberOfElectrons++;
629  }
630 
631  numCand_->Fill(siStripElectronHandle->size());
632 
634 
635  // Now loop over the electrons (ie the fitted things.)
636 
637  LogDebug("") << " About to check Electrons";
638 
639  edm::Handle<reco::ElectronCollection> electrons;
640  iEvent.getByLabel(electronProducer_, electronCollection_, electrons);
641 
642  numElectrons_->Fill(electrons->end() - electrons->begin());
643 
644  // set up vector of bool for SiStrips having or not having Electrons
645  // this causes a warning because of variable array size at compilation time ;
646  // BAD bool hasElectron_[siStripElectronHandle->end()- siStripElectronHandle->begin()] ;
647  bool* hasElectron_ = new bool[siStripElectronHandle->end() - siStripElectronHandle->begin()];
648  for (int icount = 0; icount < siStripElectronHandle->end() - siStripElectronHandle->begin(); ++icount) {
649  hasElectron_[icount] = false;
650  }
651 
652  // also set up a counter to associate the ith electron to the jth strippy
653  // Electron_to_strippy[i] = j: i-th Electron is j-th strippy
654  // BAD unsigned int Electron_to_strippy[electrons->end()- electrons->begin()];
655  unsigned int* Electron_to_strippy = new unsigned int[electrons->end() - electrons->begin()];
656  for (int icount = 0; icount < electrons->end() - electrons->begin(); ++icount) {
657  Electron_to_strippy[icount] = 0;
658  }
659 
660  unsigned int ecount = 0;
661  for (reco::ElectronCollection::const_iterator electronIter = electrons->begin(); electronIter != electrons->end();
662  ++electronIter) {
663  LogDebug("") << " Associating Electrons to Strippies ";
664  LogDebug("") << " PT is " << electronIter->track()->pt();
665 
666  reco::TrackRef tr = (*electronIter).track();
667  uint32_t id = (*electronIter->track()->recHitsBegin())->geographicalId().rawId();
668  LocalPoint pos = (*electronIter->track()->recHitsBegin())->localPosition();
669 
670  unsigned int icount = 0;
671  LogDebug("") << " About to loop over Strippies "
672  << " \n "
673  << " icount " << icount << " max " << siStripElectronHandle->end() - siStripElectronHandle->begin();
674 
675  for (reco::SiStripElectronCollection::const_iterator strippyiter = siStripElectronHandle->begin();
676  strippyiter != siStripElectronHandle->end();
677  ++strippyiter) {
678  bool hitInCommon = false;
679  // loop over rphi hits
680  for (std::vector<SiStripRecHit2D>::const_iterator hiter = strippyiter->rphiRecHits().begin();
681  hiter != strippyiter->rphiRecHits().end();
682  ++hiter) {
683  if (hiter->geographicalId().rawId() == id && (hiter->localPosition() - pos).mag() < 1e-10) {
684  hitInCommon = true;
685  break;
686  }
687  }
688 
689  for (std::vector<SiStripRecHit2D>::const_iterator hiter = strippyiter->stereoRecHits().begin();
690  hiter != strippyiter->stereoRecHits().end();
691  ++hiter) {
692  if (hiter->geographicalId().rawId() == id && (hiter->localPosition() - pos).mag() < 1e-10) {
693  hitInCommon = true;
694  break;
695  }
696  }
697  if (hitInCommon) { //this Electron belongs to this SiStripElectron.
698  hasElectron_[icount] = true;
699  Electron_to_strippy[ecount] = icount;
700  ptDiff->Fill(std::abs(electronIter->track()->pt()) - std::abs(strippyiter->pt()));
701  pDiff->Fill(std::abs(electronIter->track()->p()) - std::abs(strippyiter->p()));
702  }
703  icount++;
704  } // Sistrip loop
705  ecount++;
706  } // Electrons
707 
708  LogDebug("") << " Done looping over Electrons ";
709 
710  unsigned int counter = 0;
711  for (reco::SiStripElectronCollection::const_iterator strippyIter = siStripElectronHandle->begin();
712  strippyIter != siStripElectronHandle->end();
713  ++strippyIter) {
714  bool skipThis = !hasElectron_[counter];
715  if (skipThis) {
716  // plot stuff for SIStripElectrons that don't have fits associated
717 
718  LogDebug("") << " SiStrip Failed Electron "
719  << " \n "
720  << " p " << strippyIter->p() << " \n "
721  << " pt " << strippyIter->pt() << " \n "
722  << " SuperClust size " << strippyIter->superCluster()->clustersSize();
723 
724  pElectronFailed->Fill(std::abs(strippyIter->p()));
725  ptElectronFailed->Fill(std::abs(strippyIter->pt()));
726  sizeSuperClustersFailed->Fill(strippyIter->superCluster()->clustersSize());
727  LogDebug("") << " done filling Failed histos ";
728  // energySuperClustersFailed->Fill(strippyIter->superCluster()->energy());
729  // if(strippyIter->p()>0.) {
730  // eOverPFailed->Fill(strippyIter->superCluster()->energy()/strippyIter->p());
731  // }else {
732  // eOverPFailed->Fill(-1.0);
733  // }
734 
735  } else {
736  LogDebug("") << " SiStrip Passed Electron "
737  << " \n "
738  << " p " << strippyIter->p() << " \n "
739  << " pt " << strippyIter->pt() << " \n "
740  << " SuperClust size " << strippyIter->superCluster()->clustersSize();
741  pElectronPassed->Fill(std::abs(strippyIter->p()));
742  ptElectronPassed->Fill(std::abs(strippyIter->pt()));
743  sizeSuperClustersPassed->Fill(strippyIter->superCluster()->clustersSize());
744  LogDebug("") << " done filling passed histos ";
745  // energySuperClustersPassed->Fill(strippyIter->superCluster()->energy());
746  // if(strippyIter->p()>0.) {
747  // eOverPPassed->Fill(strippyIter->superCluster()->energy()/strippyIter->p());
748  // }else {
749  // eOverPPassed->Fill(-1.0);
750  // }
751 
752  } // skipThis
753  counter++;
754  }
755 
756  LogDebug("") << "Dump info for all electrons ";
757 
758  for (reco::ElectronCollection::const_iterator electronIter1 = electrons->begin(); electronIter1 != electrons->end();
759  ++electronIter1) {
760  reco::TrackRef tr1 = (*electronIter1).track();
761  // let's find its associated SiStripElectron and SuperCluster
762  unsigned int ecount1 = electronIter1 - electrons->begin();
763  unsigned int stripCount1 = 0;
764  reco::SiStripElectronCollection::const_iterator strippyIter1;
765  for (reco::SiStripElectronCollection::const_iterator strippyIter = siStripElectronHandle->begin();
766  strippyIter != siStripElectronHandle->end();
767  ++strippyIter) {
768  if (Electron_to_strippy[ecount1] == stripCount1) {
769  strippyIter1 = strippyIter;
770  break;
771  }
772  stripCount1++;
773  } // strippy loop
774  ecount1++;
775 
776  std::ostringstream str;
777 
778  str << " SiStripElect p , px, py, pz " << strippyIter1->p() << " " << strippyIter1->px() << " "
779  << strippyIter1->py() << " " << strippyIter1->pz() << "\n " << std::endl;
780 
781  str << " Electron p px, py, pz, = " << tr1->p() << " " << tr1->px() << " " << tr1->py() << " " << tr1->pz()
782  << "\n"
783  << std::endl;
784 
785  double EClust1 = strippyIter1->superCluster()->energy();
786  double XClust1 = strippyIter1->superCluster()->x();
787  double YClust1 = strippyIter1->superCluster()->y();
788  double ZClust1 = strippyIter1->superCluster()->z();
789 
790  double rho1 = sqrt(XClust1 * XClust1 + YClust1 * YClust1 + ZClust1 * ZClust1);
791  double costheta1 = ZClust1 / rho1;
792  double sintheta1 = sqrt(1 - costheta1 * costheta1);
793  if (ZClust1 < 0) {
794  sintheta1 = -sintheta1;
795  }
796  double cosphi1 = XClust1 / sqrt(XClust1 * XClust1 + YClust1 * YClust1);
797  double sinphi1 = YClust1 / sqrt(XClust1 * XClust1 + YClust1 * YClust1);
798 
799  str << " Ecal for electron E, px, py, pz " << EClust1 << " " << EClust1 * sintheta1 * cosphi1 << " "
800  << EClust1 * sintheta1 * sinphi1 << " " << EClust1 * costheta1 << "\n"
801  << std::endl;
802 
803  LogDebug("") << str.str();
804 
805  } // loop over electrons
806  LogDebug("") << "Done Dumping info for all electrons ";
807 
809  // LogDebug("")<< " Checking Electrons" ;
810  // LogDebug("")<< " PT is " << electronIter->track()->pt() ;
811  // reco::TrackRef tr =(*electronIter).track();
813  if (electrons->end() - electrons->begin() > 1) {
814  edm::LogInfo("") << " Two electrons in this event " << std::endl;
815  for (reco::ElectronCollection::const_iterator electronIter1 = electrons->begin();
816  electronIter1 != electrons->end() - 1;
817  ++electronIter1) {
818  reco::TrackRef tr1 = (*electronIter1).track();
819 
820  // let's find its associated SiStripElectron and SuperCluster
821  // use the Electron_to_strippy[] array
822  unsigned int ecount1 = electronIter1 - electrons->begin();
823  // loop over strippies to find the corresponding one
824  unsigned int stripCount1 = 0;
825  reco::SiStripElectronCollection::const_iterator strippyIter1;
826  for (reco::SiStripElectronCollection::const_iterator strippyIter = siStripElectronHandle->begin();
827  strippyIter != siStripElectronHandle->end();
828  ++strippyIter) {
829  if (Electron_to_strippy[ecount1] == stripCount1) {
830  strippyIter1 = strippyIter;
831  break;
832  }
833  stripCount1++;
834  } // strippy loop
835 
836  double EClust1 = strippyIter1->superCluster()->energy();
837  double XClust1 = strippyIter1->superCluster()->x();
838  double YClust1 = strippyIter1->superCluster()->y();
839  double ZClust1 = strippyIter1->superCluster()->z();
840 
841  for (reco::ElectronCollection::const_iterator electronIter2 = electronIter1 + 1;
842  electronIter2 != electrons->end();
843  ++electronIter2) {
844  reco::TrackRef tr2 = (*electronIter2).track();
845 
846  unsigned int ecount2 = electronIter2 - electrons->begin();
847  unsigned int stripCount2 = 0;
848  reco::SiStripElectronCollection::const_iterator strippyIter2;
849  for (reco::SiStripElectronCollection::const_iterator strippyIter = siStripElectronHandle->begin();
850  strippyIter != siStripElectronHandle->end();
851  ++strippyIter) {
852  if (Electron_to_strippy[ecount2] == stripCount2) {
853  strippyIter2 = strippyIter;
854  break;
855  }
856  stripCount2++;
857  } // strippy loop
858 
859  double EClust2 = strippyIter2->superCluster()->energy();
860  double XClust2 = strippyIter2->superCluster()->x();
861  double YClust2 = strippyIter2->superCluster()->y();
862  double ZClust2 = strippyIter2->superCluster()->z();
863 
864  // now get supercluster from this:
865 
866  edm::LogInfo("") << " Electron p1 = " << tr1->p() << " p1x " << tr1->px() << " p1y " << tr1->py() << " p1z "
867  << tr1->pz() << std::endl;
868 
869  edm::LogInfo("") << " Electron p2 = " << tr2->p() << " p2x " << tr2->px() << " p2y " << tr2->py() << " p2z "
870  << tr2->pz() << std::endl;
871 
872  // combine the two in an (e,e) pair
873  double Zpx = tr1->px() + tr2->px();
874  double Zpy = tr1->py() + tr2->py();
875  double Zpz = tr1->pz() + tr2->pz();
876  double Ze = std::abs(tr1->p()) + std::abs(tr2->p());
877  edm::LogInfo("") << " Z mass " << sqrt(Ze * Ze - Zpx * Zpx - Zpy * Zpy - Zpz * Zpz) << std::endl;
878 
879  // combine the SuperClusts into a Z
880  double rho1 = sqrt(XClust1 * XClust1 + YClust1 * YClust1 + ZClust1 * ZClust1);
881  double costheta1 = ZClust1 / rho1;
882  double sintheta1 = sqrt(1 - costheta1 * costheta1);
883  if (ZClust1 < 0) {
884  sintheta1 = -sintheta1;
885  }
886  double cosphi1 = XClust1 / sqrt(XClust1 * XClust1 + YClust1 * YClust1);
887  double sinphi1 = YClust1 / sqrt(XClust1 * XClust1 + YClust1 * YClust1);
888 
889  double rho2 = sqrt(XClust2 * XClust2 + YClust2 * YClust2 + ZClust2 * ZClust2);
890  double costheta2 = ZClust2 / rho2;
891  double sintheta2 = sqrt(1 - costheta2 * costheta2);
892  if (ZClust2 < 0) {
893  sintheta2 = -sintheta2;
894  }
895  double cosphi2 = XClust2 / sqrt(XClust2 * XClust2 + YClust2 * YClust2);
896  double sinphi2 = YClust2 / sqrt(XClust2 * XClust2 + YClust2 * YClust2);
897 
898  edm::LogInfo("") << "Energy of supercluster for 1st electron " << EClust1 << " "
899  << EClust1 * sintheta1 * cosphi1 << " " << EClust1 * sintheta1 * sinphi1 << " "
900  << EClust1 * costheta1 << " " << std::endl;
901 
902  edm::LogInfo("") << "Energy of supercluster for 2nd electron " << EClust2 << " "
903  << EClust2 * sintheta2 * cosphi2 << " " << EClust2 * sintheta2 * sinphi2 << " "
904  << EClust2 * costheta2 << " " << std::endl;
905 
906  // get the supercluster pair
907  double Zgpx = EClust1 * sintheta1 * cosphi1 + EClust2 * sintheta2 * cosphi2;
908  double Zgpy = EClust1 * sintheta1 * sinphi1 + EClust2 * sintheta2 * sinphi2;
909  double Zgpz = EClust1 * costheta1 + EClust2 * costheta2;
910  double ZgE = EClust1 + EClust2;
911 
912  edm::LogInfo("") << " Z mass from ECAL " << sqrt(ZgE * ZgE - Zgpx * Zgpx - Zgpy * Zgpy - Zgpz * Zgpz)
913  << std::endl;
914 
915  } //inner loop
916  } // outer loop
917  } // m(ee) loop
918 
919  delete[] hasElectron_;
920  delete[] Electron_to_strippy;
921 
925  LogDebug("") << " About to dump tracker info ";
926 
927  edm::ESHandle<TrackerGeometry> trackerHandle;
928  iSetup.get<TrackerDigiGeometryRecord>().get(trackerHandle);
929 
930  edm::Handle<SiStripRecHit2DCollection> rphiHitsHandle;
931  iEvent.getByLabel(siHitProducer_, siRphiHitCollection_, rphiHitsHandle);
932 
933  edm::Handle<SiStripRecHit2DCollection> stereoHitsHandle;
934  iEvent.getByLabel(siHitProducer_, siStereoHitCollection_, stereoHitsHandle);
935 
936  edm::Handle<SiStripMatchedRecHit2DCollection> matchedHitsHandle;
937  iEvent.getByLabel(siHitProducer_, siMatchedHitCollection_, matchedHitsHandle);
938 
940 
942  NShowers_ = 0;
943  for (reco::SuperClusterCollection::const_iterator clusterIter = clusterHandle->begin();
944  clusterIter != clusterHandle->end();
945  ++clusterIter) {
946  double energy = clusterIter->energy();
947  math::XYZPoint position = clusterIter->position();
948  if (NShowers_ < myMaxHits) {
949  EShower_[NShowers_] = energy;
950  XShower_[NShowers_] = position.x();
951  YShower_[NShowers_] = position.y();
952  ZShower_[NShowers_] = position.z();
953  ++NShowers_;
954  }
955  // Loop over all crystals in this supercluster - see
956  // RecoEcal/EgamaClusterProducers/src/EgammaSimpleAnalyzer.cc
957  // Look also at DataFormats/EgammaReco/interface/SuperCluster.h
958  }
959  numSuperClusters_->Fill(NShowers_);
961 
962  LogDebug("") << " Looping over stereo hits ";
963 
965  int myHits = 0;
966  for (SiStripRecHit2DCollection::DataContainer::const_iterator hit = stereoHitsHandle->data().begin(),
967  hitend = stereoHitsHandle->data().end();
968  hit != hitend;
969  ++hit) {
970  DetId id(hit->geographicalId());
971  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB ||
972  (hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
973  GlobalPoint position = trackerHandle->idToDet(hit->geographicalId())->surface().toGlobal(hit->localPosition());
974  //from RecoLocalTracker/SiStripClusterizer/test/TestCluster.cc
975  // cf also TrackHitAssociator.cc SiStripRecHitMatcher.cc SiStrip1DMeasurementTransformator.cc (KalmanUpdators)
976  SiStripRecHit2D const rechit = *hit;
977  // LocalPoint myposition = rechit.localPosition() ;
978  LocalError myerror = rechit.localPositionError();
979 
980  // Get layer and subdetector ID here for this hit
981  // see SiStripRecHitConverter/test/ValHit.cc
982  Int_t siLayerNum = 0;
983  Int_t siDetNum = 0;
984  string siDetName = "";
985  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB) {
986  // siLayerNum = tTopo->tibLayer(rechit->geographicalID());
987  siLayerNum = tTopo->tibLayer(id);
988  siDetNum = 1;
989  siDetName = "TIB";
990  } else if ((hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
991  siLayerNum = tTopo->tobLayer(id);
992  siDetNum = 2;
993  siDetName = "TOB";
994  // } else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TID ){
995  // // should we use side/wheel/ring/module/stereo() ?
996  // siLayerNum = tTopo->tidWheel(id);
997  // siDetNum = 3 ;
998  // siDetName = "TID" ;
999  // }else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TEC ){
1000  // //choices are side/petal/wheel/ring/module/glued/stereo
1001  // siLayerNum = tTopo->tecWheel(id);
1002  // siDetNum = 4 ;
1003  // siDetName = "TEC" ;
1004  } else {
1005  siLayerNum = -999;
1006  siDetNum = -999;
1007  siDetName = "NULL";
1008  }
1009  // LogDebug("") << siDetName << " " << siLayerNum ;
1010 
1011  const SiStripRecHit2D::ClusterRef& clust = rechit.cluster();
1012  double Signal = 0;
1013  double Noise2 = 0;
1014  int StripCount = 0;
1015  if (clust.isNonnull()) {
1016  // LogDebug("") << " barycenter " << clust->barycenter() ;
1017  // const std::vector<uint16_t> amplitudes=clust->amplitudes();
1018  const auto& amplitudes = clust->amplitudes();
1019  for (size_t i = 0; i < amplitudes.size(); i++) {
1020  Signal += amplitudes[i];
1021  //ignore for now Noise2 +=SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i)*SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i);
1022  StripCount++;
1023  }
1024  } else {
1025  LogDebug("") << " null cluster ";
1026  }
1027  // LogDebug("") << "Signal " << Signal << " Noise2 " << Noise2 << " StripCount " << StripCount ;
1028  // Dump position
1029  // LogDebug("") << " Stereo "
1030  // << "local position: "<<myposition.x()<<" "
1031  // << myposition.y()<<" "<<myposition.z()<<"\n"
1032  // << "local error: "<<myerror.xx()<<" "
1033  // << myerror.xy()<<" "<<myerror.yy() << "\n"
1034  // << "global position: " << position.x() << " "
1035  // << position.y()<<" "<< position.z()<<"\n"
1036  // << " siDetNum " << siDetNum
1037  // << " siLayerNum " << siLayerNum ;
1038 
1039  if (myHits < myMaxHits) {
1040  StereoHitX_[myHits] = position.x();
1041  StereoHitY_[myHits] = position.y();
1042  StereoHitZ_[myHits] = position.z();
1043 
1044  StereoHitR_[myHits] = position.perp();
1045  StereoHitPhi_[myHits] = position.phi();
1046  StereoHitTheta_[myHits] = position.theta();
1047 
1048  StereoHitSigX_[myHits] = sqrt(myerror.xx());
1049  StereoHitSigY_[myHits] = sqrt(myerror.yy());
1050  StereoHitCorr_[myHits] = myerror.xy() / sqrt(myerror.xx() * myerror.yy());
1051 
1052  StereoHitSignal_[myHits] = Signal;
1053  StereoHitNoise_[myHits] = Noise2;
1054  StereoHitWidth_[myHits] = StripCount;
1055 
1056  StereoDetector_[myHits] = siDetNum;
1057  StereoLayer_[myHits] = siLayerNum;
1058 
1059  ++myHits;
1060  }
1061  } // end if this is the right subdetector
1062  } // end loop over hits
1063  NStereoHits_ = myHits;
1064 
1065  numSiStereoHits_->Fill(NStereoHits_);
1066 
1069 
1070  LogDebug("") << " Looping over Mono Hits ";
1072  myHits = 0;
1073  for (SiStripRecHit2DCollection::DataContainer::const_iterator hit = rphiHitsHandle->data().begin(),
1074  hitend = rphiHitsHandle->data().end();
1075  hit != hitend;
1076  ++hit) {
1077  DetId id(hit->geographicalId());
1078 
1079  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB ||
1080  (hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
1081  GlobalPoint position = trackerHandle->idToDet(hit->geographicalId())->surface().toGlobal(hit->localPosition());
1082  //from RecoLocalTracker/SiStripClusterizer/test/TestCluster.cc
1083  // cf also TrackHitAssociator.cc SiStripRecHitMatcher.cc SiStrip1DMeasurementTransformator.cc (KalmanUpdators)
1084  SiStripRecHit2D const rechit = *hit;
1085  // LocalPoint myposition = rechit.localPosition() ;
1086  LocalError myerror = rechit.localPositionError();
1087 
1088  // Get layer and subdetector ID here for this hit
1089  // see SiStripRecHitConverter/test/ValHit.cc
1090  Int_t siLayerNum = 0;
1091  Int_t siDetNum = 0;
1092  string siDetName = "";
1093  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB) {
1094  // siLayerNum = tTopo->tibLayer(rechit->geographicalID());
1095  siLayerNum = tTopo->tibLayer(id);
1096  siDetNum = 1;
1097  siDetName = "TIB";
1098  } else if ((hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
1099  siLayerNum = tTopo->tobLayer(id);
1100  siDetNum = 2;
1101  siDetName = "TOB";
1102  // } else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TID ){
1103  // // should we use side/wheel/ring/module/stereo() ?
1104  // siLayerNum = tTopo->tidWheel(id);
1105  // siDetNum = 3 ;
1106  // siDetName = "TID" ;
1107  // }else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TEC ){
1108  // //choices are side/petal/wheel/ring/module/glued/stereo
1109  // siLayerNum = tTopo->tecWheel(id);
1110  // siDetNum = 4 ;
1111  // siDetName = "TEC"
1112  ;
1113  } else {
1114  siLayerNum = -999;
1115  siDetNum = -999;
1116  siDetName = "NULL";
1117  }
1118  // LogDebug("") << siDetName << " " << siLayerNum ;
1119  const SiStripRecHit2D::ClusterRef& clust = rechit.cluster();
1120  double Signal = 0;
1121  double Noise2 = 0;
1122  int StripCount = 0;
1123  if (clust.isNonnull()) {
1124  // LogDebug("") << " barycenter " << clust->barycenter() ;
1125  // const std::vector<uint16_t> amplitudes=clust->amplitudes();
1126  const auto& amplitudes = clust->amplitudes();
1127  for (size_t i = 0; i < amplitudes.size(); i++) {
1128  Signal += amplitudes[i];
1129  //ignore for now Noise2 +=SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i)*SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i);
1130  StripCount++;
1131  }
1132  } else {
1133  LogDebug("") << " null cluster ";
1134  }
1135  // LogDebug("") << "Signal " << Signal << " Noise2 " << Noise2 << " StripCount " << StripCount ;
1136 
1137  // Dump position info
1138  // LogDebug("") << " Mono "
1139  // << "local position: "<<myposition.x()<<" "
1140  // << myposition.y()<<" "<<myposition.z()<<"\n"
1141  // <<"local error: "<<myerror.xx()<<" "
1142  // << myerror.xy()<<" "<<myerror.yy() << "\n"
1143  // << "global position: " << position.x() << " "
1144  // << position.y()<<" "<< position.z()<<"\n"
1145  // << " siDetNum " << siDetNum
1146  // << " siLayerNum " << siLayerNum ;
1147 
1148  if (myHits < myMaxHits) {
1149  MonoHitX_[myHits] = position.x();
1150  MonoHitY_[myHits] = position.y();
1151  MonoHitZ_[myHits] = position.z();
1152 
1153  MonoHitR_[myHits] = position.perp();
1154  MonoHitPhi_[myHits] = position.phi();
1155  MonoHitTheta_[myHits] = position.theta();
1156 
1157  MonoHitSigX_[myHits] = sqrt(myerror.xx());
1158  MonoHitSigY_[myHits] = sqrt(myerror.yy());
1159  MonoHitCorr_[myHits] = myerror.xy() / sqrt(myerror.xx() * myerror.yy());
1160 
1161  MonoHitSignal_[myHits] = Signal;
1162  MonoHitNoise_[myHits] = Noise2;
1163  MonoHitWidth_[myHits] = StripCount;
1164 
1165  MonoDetector_[myHits] = siDetNum;
1166  MonoLayer_[myHits] = siLayerNum;
1167 
1168  ++myHits;
1169  } // of if(myHits < myMaxHits)
1170  // LogDebug("")<< "end of myHits < myMaxHits " ;
1171  } // end if this is the right subdetector
1172  // LogDebug("")<< "end of TIB/TOB check " ;
1173  } // end loop over hits
1174  // LogDebug("")<< " end of loop over hits " ;
1175  NMonoHits_ = myHits;
1176 
1177  numSiMonoHits_->Fill(NMonoHits_);
1178 
1181 
1182  LogDebug("") << " Loop over Matched Hits ";
1183 
1185  myHits = 0;
1186  for (SiStripMatchedRecHit2DCollection::DataContainer::const_iterator hit = matchedHitsHandle->data().begin(),
1187  hitend = matchedHitsHandle->data().end();
1188  hit != hitend;
1189  ++hit) {
1190  DetId id(hit->geographicalId());
1191  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB ||
1192  (hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
1193  GlobalPoint position = trackerHandle->idToDet(hit->geographicalId())->surface().toGlobal(hit->localPosition());
1194  SiStripMatchedRecHit2D const rechit = *hit;
1195  // LocalPoint myposition = rechit.localPosition() ;
1196  LocalError myerror = rechit.localPositionError();
1197 
1198  // Get layer and subdetector ID here for this hit
1199  // see SiStripRecHitConverter/test/ValHit.cc
1200  Int_t siLayerNum = 0;
1201  Int_t siDetNum = 0;
1202  string siDetName = "";
1203  if ((hit->geographicalId()).subdetId() == StripSubdetector::TIB) {
1204  siLayerNum = tTopo->tibLayer(id);
1205  siDetNum = 1;
1206  siDetName = "TIB";
1207  } else if ((hit->geographicalId()).subdetId() == StripSubdetector::TOB) {
1208  siLayerNum = tTopo->tobLayer(id);
1209  siDetNum = 2;
1210  siDetName = "TOB";
1211  // } else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TID ){
1212  // // should we use side/wheel/ring/module/stereo() ?
1213  // siLayerNum = tTopo->tidWheel(id);
1214  // siDetNum = 3 ;
1215  // siDetName = "TID" ;
1216  // }else if ( (hit->geographicalId()).subdetId() == StripSubdetector::TEC ){
1217  // //choices are side/petal/wheel/ring/module/glued/stereo
1218  // siLayerNum = tTopo->tecWheel(id);
1219  // siDetNum = 4 ;
1220  // siDetName = "TEC" ;
1221  } else {
1222  siLayerNum = -999;
1223  siDetNum = -999;
1224  siDetName = "NULL";
1225  }
1226  // const edm::Ref<edm::DetSetVector<SiStripCluster>, SiStripCluster, edm::refhelper::FindForDetSetVector<SiStripCluster> > clust=rechit.cluster();
1227  double Signal = 0;
1228  double Noise2 = 0;
1229  int StripCount = 0;
1230  // if(clust.isNonnull()) {
1231  // LogDebug("") << " barycenter " << clust->barycenter() ;
1232  // const std::vector<uint16_t> amplitudes=clust->amplitudes();
1233  // for(size_t i = 0 ; i<amplitudes.size(); i++ ){
1234  // Signal +=amplitudes[i] ;
1235  // //ignore for now Noise2 +=SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i)*SiStripNoiseService_.getNoise(detid,clust->firstStrip()+i);
1236  // StripCount++;
1237  // }
1238  // } else {
1239  // LogDebug("") << " null cluster " ;
1240  // }
1241  // LogDebug("") << "Signal " << Signal << " Noise2 " << Noise2 << " StripCount " << StripCount ;
1242 
1243  // Dump position info
1244  // LogDebug("") << " Matched "
1245  // << "local position: "<<myposition.x()<<" "
1246  // << myposition.y()<<" "<<myposition.z()<<"\n"
1247  // << "local error: "<<myerror.xx()<<" "
1248  // << myerror.xy()<<" "<<myerror.yy() << "\n"
1249  // << "global position: " << position.x() << " "
1250  // << position.y()<<" "<< position.z()<<"\n"
1251  // << " siDetNum " << siDetNum
1252  // << " siLayerNum " << siLayerNum ;
1253 
1254  if (myHits < myMaxHits) {
1255  MatchedHitX_[myHits] = position.x();
1256  MatchedHitY_[myHits] = position.y();
1257  MatchedHitZ_[myHits] = position.z();
1258 
1259  MatchedHitR_[myHits] = position.perp();
1260  MatchedHitPhi_[myHits] = position.phi();
1261  MatchedHitTheta_[myHits] = position.theta();
1262 
1263  MatchedHitSigX_[myHits] = sqrt(myerror.xx());
1264  MatchedHitSigY_[myHits] = sqrt(myerror.yy());
1265  MatchedHitCorr_[myHits] = myerror.xy() / sqrt(myerror.xx() * myerror.yy());
1266 
1267  MatchedHitSignal_[myHits] = Signal;
1268  MatchedHitNoise_[myHits] = Noise2;
1269  MatchedHitWidth_[myHits] = StripCount;
1270 
1271  MatchedDetector_[myHits] = siDetNum;
1272  MatchedLayer_[myHits] = siLayerNum;
1273 
1274  ++myHits;
1275  }
1276  } // end if this is the right subdetector (TIB/TOB)
1277  } // end loop over hits
1278  NMatchedHits_ = myHits;
1279 
1280  numSiMatchedHits_->Fill(NMatchedHits_);
1281 
1283 
1286  LogDebug("") << "Writing to myTree with " << NShowers_ << " Showers " << NStereoHits_ << " Si StereoHits "
1287  << NMonoHits_ << " Si MonoHits " << NMatchedHits_ << " Si MatchedHits ";
1288 
1289  myTree_->Fill();
1290 
1291 } // end of Analyzer
1292 
1293 void SiStripElectronAnalyzer::endJob() {
1294  LogDebug("") << "Entering endJob ";
1295  file_->cd();
1296  numCand_->Write();
1297  numElectrons_->Write();
1298  numSuperClusters_->Write();
1299 
1300  energySuperClusters_->Write();
1301  sizeSuperClusters_->Write();
1302  emaxSuperClusters_->Write();
1303  phiWidthSuperClusters_->Write();
1304 
1305  energySuperClustersEl_->Write();
1306  sizeSuperClustersEl_->Write();
1307  emaxSuperClustersEl_->Write();
1308  phiWidthSuperClustersEl_->Write();
1309 
1310  ptDiff->Write();
1311  pDiff->Write();
1312  pElectronFailed->Write();
1313  ptElectronFailed->Write();
1314  pElectronPassed->Write();
1315  ptElectronPassed->Write();
1316  sizeSuperClustersPassed->Write();
1317  sizeSuperClustersFailed->Write();
1318  // energySuperClustersPassed->Write();
1319  // energySuperClustersFailed->Write();
1320  // eOverPPassed->Write();
1321  // eOverPFailed->Write();
1322 
1323  numSiStereoHits_->Write();
1324  numSiMonoHits_->Write();
1325  numSiMatchedHits_->Write();
1326 
1327  // disable for large dataset
1328  LogDebug("") << " Writing out ntuple is disabled for now ";
1329  myTree_->Write();
1330 
1331  file_->Close();
1332 }
1333 
1334 //
1335 // const member functions
1336 //
1337 
1338 //
1339 // static member functions
1340 //
myMaxHits
#define myMaxHits
Definition: SiStripElectronAnalyzer.cc:55
LocalError::xx
float xx() const
Definition: LocalError.h:22
edm::Ref::isNonnull
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
mps_fire.i
i
Definition: mps_fire.py:428
PV3DBase::perp
T perp() const
Definition: PV3DBase.h:69
gpuClustering::id
uint16_t *__restrict__ id
Definition: gpuClusterChargeCut.h:20
SiStripElectronAnalyzer::unwrapPhi
double unwrapPhi(double phi) const
Definition: SiStripElectronAnalyzer.cc:68
mag
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
Definition: Basic3DVectorLD.h:127
DEFINE_FWK_MODULE
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
SiStripElectronAnalyzer::SiStripElectronAnalyzer
SiStripElectronAnalyzer(const edm::ParameterSet &)
Definition: SiStripElectronAnalyzer.cc:206
init
int init
Definition: HydjetWrapper.h:64
PV3DBase::phi
Geom::Phi< T > phi() const
Definition: PV3DBase.h:66
PV3DBase::y
T y() const
Definition: PV3DBase.h:60
HI_PhotonSkim_cff.electrons
tuple electrons
Definition: HI_PhotonSkim_cff.py:77
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
PV3DBase::theta
Geom::Theta< T > theta() const
Definition: PV3DBase.h:72
LocalError::xy
float xy() const
Definition: LocalError.h:23
iEvent
int iEvent
Definition: GenABIO.cc:224
LocalError::yy
float yy() const
Definition: LocalError.h:24
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
PV3DBase::z
T z() const
Definition: PV3DBase.h:61
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
SiStripRecHit2D::cluster
ClusterRef cluster() const
Definition: SiStripRecHit2D.h:22
StripSubdetector::TOB
static constexpr auto TOB
Definition: StripSubdetector.h:18
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
M_PI
#define M_PI
Definition: BXVectorInputProducer.cc:49
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
edmNew::DetSetVector::const_iterator
boost::transform_iterator< IterHelp, const_IdIter > const_iterator
Definition: DetSetVectorNew.h:197
BaseTrackerRecHit::localPositionError
LocalError localPositionError() const override
Definition: BaseTrackerRecHit.h:61
DetId
Definition: DetId.h:17
StripSubdetector::TIB
static constexpr auto TIB
Definition: StripSubdetector.h:16
math::XYZPoint
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
counter
static std::atomic< unsigned int > counter
Definition: SharedResourceNames.cc:18
position
static int position[264][3]
Definition: ReadPGInfo.cc:289
edm::EventSetup::get
T get() const
Definition: EventSetup.h:88
SiStripElectronAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: SiStripElectronAnalyzer.cc:435
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
str
#define str(s)
Definition: TestProcessor.cc:53
PV3DBase::x
T x() const
Definition: PV3DBase.h:59
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:233
Generated for CMSSW Reference Manual by doxygen 1.8.5