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