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