CMS 3D CMS Logo

HIPTwoBodyDecayAnalyzer.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <memory>
4 #include <utility>
5 #include <cmath>
6 #include <string>
7 #include "TLorentzVector.h"
8 #include "TTree.h"
9 #include "TH1F.h"
10 #include "TMath.h"
11 
12 // user include files
15 
24 
31 
32 class HIPTwoBodyDecayAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
33 public:
35  ~HIPTwoBodyDecayAnalyzer() override;
36 
41 
42  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
43 
44  TTree* tree;
45  std::vector<std::pair<std::string, float*>> floatBranches;
46  std::vector<std::pair<std::string, int*>> intBranches;
47  std::vector<std::pair<std::string, short*>> shortBranches;
48 
49 private:
50  // es consumes
52 
54 
55  void beginJob() override;
56  void analyze(const edm::Event&, const edm::EventSetup&) override;
57  void endJob() override;
58 
59  void bookAllBranches();
61  BranchType searchArray(std::string branchname, int& position);
62  template <typename varType>
64  int varposition = -1;
65  BranchType varbranchtype = searchArray(bname, varposition);
66  if (varposition == -1)
67  std::cerr << "HIPTwoBodyDecayAnalyzer::setVal -> Could not find the branch called " << bname << "!" << std::endl;
68  else if (varbranchtype == BranchType_short_t)
69  *(shortBranches.at(varposition).second) = value;
70  else if (varbranchtype == BranchType_int_t)
71  *(intBranches.at(varposition).second) = value;
72  else if (varbranchtype == BranchType_float_t)
73  *(floatBranches.at(varposition).second) = value;
74  else
75  std::cerr << "HIPTwoBodyDecayAnalyzer::setVal -> Could not find the type " << varbranchtype
76  << " for branch called " << bname << "!" << std::endl;
77  }
78  void cleanBranches();
79  void initializeBranches();
80  bool actuateBranches();
81 
82  void analyzeTrackCollection(std::string strTrackType,
83  const TransientTrackBuilder& theTTBuilder,
85  bool verbose = false);
88  bool& fitOk);
89 };
90 
92  : ttkbuilderToken_(esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))) {
93  usesResource(TFileService::kSharedResource);
94 
95  alcareco_trackCollToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("alcarecotracks"));
96  refit1_trackCollToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("refit1tracks"));
97  ctf_trackCollToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("refit2tracks"));
98  final_trackCollToken_ = consumes<reco::TrackCollection>(iConfig.getParameter<edm::InputTag>("finaltracks"));
99 
101 
102  tree = fs->make<TTree>("TestTree", "");
103  bookAllBranches();
104 }
105 
107  for (unsigned short el = 0; el < shortBranches.size(); el++) {
108  if (branchname == shortBranches.at(el).first) {
109  position = el;
110  return BranchType_short_t;
111  }
112  }
113  for (unsigned int el = 0; el < intBranches.size(); el++) {
114  if (branchname == intBranches.at(el).first) {
115  position = el;
116  return BranchType_int_t;
117  }
118  }
119  for (unsigned int el = 0; el < floatBranches.size(); el++) {
120  if (branchname == floatBranches.at(el).first) {
121  position = el;
122  return BranchType_float_t;
123  }
124  }
125  return BranchType_unknown_t;
126 }
128  for (unsigned short el = 0; el < shortBranches.size(); el++) {
129  if (shortBranches.at(el).second != nullptr)
130  delete shortBranches.at(el).second;
131  shortBranches.at(el).second = nullptr;
132  }
133  shortBranches.clear();
134  for (unsigned int el = 0; el < intBranches.size(); el++) {
135  if (intBranches.at(el).second != nullptr)
136  delete intBranches.at(el).second;
137  intBranches.at(el).second = nullptr;
138  }
139  intBranches.clear();
140  for (unsigned int el = 0; el < floatBranches.size(); el++) {
141  if (floatBranches.at(el).second != nullptr)
142  delete floatBranches.at(el).second;
143  floatBranches.at(el).second = nullptr;
144  }
145  floatBranches.clear();
146 }
148  for (unsigned short el = 0; el < shortBranches.size(); el++) {
149  if (shortBranches.at(el).second != nullptr)
150  *(shortBranches.at(el).second) = 0;
151  }
152  for (unsigned int el = 0; el < intBranches.size(); el++) {
153  if (intBranches.at(el).second != nullptr)
154  *(intBranches.at(el).second) = 0;
155  }
156  for (unsigned int el = 0; el < floatBranches.size(); el++) {
157  if (floatBranches.at(el).second != nullptr)
158  *(floatBranches.at(el).second) = 0;
159  }
160 }
161 
163  const int nTrackTypes = 4;
164  std::vector<std::string> strTrackTypes;
165  strTrackTypes.reserve(nTrackTypes);
166  strTrackTypes.push_back("alcareco");
167  strTrackTypes.push_back("refit1");
168  strTrackTypes.push_back("refit2");
169  strTrackTypes.push_back("final");
170  for (unsigned int it = 0; it < nTrackTypes; it++) {
171  std::string& strTrackType = strTrackTypes[it];
172  bookBranch(strTrackType + "_present", BranchType_short_t);
173  bookBranch(strTrackType + "_ZVtxFitOk", BranchType_short_t);
174  bookBranch(strTrackType + "_ZMass", BranchType_float_t);
175  bookBranch(strTrackType + "_ZPt", BranchType_float_t);
176  bookBranch(strTrackType + "_ZPz", BranchType_float_t);
177  bookBranch(strTrackType + "_ZPhi", BranchType_float_t);
178  bookBranch(strTrackType + "_ZVertex_x", BranchType_float_t);
179  bookBranch(strTrackType + "_ZVertex_y", BranchType_float_t);
180  bookBranch(strTrackType + "_ZVertex_z", BranchType_float_t);
181  bookBranch(strTrackType + "_ZVertex_NormChi2", BranchType_float_t);
182  bookBranch(strTrackType + "_MuPlusVertex_x", BranchType_float_t);
183  bookBranch(strTrackType + "_MuPlusVertex_y", BranchType_float_t);
184  bookBranch(strTrackType + "_MuPlusVertex_z", BranchType_float_t);
185  bookBranch(strTrackType + "_MuMinusPt_AfterZVtxFit", BranchType_float_t);
186  bookBranch(strTrackType + "_MuMinusPz_AfterZVtxFit", BranchType_float_t);
187  bookBranch(strTrackType + "_MuMinusPhi_AfterZVtxFit", BranchType_float_t);
188  bookBranch(strTrackType + "_MuPlusPt_AfterZVtxFit", BranchType_float_t);
189  bookBranch(strTrackType + "_MuPlusPz_AfterZVtxFit", BranchType_float_t);
190  bookBranch(strTrackType + "_MuPlusPhi_AfterZVtxFit", BranchType_float_t);
191  bookBranch(strTrackType + "_ZMass_AfterZVtxFit", BranchType_float_t);
192  bookBranch(strTrackType + "_ZPt_AfterZVtxFit", BranchType_float_t);
193  bookBranch(strTrackType + "_ZPz_AfterZVtxFit", BranchType_float_t);
194  bookBranch(strTrackType + "_ZPhi_AfterZVtxFit", BranchType_float_t);
195  bookBranch(strTrackType + "_MuMinusPt", BranchType_float_t);
196  bookBranch(strTrackType + "_MuMinusPz", BranchType_float_t);
197  bookBranch(strTrackType + "_MuMinusPhi", BranchType_float_t);
198  bookBranch(strTrackType + "_MuMinusVertex_x", BranchType_float_t);
199  bookBranch(strTrackType + "_MuMinusVertex_y", BranchType_float_t);
200  bookBranch(strTrackType + "_MuMinusVertex_z", BranchType_float_t);
201  bookBranch(strTrackType + "_MuPlusPt", BranchType_float_t);
202  bookBranch(strTrackType + "_MuPlusPz", BranchType_float_t);
203  bookBranch(strTrackType + "_MuPlusPhi", BranchType_float_t);
204  }
205  actuateBranches();
206 }
208  if (btype == BranchType_float_t)
209  floatBranches.emplace_back(bname, new float);
210  else if (btype == BranchType_int_t)
211  intBranches.emplace_back(bname, new int);
212  else if (btype == BranchType_short_t)
213  shortBranches.emplace_back(bname, new short);
214  else {
215  std::cerr << "HIPTwoBodyDecayAnalyzer::bookBranch: No support for type " << btype << " for the branch " << bname
216  << " is available." << std::endl;
217  return false;
218  }
219  return true;
220 }
222  bool success = true;
223  std::cout << "Begin HIPTwoBodyDecayAnalyzer::actuateBranches" << std::endl;
224  std::cout << "Number of short branches: " << shortBranches.size() << std::endl;
225  std::cout << "Number of int branches: " << intBranches.size() << std::endl;
226  std::cout << "Number of float branches: " << floatBranches.size() << std::endl;
227  if (tree != nullptr) {
228  for (unsigned short el = 0; el < shortBranches.size(); el++) {
229  std::cout << "Actuating branch " << shortBranches.at(el).first << " at address " << shortBranches.at(el).second
230  << std::endl;
231  if (!tree->GetBranchStatus(shortBranches.at(el).first.c_str()))
232  tree->Branch(shortBranches.at(el).first.c_str(), shortBranches.at(el).second);
233  else
234  std::cout << "Failed!" << std::endl;
235  }
236  for (unsigned int el = 0; el < intBranches.size(); el++) {
237  std::cout << "Actuating branch " << intBranches.at(el).first.c_str() << " at address "
238  << intBranches.at(el).second << std::endl;
239  if (!tree->GetBranchStatus(intBranches.at(el).first.c_str()))
240  tree->Branch(intBranches.at(el).first.c_str(), intBranches.at(el).second);
241  else
242  std::cout << "Failed!" << std::endl;
243  }
244  for (unsigned int el = 0; el < floatBranches.size(); el++) {
245  std::cout << "Actuating branch " << floatBranches.at(el).first.c_str() << " at address "
246  << floatBranches.at(el).second << std::endl;
247  if (!tree->GetBranchStatus(floatBranches.at(el).first.c_str()))
248  tree->Branch(floatBranches.at(el).first.c_str(), floatBranches.at(el).second);
249  else
250  std::cout << "Failed!" << std::endl;
251  }
252  } else
253  success = false;
254  if (!success)
255  std::cerr << "HIPTwoBodyDecayAnalyzer::actuateBranch: Failed to actuate the branches!" << std::endl;
256  return success;
257 }
259 
261  using namespace edm;
262  using namespace reco;
263  using reco::TrackCollection;
264 
265  edm::Handle<reco::TrackCollection> alcarecotracks;
266  iEvent.getByToken(alcareco_trackCollToken_, alcarecotracks);
268  iEvent.getByToken(refit1_trackCollToken_, refit1tracks);
270  iEvent.getByToken(ctf_trackCollToken_, ctftracks);
272  iEvent.getByToken(final_trackCollToken_, finaltracks);
273 
274  const auto& theTTBuilder = iSetup.getData(ttkbuilderToken_);
276 
277  analyzeTrackCollection("alcareco", theTTBuilder, alcarecotracks);
278  analyzeTrackCollection("refit1", theTTBuilder, refit1tracks);
279  analyzeTrackCollection("refit2", theTTBuilder, ctftracks);
280  analyzeTrackCollection("final", theTTBuilder, finaltracks);
281 
282  tree->Fill();
283 }
284 
287 
290  desc.setUnknown();
291  descriptions.addDefault(desc);
292 }
293 
295  const TransientTrackBuilder& theTTBuilder,
297  bool verbose) {
298  if (verbose)
299  std::cout << "Starting to process the track collection for " << strTrackType << std::endl;
300 
301  using namespace edm;
302  using namespace reco;
303  using reco::TrackCollection;
304 
305  if (!hTrackColl.isValid()) {
306  if (verbose)
307  std::cout << "Track collection is invalid." << std::endl;
308  return;
309  }
310  if (hTrackColl->size() < 2) {
311  if (verbose)
312  std::cout << "Track collection size<2." << std::endl;
313  return;
314  }
315 
316  unsigned int itrk = 0;
317  unsigned int j = 0;
318  int totalcharge = 0;
319  bool isValidPair = true;
320  bool ZVtxOk = false;
321  TLorentzVector trackMom[2];
322  TLorentzVector trackMomAfterZVtxFit[2];
323  TVector3 trackVtx[2];
324 
325  for (unsigned int jtrk = 0; jtrk < 2; jtrk++) {
326  trackMom[jtrk].SetXYZT(0, 0, 0, 0);
327  trackVtx[jtrk].SetXYZ(0, 0, 0);
328  }
329  for (reco::TrackCollection::const_iterator track = hTrackColl->begin(); track != hTrackColl->end(); ++track) {
330  int charge = track->charge();
331  totalcharge += charge;
332  if (j == 0) {
333  itrk = (charge > 0 ? 1 : 0);
334  } else
335  itrk = 1 - itrk;
336  trackMom[itrk].SetPtEtaPhiM(track->pt(), track->eta(), track->phi(), 0.105);
337  trackVtx[itrk].SetXYZ(track->vx(), track->vy(), track->vz());
338  j++;
339  if (j == 2)
340  break;
341  }
342 
343  isValidPair = (totalcharge == 0 && trackMom[0].P() != 0. && trackMom[1].P() != 0.);
344  if (verbose && !isValidPair)
345  std::cout << "Track collection does not contain a valid std::pair." << std::endl;
346  setVal(strTrackType + "_present", (isValidPair ? (short)1 : (short)0));
347  if (isValidPair) {
348  TLorentzVector ZMom = trackMom[0] + trackMom[1];
349  setVal(strTrackType + "_ZPt", (float)ZMom.Pt());
350  setVal(strTrackType + "_ZPz", (float)ZMom.Pz());
351  setVal(strTrackType + "_ZPhi", (float)ZMom.Phi());
352  setVal(strTrackType + "_ZMass", (float)ZMom.M());
353 
354  reco::Vertex ZVtx = fitDimuonVertex(theTTBuilder, hTrackColl, ZVtxOk);
355  if (ZVtxOk) {
356  setVal(strTrackType + "_ZVertex_x", (float)ZVtx.x());
357  setVal(strTrackType + "_ZVertex_y", (float)ZVtx.y());
358  setVal(strTrackType + "_ZVertex_z", (float)ZVtx.z());
359  setVal(strTrackType + "_ZVertex_NormChi2", (float)ZVtx.normalizedChi2());
360 
361  // Recalculate track momenta with this vertex as reference
362  j = 0;
363  for (reco::TrackCollection::const_iterator track = hTrackColl->begin(); track != hTrackColl->end(); ++track) {
364  TransientTrack t_track = theTTBuilder.build(&(*track));
365  AnalyticalImpactPointExtrapolator extrapolator(t_track.field());
366  TrajectoryStateOnSurface closestIn3DSpaceState =
367  extrapolator.extrapolate(t_track.impactPointState(), RecoVertex::convertPos(ZVtx.position()));
368  GlobalVector mom = closestIn3DSpaceState.globalMomentum();
369  int charge = track->charge();
370  totalcharge += charge;
371  if (j == 0) {
372  itrk = (charge > 0 ? 1 : 0);
373  } else
374  itrk = 1 - itrk;
375  trackMomAfterZVtxFit[itrk].SetXYZT(mom.x(), mom.y(), mom.z(), sqrt(pow(0.105, 2) + pow(mom.mag(), 2)));
376  j++;
377  if (j == 2)
378  break;
379  }
380  if (totalcharge != 0)
381  std::cerr
382  << "HIPTwoBodyDecayAnalyzer::analyzeTrackCollection: Something went wrong! The total charge is no longer 0!"
383  << std::endl;
384  for (unsigned int jtrk = 0; jtrk < 2; jtrk++) {
385  std::string strMuCore = (jtrk == 0 ? "MuMinus" : "MuPlus");
386  setVal(strTrackType + "_" + strMuCore + "Pt_AfterZVtxFit", (float)trackMomAfterZVtxFit[jtrk].Pt());
387  setVal(strTrackType + "_" + strMuCore + "Pz_AfterZVtxFit", (float)trackMomAfterZVtxFit[jtrk].Pz());
388  setVal(strTrackType + "_" + strMuCore + "Phi_AfterZVtxFit", (float)trackMomAfterZVtxFit[jtrk].Phi());
389  }
390  TLorentzVector ZMom_AfterZVtxFit = trackMomAfterZVtxFit[0] + trackMomAfterZVtxFit[1];
391  setVal(strTrackType + "_ZPt_AfterZVtxFit", (float)ZMom_AfterZVtxFit.Pt());
392  setVal(strTrackType + "_ZPz_AfterZVtxFit", (float)ZMom_AfterZVtxFit.Pz());
393  setVal(strTrackType + "_ZPhi_AfterZVtxFit", (float)ZMom_AfterZVtxFit.Phi());
394  setVal(strTrackType + "_ZMass_AfterZVtxFit", (float)ZMom_AfterZVtxFit.M());
395  } else
396  std::cerr << "HIPTwoBodyDecayAnalyzer::analyzeTrackCollection: Z vertex fit failed for track collection "
397  << strTrackType << std::endl;
398  }
399  setVal(strTrackType + "_ZVtxFitOk", (ZVtxOk ? (short)1 : (short)0));
400  for (unsigned int jtrk = 0; jtrk < 2; jtrk++) {
401  std::string strMuCore = (jtrk == 0 ? "MuMinus" : "MuPlus");
402  setVal(strTrackType + "_" + strMuCore + "Pt", (float)trackMom[jtrk].Pt());
403  setVal(strTrackType + "_" + strMuCore + "Pz", (float)trackMom[jtrk].Pz());
404  setVal(strTrackType + "_" + strMuCore + "Phi", (float)trackMom[jtrk].Phi());
405  setVal(strTrackType + "_" + strMuCore + "Vertex_x", (float)trackVtx[jtrk].X());
406  setVal(strTrackType + "_" + strMuCore + "Vertex_y", (float)trackVtx[jtrk].Y());
407  setVal(strTrackType + "_" + strMuCore + "Vertex_z", (float)trackVtx[jtrk].Z());
408  }
409 }
410 
413  bool& fitOk) {
414  using namespace edm;
415  using namespace reco;
416 
417  std::vector<TransientTrack> t_tks;
418  for (TrackCollection::const_iterator track = hTrackColl->begin(); track != hTrackColl->end(); ++track) {
419  TransientTrack tt = theTTBuilder.build(&(*track));
420  t_tks.push_back(tt);
421  }
422 
423  // Kalman vertex fit without constraint
424  KalmanVertexFitter vtxFitter;
425  TransientVertex stdVertex = vtxFitter.vertex(t_tks);
426  fitOk = stdVertex.isValid();
427  if (fitOk) {
428  reco::Vertex stdRecoVertex(Vertex::Point(stdVertex.position()),
429  stdVertex.positionError().matrix(),
430  stdVertex.totalChiSquared(),
431  stdVertex.degreesOfFreedom(),
432  0);
433  return stdRecoVertex;
434  } else {
435  reco::Vertex stdRecoVertex;
436  return stdRecoVertex;
437  }
438 }
439 
static const std::string kSharedResource
Definition: TFileService.h:76
reco::Vertex::Point convertPos(const GlobalPoint &p)
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
float totalChiSquared() const
edm::EDGetTokenT< reco::TrackCollection > refit1_trackCollToken_
const edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord > ttkbuilderToken_
GlobalPoint position() const
T z() const
Definition: PV3DBase.h:61
GlobalError positionError() const
bool verbose
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
bool bookBranch(std::string bname, BranchType btype)
#define X(str)
Definition: MuonsGrabber.cc:38
std::vector< std::pair< std::string, short * > > shortBranches
float degreesOfFreedom() const
void analyze(const edm::Event &, const edm::EventSetup &) override
reco::Vertex fitDimuonVertex(const TransientTrackBuilder &theTTBuilder, edm::Handle< reco::TrackCollection > &hTrackColl, bool &fitOk)
void setVal(std::string bname, varType value)
CachingVertex< 5 > vertex(const std::vector< reco::TransientTrack > &tracks) const override
reco::TransientTrack build(const reco::Track *p) const
T x() const
Definition: PV3DBase.h:59
T y() const
Definition: PV3DBase.h:60
int iEvent
Definition: GenABIO.cc:224
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void addDefault(ParameterSetDescription const &psetDescription)
edm::EDGetTokenT< reco::TrackCollection > alcareco_trackCollToken_
bool isValid() const
T sqrt(T t)
Definition: SSEVec.h:19
const MagneticField * field() const
T mag() const
Definition: PV3DBase.h:64
math::XYZPoint Point
point in the space
Definition: Vertex.h:40
const AlgebraicSymMatrix33 matrix() const
std::vector< std::pair< std::string, float * > > floatBranches
bool getData(T &iHolder) const
Definition: EventSetup.h:122
Definition: value.py:1
std::vector< std::pair< std::string, int * > > intBranches
HIPTwoBodyDecayAnalyzer(const edm::ParameterSet &)
BranchType searchArray(std::string branchname, int &position)
bool isValid() const
Definition: HandleBase.h:70
fixed size matrix
HLT enums.
static int position[264][3]
Definition: ReadPGInfo.cc:289
Definition: tree.py:1
edm::EDGetTokenT< reco::TrackCollection > ctf_trackCollToken_
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
void analyzeTrackCollection(std::string strTrackType, const TransientTrackBuilder &theTTBuilder, edm::Handle< reco::TrackCollection > &hTrackColl, bool verbose=false)
TrajectoryStateOnSurface impactPointState() const
edm::EDGetTokenT< reco::TrackCollection > final_trackCollToken_