CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
PhotonsWithConversionsAnalyzer.cc
Go to the documentation of this file.
1 #include "CommonTools/UtilAlgos/interface/TFileService.h"
2 #include "DataFormats/Common/interface/Handle.h"
3 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
4 #include "DataFormats/EgammaCandidates/interface/ConversionFwd.h"
5 #include "DataFormats/EgammaCandidates/interface/Photon.h"
6 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
7 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
8 #include "DataFormats/TrackReco/interface/Track.h"
9 #include "DataFormats/TrackReco/interface/TrackExtra.h"
10 #include "DataFormats/TrackReco/interface/TrackFwd.h"
11 #include "FWCore/Framework/interface/ESHandle.h"
12 #include "FWCore/Framework/interface/Event.h"
13 #include "FWCore/Framework/interface/EventSetup.h"
14 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
15 #include "FWCore/MessageLogger/interface/MessageLogger.h"
16 #include "FWCore/ServiceRegistry/interface/Service.h"
17 #include "FWCore/Utilities/interface/Exception.h"
18 #include "MagneticField/Engine/interface/MagneticField.h"
19 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
20 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruth.h"
21 #include "RecoEgamma/EgammaMCTools/interface/PhotonMCTruthFinder.h"
22 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
23 #include "SimDataFormats/Track/interface/SimTrack.h"
24 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
25 #include "SimDataFormats/Vertex/interface/SimVertex.h"
26 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
27 
28 #include "TFile.h"
29 #include "TH1.h"
30 #include "TH2.h"
31 #include "TTree.h"
32 #include "TVector3.h"
33 #include "TProfile.h"
34 
35 #include <iostream>
36 #include <map>
37 #include <vector>
38 
39 class PhotonsWithConversionsAnalyzer : public edm::one::EDAnalyzer<> {
40 public:
41  //
42  explicit PhotonsWithConversionsAnalyzer(const edm::ParameterSet&);
43  ~PhotonsWithConversionsAnalyzer() override;
44 
45  void analyze(const edm::Event&, const edm::EventSetup&) override;
46  void beginJob() override;
47  void endJob() override;
48 
49 private:
50  float etaTransformation(float a, float b);
51 
52  //
53  PhotonMCTruthFinder* thePhotonMCTruthFinder_;
54 
55  std::string fOutputFileName_;
56  TFile* fOutputFile_;
57 
58  int nEvt_;
59  int nMCPho_;
60  int nMatched_;
61 
62  std::string HepMCLabel;
63  std::string SimTkLabel;
64  std::string SimVtxLabel;
65  std::string SimHitLabel;
66 
67  std::string photonCollectionProducer_;
68  std::string photonCollection_;
69 
70  TH1F* h_ErecoEMC_;
71  TH1F* h_deltaPhi_;
72  TH1F* h_deltaEta_;
73 
75  TH1F* h_MCphoE_;
76  TH1F* h_MCphoPhi_;
77  TH1F* h_MCphoEta_;
78 
80  TH1F* h_MCConvE_;
81  TH1F* h_MCConvPt_;
82  TH1F* h_MCConvEta_;
83 
84  // SC from reco photons
85  TH1F* h_scE_;
86  TH1F* h_scEt_;
87  TH1F* h_scEta_;
88  TH1F* h_scPhi_;
89  //
90  TH1F* h_phoE_;
91  TH1F* h_phoEta_;
92  TH1F* h_phoPhi_;
93  //
94  // All tracks from reco photons
95  TH2F* h2_tk_nHitsVsR_;
96  //
97  TH2F* h2_tk_inPtVsR_;
98 };
99 
100 using namespace std;
101 
102 PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer(const edm::ParameterSet& pset) {
103  photonCollectionProducer_ = pset.getParameter<std::string>("phoProducer");
104  photonCollection_ = pset.getParameter<std::string>("photonCollection");
105  //
106 }
107 
108 PhotonsWithConversionsAnalyzer::~PhotonsWithConversionsAnalyzer() { delete thePhotonMCTruthFinder_; }
109 
110 void PhotonsWithConversionsAnalyzer::beginJob() {
111  nEvt_ = 0;
112  nMCPho_ = 0;
113  nMatched_ = 0;
114 
115  thePhotonMCTruthFinder_ = new PhotonMCTruthFinder();
116 
117  edm::Service<TFileService> fs;
118 
120  h_ErecoEMC_ = fs->make<TH1F>("deltaE", " delta(reco-mc) energy", 100, 0., 2.);
121  h_deltaPhi_ = fs->make<TH1F>("deltaPhi", " delta(reco-mc) phi", 100, -0.1, 0.1);
122  h_deltaEta_ = fs->make<TH1F>("deltaEta", " delta(reco-mc) eta", 100, -0.05, 0.05);
123 
125  h_MCphoE_ = fs->make<TH1F>("MCphoE", "MC photon energy", 100, 0., 100.);
126  h_MCphoPhi_ = fs->make<TH1F>("MCphoPhi", "MC photon phi", 40, -3.14, 3.14);
127  h_MCphoEta_ = fs->make<TH1F>("MCphoEta", "MC photon eta", 40, -3., 3.);
128 
130  h_MCConvE_ = fs->make<TH1F>("MCConvE", "MC converted photon energy", 100, 0., 100.);
131  h_MCConvPt_ = fs->make<TH1F>("MCConvPt", "MC converted photon pt", 100, 0., 100.);
132  h_MCConvEta_ = fs->make<TH1F>("MCConvEta", "MC converted photon eta", 50, 0., 2.5);
133 
135  h_scE_ = fs->make<TH1F>("scE", "SC Energy ", 100, 0., 200.);
136  h_scEt_ = fs->make<TH1F>("scEt", "SC Et ", 100, 0., 200.);
137  h_scEta_ = fs->make<TH1F>("scEta", "SC Eta ", 40, -3., 3.);
138  h_scPhi_ = fs->make<TH1F>("scPhi", "SC Phi ", 40, -3.14, 3.14);
139  //
140  h_phoE_ = fs->make<TH1F>("phoE", "Photon Energy ", 100, 0., 200.);
141  h_phoEta_ = fs->make<TH1F>("phoEta", "Photon Eta ", 40, -3., 3.);
142  h_phoPhi_ = fs->make<TH1F>("phoPhi", "Photon Phi ", 40, -3.14, 3.14);
143 
144  // Recontructed tracks from converted photon candidates
145  h2_tk_nHitsVsR_ = fs->make<TH2F>("tknHitsVsR", "Tracks Hits vs R ", 12, 0., 120., 20, 0.5, 20.5);
146  h2_tk_inPtVsR_ = fs->make<TH2F>("tkInPtvsR", "Tracks inner Pt vs R ", 12, 0., 120., 100, 0., 100.);
147 
148  return;
149 }
150 
151 float PhotonsWithConversionsAnalyzer::etaTransformation(float EtaParticle, float Zvertex) {
152  //---Definitions
153  const float PI = 3.1415927;
154  //UNUSED const float TWOPI = 2.0*PI;
155 
156  //---Definitions for ECAL
157  const float R_ECAL = 136.5;
158  const float Z_Endcap = 328.0;
159  const float etaBarrelEndcap = 1.479;
160 
161  //---ETA correction
162 
163  float Theta = 0.0;
164  float ZEcal = R_ECAL * sinh(EtaParticle) + Zvertex;
165 
166  if (ZEcal != 0.0)
167  Theta = atan(R_ECAL / ZEcal);
168  if (Theta < 0.0)
169  Theta = Theta + PI;
170  float ETA = -log(tan(0.5 * Theta));
171 
172  if (fabs(ETA) > etaBarrelEndcap) {
173  float Zend = Z_Endcap;
174  if (EtaParticle < 0.0)
175  Zend = -Zend;
176  float Zlen = Zend - Zvertex;
177  float RR = Zlen / sinh(EtaParticle);
178  Theta = atan(RR / Zend);
179  if (Theta < 0.0)
180  Theta = Theta + PI;
181  ETA = -log(tan(0.5 * Theta));
182  }
183  //---Return the result
184  return ETA;
185  //---end
186 }
187 
188 void PhotonsWithConversionsAnalyzer::analyze(const edm::Event& e, const edm::EventSetup&) {
189  using namespace edm;
190  const float etaPhiDistance = 0.01;
191  // Fiducial region
192  //UNUSED const float TRK_BARL =0.9;
193  //UNUSED const float BARL = 1.4442; // DAQ TDR p.290
194  //UNUSED const float END_LO = 1.566;
195  //UNUSED const float END_HI = 2.5;
196  // Electron mass
197  //UNUSED const Float_t mElec= 0.000511;
198 
199  nEvt_++;
200  LogInfo("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id()
201  << " Global Counter " << nEvt_ << "\n";
202  // LogDebug("ConvertedPhotonAnalyzer") << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ <<"\n";
203  std::cout << "ConvertedPhotonAnalyzer Analyzing event number: " << e.id() << " Global Counter " << nEvt_ << "\n";
204 
206  Handle<reco::PhotonCollection> photonHandle;
207  e.getByLabel(photonCollectionProducer_, photonCollection_, photonHandle);
208  const reco::PhotonCollection photonCollection = *(photonHandle.product());
209  std::cout << "ConvertedPhotonAnalyzer Photons with conversions collection size " << photonCollection.size() << "\n";
210 
212  std::cout << " ConvertedPhotonAnalyzer Looking for MC truth "
213  << "\n";
214 
215  //get simtrack info
216  std::vector<SimTrack> theSimTracks;
217  std::vector<SimVertex> theSimVertices;
218 
219  edm::Handle<SimTrackContainer> SimTk;
220  edm::Handle<SimVertexContainer> SimVtx;
221  e.getByLabel("g4SimHits", SimTk);
222  e.getByLabel("g4SimHits", SimVtx);
223 
224  theSimTracks.insert(theSimTracks.end(), SimTk->begin(), SimTk->end());
225  theSimVertices.insert(theSimVertices.end(), SimVtx->begin(), SimVtx->end());
226 
227  std::vector<PhotonMCTruth> mcPhotons = thePhotonMCTruthFinder_->find(theSimTracks, theSimVertices);
228  std::cout << " ConvertedPhotonAnalyzer mcPhotons size " << mcPhotons.size() << std::endl;
229 
230  // Loop over simulated photons
231  //UNUSED int iDet=0;
232  //UNUSED int iRadius=-1;
233  //UNUSED int indPho=0;
234 
235  for (std::vector<PhotonMCTruth>::const_iterator mcPho = mcPhotons.begin(); mcPho != mcPhotons.end(); mcPho++) {
236  float mcPhi = (*mcPho).fourMomentum().phi();
237  float mcEta = (*mcPho).fourMomentum().pseudoRapidity();
238  mcEta = etaTransformation(mcEta, (*mcPho).primaryVertex().z());
239 
240  if ((*mcPho).fourMomentum().et() < 20)
241  continue;
242  // if ( ! ( fabs(mcEta) <= BARL || ( fabs(mcEta) >= END_LO && fabs(mcEta) <=END_HI ) ) ) {
243  // continue;
244  //} // all ecal fiducial region
245 
246  h_MCphoE_->Fill((*mcPho).fourMomentum().e());
247  h_MCphoEta_->Fill((*mcPho).fourMomentum().eta());
248  h_MCphoPhi_->Fill((*mcPho).fourMomentum().phi());
249 
250  // keep only visible conversions
251  if ((*mcPho).isAConversion() == 0)
252  continue;
253 
254  nMCPho_++;
255 
256  h_MCConvEta_->Fill(fabs((*mcPho).fourMomentum().pseudoRapidity()) - 0.001);
257 
258  bool REJECTED;
259 
261  //std::cout << " ConvertedPhotonAnalyzer Starting loop over photon candidates " << "\n";
262  for (reco::PhotonCollection::const_iterator iPho = photonCollection.begin(); iPho != photonCollection.end();
263  iPho++) {
264  REJECTED = false;
265 
266  // std::cout << " ConvertedPhotonAnalyzer Reco SC energy " << (*iPho).superCluster()->energy() << "\n";
267 
268  float phiClu = (*iPho).superCluster()->phi();
269  float etaClu = (*iPho).superCluster()->eta();
270  float deltaPhi = phiClu - mcPhi;
271  float deltaEta = etaClu - mcEta;
272 
273  if (deltaPhi > Geom::pi())
274  deltaPhi -= Geom::twoPi();
275  if (deltaPhi < -Geom::pi())
276  deltaPhi += Geom::twoPi();
277  deltaPhi = pow(deltaPhi, 2);
278  deltaEta = pow(deltaEta, 2);
279  float delta = deltaPhi + deltaEta;
280  if (delta >= etaPhiDistance)
281  REJECTED = true;
282 
283  // if ( ! ( fabs(etaClu) <= BARL || ( fabs(etaClu) >= END_LO && fabs(etaClu) <=END_HI ) ) ) REJECTED=true;
284 
285  if (REJECTED)
286  continue;
287  std::cout << " MATCHED " << std::endl;
288  nMatched_++;
289 
290  // std::cout << " ConvertedPhotonAnalyzer Matching candidate " << std::endl;
291 
292  // std::cout << " ConvertedPhotonAnalyzer Photons isAconversion " << (*mcPho).isAConversion() << " mcMatchingPhoton energy " << (*mcPho).fourMomentum().e() << " ConvertedPhotonAnalyzer conversion vertex R " << (*mcPho).vertex().perp() << " Z " << (*mcPho).vertex().z() << std::endl;
293 
294  h_ErecoEMC_->Fill((*iPho).superCluster()->energy() / (*mcPho).fourMomentum().e());
295  h_deltaPhi_->Fill((*iPho).superCluster()->position().phi() - mcPhi);
296  h_deltaEta_->Fill((*iPho).superCluster()->position().eta() - mcEta);
297 
298  h_scE_->Fill((*iPho).superCluster()->energy());
299  h_scEt_->Fill((*iPho).superCluster()->energy() / cosh((*iPho).superCluster()->position().eta()));
300  h_scEta_->Fill((*iPho).superCluster()->position().eta());
301  h_scPhi_->Fill((*iPho).superCluster()->position().phi());
302 
303  h_phoE_->Fill((*iPho).energy());
304  h_phoEta_->Fill((*iPho).eta());
305  h_phoPhi_->Fill((*iPho).phi());
306 
307  if (!(*iPho).hasConversionTracks())
308  continue;
309  // std::cout << " This photons has " << (*iPho).conversions().size() << " conversions candidates " << std::endl;
310  reco::ConversionRefVector conversions = (*iPho).conversions();
311  //std::vector<reco::ConversionRef> conversions = (*iPho).conversions();
312 
313  for (unsigned int i = 0; i < conversions.size(); i++) {
314  //std::cout << " Conversion candidate Energy " << (*iPho).energy() << " number of tracks " << conversions[i]->nTracks() << std::endl;
315  std::vector<edm::RefToBase<reco::Track> > tracks = conversions[i]->tracks();
316 
317  for (unsigned int i = 0; i < tracks.size(); i++) {
318  // std::cout << " ConvertedPhotonAnalyzer Reco Track charge " << tracks[i]->charge() << " Num of RecHits " << tracks[i]->recHitsSize() << " inner momentum " << sqrt ( tracks[i]->innerMomentum().Mag2() ) << "\n";
319 
320  h2_tk_nHitsVsR_->Fill((*mcPho).vertex().perp(), tracks[i]->recHitsSize());
321  h2_tk_inPtVsR_->Fill((*mcPho).vertex().perp(), sqrt(tracks[i]->innerMomentum().Mag2()));
322  }
323  }
324 
325  }
326 
327  }
328 }
329 
330 void PhotonsWithConversionsAnalyzer::endJob() {
331  // fOutputFile_->Write() ;
332  // fOutputFile_->Close() ;
333 
334  edm::LogInfo("ConvertedPhotonAnalyzer") << "Analyzed " << nEvt_ << "\n";
335  // std::cout << "::endJob Analyzed " << nEvt_ << " events " << " with total " << nPho_ << " Photons " << "\n";
336  std::cout << "ConvertedPhotonAnalyzer::endJob Analyzed " << nEvt_ << " events "
337  << "\n";
338 
339  return;
340 }
log
static std::vector< std::string > checklist log
Definition: LHERunInfoProduct.cc:179
mps_fire.i
i
Definition: mps_fire.py:428
TrackValidation_cff.pset
tuple pset
Definition: TrackValidation_cff.py:591
TFileService::make
T * make(const Args &...args) const
make new ROOT object
Definition: TFileService.h:64
tracks
auto const & tracks
cannot be loose
Definition: CAHitNtupletGeneratorKernelsImpl.h:146
R_ECAL
static constexpr float R_ECAL
Definition: ECALPositionCalculator.cc:10
srCondWrite_cfg.deltaPhi
tuple deltaPhi
Definition: srCondWrite_cfg.py:54
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
spr::deltaEta
static const double deltaEta
Definition: CaloConstants.h:8
PhotonsWithConversionsAnalyzer::PhotonsWithConversionsAnalyzer
PhotonsWithConversionsAnalyzer(const edm::ParameterSet &)
Definition: PhotonsWithConversionsAnalyzer.cc:102
ETA
#define ETA
Definition: GenericBenchmark.cc:28
conversions_cfi.conversions
tuple conversions
Definition: conversions_cfi.py:8
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
funct::tan
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
etaBarrelEndcap
static constexpr float etaBarrelEndcap
Definition: ECALPositionCalculator.cc:12
PI
#define PI
Definition: QcdUeDQM.h:37
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:500
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:125
reco::PhotonCollection
std::vector< Photon > PhotonCollection
collectin of Photon objects
Definition: PhotonFwd.h:9
edm::ParameterSet::getParameter
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
b
double b
Definition: hdecay.h:118
PhotonsWithConversionsAnalyzer::analyze
void analyze(const edm::Event &, const edm::EventSetup &) override
Definition: PhotonsWithConversionsAnalyzer.cc:188
edm::EventBase::id
edm::EventID id() const
Definition: EventBase.h:59
a
double a
Definition: hdecay.h:119
ZEcal
static constexpr float ZEcal
Definition: L1TkEmParticleProducer.cc:37
edm::RefVector::size
size_type size() const
Size of the RefVector.
Definition: RefVector.h:102
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:144
Geom::pi
constexpr double pi()
Definition: Pi.h:31
Geom::twoPi
constexpr double twoPi()
Definition: Pi.h:32
funct::pow
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
Z_Endcap
static constexpr float Z_Endcap
Definition: ECALPositionCalculator.cc:11
Generated for CMSSW Reference Manual by doxygen 1.8.5