ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/codeHeavyBaryons/LambdaB/src/LambdaAnalyzer.cc
Revision: 1.2
Committed: Fri Dec 10 10:29:56 2010 UTC (14 years, 5 months ago) by mirena
Content type: text/plain
Branch: MAIN
CVS Tags: HEAD
Changes since 1.1: +1 -1 lines
State: FILE REMOVED
Log Message:
*** empty log message ***

File Contents

# Content
1 // -*- C++ -*-
2 //
3 // Package: LambdaAnalyzer
4 // Class: LambdaAnalyzer
5 //
6 /**\class LambdaAnalyzer LambdaAnalyzer.cc Lambda/LambdaAnalyzer/src/LambdaAnalyzer.cc
7
8 Description: [one line class summary]
9
10 Implementation:
11 [Notes on implementation]
12 */
13 //
14 // Original Author: Mirena Ivova Rikova,15 R-025,+41227678361,
15 // Created: Wed Aug 4 12:34:06 CEST 2010
16 // $Id: LambdaAnalyzer.cc,v 1.1 2010/11/28 17:41:45 mirena Exp $
17 //
18 //
19
20
21 // system include files
22 #include <memory>
23
24 // user include files
25 #include "FWCore/Framework/interface/Frameworkfwd.h"
26 #include "FWCore/Framework/interface/EDAnalyzer.h"
27
28 #include "FWCore/Framework/interface/Event.h"
29 #include "FWCore/Framework/interface/MakerMacros.h"
30
31 #include "FWCore/ParameterSet/interface/ParameterSet.h"
32
33 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
34
35 #include "DataFormats/TrackReco/interface/Track.h"
36 #include "DataFormats/Candidate/interface/VertexCompositeCandidateFwd.h"
37 #include "DataFormats/Candidate/interface/VertexCompositeCandidate.h"
38 #include "FWCore/ServiceRegistry/interface/Service.h"
39 #include "CommonTools/UtilAlgos/interface/TFileService.h"
40 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidate.h"
41 #include "DataFormats/RecoCandidate/interface/RecoChargedCandidateFwd.h"
42 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
43
44 #include "DataFormats/VertexReco/interface/Vertex.h"
45 #include "DataFormats/VertexReco/interface/VertexFwd.h"
46
47 #include "DataFormats/PatCandidates/interface/Muon.h"
48 #include "DataFormats/MuonReco/interface/MuonFwd.h"
49 #include "DataFormats/MuonReco/interface/Muon.h"
50
51 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
52 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
53 //#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
54 //#include "MagneticField/Engine/interface/MagneticField.h"
55
56 #include "RecoVertex/KalmanVertexFit/interface/KalmanVertexFitter.h"
57 #include "PhysicsTools/CandUtils/interface/AddFourMomenta.h"
58
59 //try trajectory backwards extrapolation
60 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
61 #include "TrackingTools/GeomPropagators/interface/Propagator.h"
62 //triggers
63 #include "DataFormats/Common/interface/TriggerResults.h"
64 #include "FWCore/Common/interface/TriggerNames.h"
65
66 //triggers
67 #include "DataFormats/HLTReco/interface/TriggerEvent.h"
68 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
69 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
70 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
71
72 //Vertex
73 #include "RecoVertex/VertexTools/interface/VertexDistanceXY.h"
74 #include "RecoVertex/VertexTools/interface/VertexDistance3D.h"
75 #include "TVector3.h"
76
77 //Kinematic vertex fit
78
79 #include "RecoVertex/KinematicFit/interface/KinematicParticleVertexFitter.h"
80 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicParticleFactoryFromTransientTrack.h"
81 #include "RecoVertex/KinematicFit/interface/MassKinematicConstraint.h"
82 #include "RecoVertex/KinematicFit/interface/KinematicParticleFitter.h"
83 #include "RecoVertex/KinematicFitPrimitives/interface/MultiTrackKinematicConstraint.h"
84 #include "RecoVertex/KinematicFit/interface/KinematicConstrainedVertexFitter.h"
85 #include "RecoVertex/KinematicFit/interface/TwoTrackMassKinematicConstraint.h"
86
87
88
89 #include "TH1F.h"
90 #include "TH2F.h"
91 #include "TTree.h"
92
93
94
95 //
96 // class declaration
97 //
98
99 class LambdaAnalyzer : public edm::EDAnalyzer {
100 public:
101 explicit LambdaAnalyzer(const edm::ParameterSet&);
102 ~LambdaAnalyzer();
103
104
105 private:
106 virtual void beginJob() ;
107 virtual void analyze(const edm::Event&, const edm::EventSetup&);
108 virtual void endJob() ;
109
110 virtual bool pass_HLT_Mu0_TkMu0_Jpsi(const edm::Event&);
111 virtual bool pass_HLT_Mu0_TkMu0_OST_Jpsi(const edm::Event&);
112 virtual bool pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v1(const edm::Event&);
113 virtual bool pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v2(const edm::Event&);
114 virtual bool pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v3(const edm::Event&);
115 virtual bool pass_HLT_DoubleMu0_Quarkonium_v1(const edm::Event&);
116
117 virtual void push_anything_back(); //fill tree with default values for _after_fit_ parameters in case of not succeessful kin fit
118
119
120
121
122 unsigned int tmpNvtx;
123
124 TTree* tree_;
125 unsigned int nLambdaB;
126
127 //reco particles info
128 float run;
129 float lumiBlock;
130 float event;
131 // float priVtxX, priVtxY, priVtxZ, priVtxXE, priVtxYE, priVtxZE, priVtxCL;
132 // float priVtxBSx, priVtxBSy, priVtxBSz, priVtxBSxE, priVtxBSyE, priVtxBSzE, priVtxBSCL;
133 float beamSpotX, beamSpotY, beamSpotZ;
134 unsigned int hlt_mu3, hlt_mu5, hlt_mu7, hlt_mu9, hlt_2mu0, hlt_2mu3, hlt_2mu0L2, hlt_2mu3JPsi, hlt_BJPsiMuMu;
135 unsigned int hlt_mu0trk0, hlt_mu3trk0, hlt_mu0trkmu0, hlt_mu3trkmu0, hlt_L1muOpen, hlt_L12muOpen, hlt_L12muOpenTight;
136 unsigned int hlt_mu0_L1muOpen, hlt_mu3_L1muOpen, hlt_mu5_L1muOpen, hlt_mu5trk0, hlt_mu0trkmu0_nocharge;
137 unsigned int hlt_mu0trkmu0_OST, hlt_2mu0_Qv1, hlt_2mu0_QLSv1, hlt_mu0trkmu0_OST_tightV1, hlt_mu0trkmu0_OST_tightV2, hlt_mu0trkmu0_OST_tightV3;
138
139 vector<float> *priVtxX;
140 vector<float> *priVtxY;
141 vector<float> *priVtxZ;
142 vector<float> *priVtxXE;
143 vector<float> *priVtxYE;
144 vector<float> *priVtxZE;
145 vector<float> *priVtxCL;
146 vector<float> *priVtxBSx;
147 vector<float> *priVtxBSy;
148 vector<float> *priVtxBSz;
149 vector<float> *priVtxBSxE;
150 vector<float> *priVtxBSyE;
151 vector<float> *priVtxBSzE;
152 vector<float> *priVtxBSCL;
153
154
155
156
157
158
159
160 //before kinematic fit
161 vector<float> *JpsiVtxX;
162 vector<float> *JpsiVtxY;
163 vector<float> *JpsiVtxZ;
164 vector<float> *LambdaVtxX;
165 vector<float> *LambdaVtxY;
166 vector<float> *LambdaVtxZ;
167 vector<float> *JpsiVtxProb;
168 vector<float> *JpsiVtxCL;
169 vector<float> *LambdaVtxProb;
170 vector<float> *LambdaVtxCL;
171 vector<float> *LambdaBMass;
172 vector<float> *LambdaBPx;
173 vector<float> *LambdaBPy;
174 vector<float> *LambdaBPz;
175 vector<float> *JpsiPx;
176 vector<float> *JpsiPy;
177 vector<float> *JpsiPz;
178 vector<float> *LambdaPx;
179 vector<float> *LambdaPy;
180 vector<float> *LambdaPz;
181 vector<float> *JpsiMass;
182 vector<float> *LambdaMass;
183 vector<float> *mupPx;
184 vector<float> *mupPy;
185 vector<float> *mupPz;
186 vector<float> *mumPx;
187 vector<float> *mumPy;
188 vector<float> *mumPz;
189 vector<float> *trackPPx;
190 vector<float> *trackPPy;
191 vector<float> *trackPPz;
192 vector<float> *trackMPx;
193 vector<float> *trackMPy;
194 vector<float> *trackMPz;
195 vector<float> *JpsiLambdaVtxDistXY;
196 vector<float> *JpsiLambdaVtxDist3D;
197 vector<float> *openingAngleLambdaB;
198 vector<float> *openingAngleLambda;
199 vector<float> *LxyPV;
200 vector<float> *LxyBS;
201 vector<float> *L3DPV;
202 vector<float> *L3DBS;
203 vector<float> *ctauLambdaBxyPV;
204 vector<float> *ctauLambdaBxyBS;
205 vector<float> *ctauLambdaB3DPV;
206 vector<float> *ctauLambdaB3DBS;
207 vector<float> *dRmup;
208 vector<float> *dRmum;
209 vector<float> *dRtrackp;
210 vector<float> *dRtrackm;
211 vector<float> *dRLambdaVtx;
212 vector<float> *mupChi2;
213 vector<float> *mumChi2;
214 vector<float> *mupNDOF;
215 vector<float> *mumNDOF;
216 vector<float> *mupD0;
217 vector<float> *mumD0;
218 vector<float> *mupNvalidHits;
219 vector<float> *mumNvalidHits;
220 vector<float> *mupIsGlobal;
221 vector<float> *mumIsGlobal;
222 vector<float> *mupIsTracker;
223 vector<float> *mumIsTracker;
224 vector<float> *ctauLambdaXY;
225
226 //after kinematic fit
227 vector<float> *JpsiVtxX_afterKinFit;
228 vector<float> *JpsiVtxY_afterKinFit;
229 vector<float> *JpsiVtxZ_afterKinFit;
230 vector<float> *LambdaVtxX_afterKinFit;
231 vector<float> *LambdaVtxY_afterKinFit;
232 vector<float> *LambdaVtxZ_afterKinFit;
233 vector<float> *JpsiVtxProb_afterKinFit;
234 vector<float> *JpsiVtxCL_afterKinFit;
235 vector<float> *LambdaVtxProb_afterKinFit;
236 vector<float> *LambdaVtxCL_afterKinFit;
237 vector<float> *LambdaBMass_afterKinFit;
238 vector<float> *LambdaBPx_afterKinFit;
239 vector<float> *LambdaBPy_afterKinFit;
240 vector<float> *LambdaBPz_afterKinFit;
241 vector<float> *JpsiPx_afterKinFit;
242 vector<float> *JpsiPy_afterKinFit;
243 vector<float> *JpsiPz_afterKinFit;
244 vector<float> *LambdaPx_afterKinFit;
245 vector<float> *LambdaPy_afterKinFit;
246 vector<float> *LambdaPz_afterKinFit;
247 vector<float> *JpsiMass_afterKinFit;
248 vector<float> *LambdaMass_afterKinFit;
249 vector<float> *mupPx_afterKinFit;
250 vector<float> *mupPy_afterKinFit;
251 vector<float> *mupPz_afterKinFit;
252 vector<float> *mumPx_afterKinFit;
253 vector<float> *mumPy_afterKinFit;
254 vector<float> *mumPz_afterKinFit;
255 vector<float> *trackPPx_afterKinFit;
256 vector<float> *trackPPy_afterKinFit;
257 vector<float> *trackPPz_afterKinFit;
258 vector<float> *trackMPx_afterKinFit;
259 vector<float> *trackMPy_afterKinFit;
260 vector<float> *trackMPz_afterKinFit;
261 vector<float> *JpsiLambdaVtxDistXY_afterKinFit;
262 vector<float> *JpsiLambdaVtxDist3D_afterKinFit;
263 vector<float> *openingAngleLambdaB_afterKinFit;
264 vector<float> *openingAngleLambda_afterKinFit;
265 vector<float> *LxyPV_afterKinFit;
266 vector<float> *LxyBS_afterKinFit;
267 vector<float> *L3DPV_afterKinFit;
268 vector<float> *L3DBS_afterKinFit;
269 vector<float> *ctauLambdaBxyPV_afterKinFit;
270 vector<float> *ctauLambdaBxyBS_afterKinFit;
271 vector<float> *ctauLambdaB3DPV_afterKinFit;
272 vector<float> *ctauLambdaB3DBS_afterKinFit;
273 vector<float> *dRmup_afterKinFit;
274 vector<float> *dRmum_afterKinFit;
275 vector<float> *dRtrackp_afterKinFit;
276 vector<float> *dRtrackm_afterKinFit;
277 vector<float> *dRLambdaVtx_afterKinFit;
278 vector<float> *mupChi2_afterKinFit;
279 vector<float> *mumChi2_afterKinFit;
280 vector<float> *mupNDOF_afterKinFit;
281 vector<float> *mumNDOF_afterKinFit;
282 vector<float> *mupD0_afterKinFit;
283 vector<float> *mumD0_afterKinFit;
284 vector<float> *mupNvalidHits_afterKinFit;
285 vector<float> *mumNvalidHits_afterKinFit;
286 vector<float> *mupIsGlobal_afterKinFit;
287 vector<float> *mumIsGlobal_afterKinFit;
288 vector<float> *mupIsTracker_afterKinFit;
289 vector<float> *mumIsTracker_afterKinFit;
290 vector<float> *ctauLambdaXY_afterKinFit;
291 vector<float> *LambdaBX_afterKinFit;
292 vector<float> *LambdaBY_afterKinFit;
293 vector<float> *LambdaBZ_afterKinFit;
294 vector<float> *LambdaBVtxX_afterKinFit;
295 vector<float> *LambdaBVtxY_afterKinFit;
296 vector<float> *LambdaBVtxZ_afterKinFit;
297 vector<float> *LambdaBNDOF_afterKinFit;
298 vector<float> *LambdaBchiSquared_afterKinFit;
299 vector<float> *LambdaBVtxProb_afterKinFit;
300 vector<float> *LambdaBPt_afterKinFit;
301 vector<float> *LambdaBEta_afterKinFit;
302 vector<float> *LambdaBPhi_afterKinFit;
303
304
305 int genLambdaJpsi;
306
307
308 //MC truth info
309 float truePriVtxX, truePriVtxY, truePriVtxZ;
310 float trueLambdaBPx, trueLambdaBPy,trueLambdaBPz;
311 float trueLambdaBDecayVtxX, trueLambdaBDecayVtxY, trueLambdaBDecayVtxZ;
312 float trueJpsiPx, trueJpsiPy, trueJpsiPz;
313 float trueJpsiDecayVtxX, trueJpsiDecayVtxY, trueJpsiDecayVtxZ;
314 float trueMumPx, trueMumPy, trueMumPz;
315 float trueMupPx, trueMupPy, trueMupPz;
316 float trueLambdaPx, trueLambdaPy, trueLambdaPz;
317 float trueLambdaDecayVtxX,trueLambdaDecayVtxY, trueLambdaDecayVtxZ;
318 float trueTrackPPx, trueTrackPPy, trueTrackPPz;
319 float trueTrackMPx, trueTrackMPy, trueTrackMPz;
320 float trueLxyPV;
321 float trueLxyBS;
322 float trueL3DPV;
323 float trueL3DBS;
324 float truectauLambdaBxyPV;
325 float truectauLambdaB3DPV;
326 float trueLambdaBMass;
327
328 vector<int> truthMatch, truthLambda, truthJpsi;
329
330
331 //parameters from config file
332 bool MCtruth_;
333 bool excludeKshortTracks_;
334 edm::InputTag triggerTag_;
335 edm::InputTag triggerEventTag_;
336
337
338 // ----------member data ---------------------------
339 };
340
341 //
342 // constants, enums and typedefs
343 //
344
345 //
346 // static data member definitions
347 //
348
349 //
350 // constructors and destructor
351 //
352 LambdaAnalyzer::LambdaAnalyzer(const edm::ParameterSet& iConfig)
353 :
354 tree_(0),
355 nLambdaB(0),
356 run(0),
357 lumiBlock(0),
358 event(0),
359 priVtxX(0),
360 priVtxY(0),
361 priVtxZ(0),
362 priVtxXE(0),
363 priVtxYE(0),
364 priVtxZE(0),
365 priVtxCL(0),
366 priVtxBSx(0),
367 priVtxBSy(0),
368 priVtxBSz(0),
369 priVtxBSxE(0),
370 priVtxBSyE(0),
371 priVtxBSzE(0),
372 priVtxBSCL(0),
373 beamSpotX(0),
374 beamSpotY(0),
375 beamSpotZ(0),
376 hlt_mu3(0),
377 hlt_mu5(0),
378 hlt_mu7(0),
379 hlt_mu9(0),
380 hlt_2mu0(0),
381 hlt_2mu3(0),
382 hlt_2mu0L2(0),
383 hlt_2mu3JPsi(0),
384 hlt_BJPsiMuMu(0),
385 hlt_mu0trk0(0),
386 hlt_mu3trk0(0),
387 hlt_mu0trkmu0(0),
388 hlt_mu3trkmu0(0),
389 hlt_L1muOpen(0),
390 hlt_L12muOpen(0),
391 hlt_L12muOpenTight(0),
392 hlt_mu0_L1muOpen(0),
393 hlt_mu3_L1muOpen(0),
394 hlt_mu5_L1muOpen(0),
395 hlt_mu5trk0(0),
396 hlt_mu0trkmu0_nocharge(0),
397 hlt_mu0trkmu0_OST(0),
398 hlt_2mu0_Qv1(0),
399 hlt_2mu0_QLSv1(0),
400 hlt_mu0trkmu0_OST_tightV1(0),
401 hlt_mu0trkmu0_OST_tightV2(0),
402 hlt_mu0trkmu0_OST_tightV3(0),
403 JpsiVtxX(0),
404 JpsiVtxY(0),
405 JpsiVtxZ(0),
406 LambdaVtxX(0),
407 LambdaVtxY(0),
408 LambdaVtxZ(0),
409 JpsiVtxProb(0),
410 JpsiVtxCL(0),
411 LambdaVtxProb(0),
412 LambdaVtxCL(0),
413 LambdaBMass(0),
414 LambdaBPx(0),
415 LambdaBPy(0),
416 LambdaBPz(0),
417 JpsiPx(0),
418 JpsiPy(0),
419 JpsiPz(0),
420 LambdaPx(0),
421 LambdaPy(0),
422 LambdaPz(0),
423 JpsiMass(0),
424 LambdaMass(0),
425 mupPx(0),
426 mupPy(0),
427 mupPz(0),
428 mumPx(0),
429 mumPy(0),
430 mumPz(0),
431 trackPPx(0),
432 trackPPy(0),
433 trackPPz(0),
434 trackMPx(0),
435 trackMPy(0),
436 trackMPz(0),
437 JpsiLambdaVtxDistXY(0),
438 JpsiLambdaVtxDist3D(0),
439 openingAngleLambdaB(0),
440 openingAngleLambda(0),
441 LxyPV(0),
442 LxyBS(0),
443 L3DPV(0),
444 L3DBS(0),
445 ctauLambdaBxyPV(0),
446 ctauLambdaBxyBS(0),
447 ctauLambdaB3DPV(0),
448 ctauLambdaB3DBS(0),
449 ctauLambdaXY(0),
450 dRmup(0),
451 dRmum(0),
452 dRtrackp(0),
453 dRtrackm(0),
454 dRLambdaVtx(0),
455 mupChi2(0),
456 mumChi2(0),
457 mupNDOF(0),
458 mumNDOF(0),
459 mupD0(0),
460 mumD0(0),
461 mupNvalidHits(0),
462 mumNvalidHits(0),
463 mupIsGlobal(0),
464 mumIsGlobal(0),
465 mupIsTracker(0),
466 mumIsTracker(0),
467 JpsiVtxX_afterKinFit(0),
468 JpsiVtxY_afterKinFit(0),
469 JpsiVtxZ_afterKinFit(0),
470 LambdaVtxX_afterKinFit(0),
471 LambdaVtxY_afterKinFit(0),
472 LambdaVtxZ_afterKinFit(0),
473 JpsiVtxProb_afterKinFit(0),
474 JpsiVtxCL_afterKinFit(0),
475 LambdaVtxProb_afterKinFit(0),
476 LambdaVtxCL_afterKinFit(0),
477 LambdaBMass_afterKinFit(0),
478 LambdaBPx_afterKinFit(0),
479 LambdaBPy_afterKinFit(0),
480 LambdaBPz_afterKinFit(0),
481 JpsiPx_afterKinFit(0),
482 JpsiPy_afterKinFit(0),
483 JpsiPz_afterKinFit(0),
484 LambdaPx_afterKinFit(0),
485 LambdaPy_afterKinFit(0),
486 LambdaPz_afterKinFit(0),
487 JpsiMass_afterKinFit(0),
488 LambdaMass_afterKinFit(0),
489 mupPx_afterKinFit(0),
490 mupPy_afterKinFit(0),
491 mupPz_afterKinFit(0),
492 mumPx_afterKinFit(0),
493 mumPy_afterKinFit(0),
494 mumPz_afterKinFit(0),
495 trackPPx_afterKinFit(0),
496 trackPPy_afterKinFit(0),
497 trackPPz_afterKinFit(0),
498 trackMPx_afterKinFit(0),
499 trackMPy_afterKinFit(0),
500 trackMPz_afterKinFit(0),
501 JpsiLambdaVtxDistXY_afterKinFit(0),
502 JpsiLambdaVtxDist3D_afterKinFit(0),
503 openingAngleLambdaB_afterKinFit(0),
504 openingAngleLambda_afterKinFit(0),
505 LxyPV_afterKinFit(0),
506 LxyBS_afterKinFit(0),
507 L3DPV_afterKinFit(0),
508 L3DBS_afterKinFit(0),
509 ctauLambdaBxyPV_afterKinFit(0),
510 ctauLambdaBxyBS_afterKinFit(0),
511 ctauLambdaB3DPV_afterKinFit(0),
512 ctauLambdaB3DBS_afterKinFit(0),
513 ctauLambdaXY_afterKinFit(0),
514 LambdaBX_afterKinFit(0),
515 LambdaBY_afterKinFit(0),
516 LambdaBZ_afterKinFit(0),
517 LambdaBVtxX_afterKinFit(0),
518 LambdaBVtxY_afterKinFit(0),
519 LambdaBVtxZ_afterKinFit(0),
520 LambdaBNDOF_afterKinFit(0),
521 LambdaBchiSquared_afterKinFit(0),
522 LambdaBVtxProb_afterKinFit(0),
523 LambdaBPt_afterKinFit(0),
524 LambdaBEta_afterKinFit(0),
525 LambdaBPhi_afterKinFit(0),
526 dRmup_afterKinFit(0),
527 dRmum_afterKinFit(0),
528 dRtrackp_afterKinFit(0),
529 dRtrackm_afterKinFit(0),
530 dRLambdaVtx_afterKinFit(0),
531 mupChi2_afterKinFit(0),
532 mumChi2_afterKinFit(0),
533 mupNDOF_afterKinFit(0),
534 mumNDOF_afterKinFit(0),
535 mupD0_afterKinFit(0),
536 mumD0_afterKinFit(0),
537 mupNvalidHits_afterKinFit(0),
538 mumNvalidHits_afterKinFit(0),
539 mupIsGlobal_afterKinFit(0),
540 mumIsGlobal_afterKinFit(0),
541 mupIsTracker_afterKinFit(0),
542 mumIsTracker_afterKinFit(0),
543 genLambdaJpsi(0),
544 truePriVtxX(0),
545 truePriVtxY(0),
546 truePriVtxZ(0),
547 trueLambdaBPx(0),
548 trueLambdaBPy(0),
549 trueLambdaBPz(0),
550 trueLambdaBDecayVtxX(0),
551 trueLambdaBDecayVtxY(0),
552 trueLambdaBDecayVtxZ(0),
553 trueJpsiPx(0),
554 trueJpsiPy(0),
555 trueJpsiPz(0),
556 trueJpsiDecayVtxX(0),
557 trueJpsiDecayVtxY(0),
558 trueJpsiDecayVtxZ(0),
559 trueMumPx(0),
560 trueMumPy(0),
561 trueMumPz(0),
562 trueMupPx(0),
563 trueMupPy(0),
564 trueMupPz(0),
565 trueLambdaPx(0),
566 trueLambdaPy(0),
567 trueLambdaPz(0),
568 trueLambdaDecayVtxX(0),
569 trueLambdaDecayVtxY(0),
570 trueLambdaDecayVtxZ(0),
571 trueTrackPPx(0),
572 trueTrackPPy(0),
573 trueTrackPPz(0),
574 trueTrackMPx(0),
575 trueTrackMPy(0),
576 trueTrackMPz(0),
577 trueLxyPV(0),
578 trueLxyBS(0),
579 trueL3DPV(0),
580 trueL3DBS(0),
581 truectauLambdaBxyPV(0),
582 truectauLambdaB3DPV(0),
583 trueLambdaBMass(0),
584 truthMatch(0), truthLambda(0), truthJpsi(0)
585
586 {
587 //now do what ever initialization is needed
588 edm::Service<TFileService> fs;
589 //then declare any histograms needed
590 MCtruth_ = iConfig.getParameter<bool>("MCtruth");
591 excludeKshortTracks_ = iConfig.getParameter<bool>("excludeKshortTracks");
592 triggerTag_ = iConfig.getParameter<edm::InputTag>("TriggerTag");
593 triggerEventTag_ = iConfig.getParameter<edm::InputTag>("TriggerEventTag");
594
595
596 }
597
598 LambdaAnalyzer::~LambdaAnalyzer()
599 {
600
601 // do anything here that needs to be done at desctruction time
602 // (e.g. close files, deallocate resources etc.)
603
604 }
605
606
607 //
608 // member functions
609 //
610
611 // ------------ method called to for each event ------------
612 void
613 LambdaAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
614 {
615
616 run = iEvent.id().run();
617 lumiBlock = iEvent.luminosityBlock();
618 event = iEvent.id().event();
619
620
621 float pi = 3.14159265;
622
623 using namespace edm;
624 using namespace std;
625 using namespace reco;
626
627 Handle<reco::VertexCompositeCandidateCollection> theKshorts;
628 iEvent.getByLabel( InputTag("generalV0Candidates:Kshort"), theKshorts );
629
630 Handle<reco::VertexCompositeCandidateCollection> theLambdas;
631 iEvent.getByLabel( InputTag("generalV0Candidates:Lambda"), theLambdas );
632
633 Handle<reco::MuonCollection> theMuons;
634 iEvent.getByLabel( InputTag("muons"), theMuons);
635
636 Handle<reco::BeamSpot> theBeamSpot;
637 iEvent.getByLabel(string("offlineBeamSpot"), theBeamSpot );
638
639 //get primary vertex
640 Handle<reco::VertexCollection> thePrimaryVertex;
641 iEvent.getByLabel("offlinePrimaryVertices", thePrimaryVertex);
642
643
644
645 reco::VertexCollection::const_iterator bestVtx = thePrimaryVertex->begin();
646
647
648 //get primary vertex with beam spot constraint
649 Handle<reco::VertexCollection>thePrimaryVertexBS;
650 iEvent.getByLabel("offlinePrimaryVerticesWithBS",thePrimaryVertexBS);
651
652
653 reco::VertexCollection::const_iterator bestVtxBS = thePrimaryVertexBS->begin();
654
655
656
657 //cout<<"event has "<<theLambdas->size()<<" Lambdas, "<<theKshorts->size()<<" Kshorts and "<<theMuons->size()<<" muons"<<endl;
658
659 //HLT results
660 edm::Handle<edm::TriggerResults> hltresults;
661 iEvent.getByLabel(triggerTag_, hltresults);
662 const edm::TriggerNames & triggerNames_ = iEvent.triggerNames(*hltresults);
663 int ntrigs=hltresults->size();
664 if (!MCtruth_) {
665 for (int itrig=0; itrig< ntrigs; itrig++) {
666 TString trigName = triggerNames_.triggerName(itrig);
667 int hltflag = (*hltresults)[itrig].accept();
668 if (trigName=="HLT_DoubleMu3_BJPsi") hlt_BJPsiMuMu = hltflag;
669 if (trigName=="HLT_Mu3") hlt_mu3 = hltflag;
670 if (trigName=="HLT_Mu5") hlt_mu5 = hltflag;
671 if (trigName=="HLT_Mu7") hlt_mu7 = hltflag;
672 if (trigName=="HLT_Mu9") hlt_mu9 = hltflag;
673 if (trigName=="HLT_DoubleMu0") hlt_2mu0 = hltflag;
674 if (trigName=="HLT_DoubleMu3") hlt_2mu3 = hltflag;
675 if (trigName=="HLT_DoubleMu3_Jpsi") hlt_2mu3JPsi = hltflag;
676 if (trigName=="HLT_Mu0_Track0_Jpsi") hlt_mu0trk0 = hltflag;
677 if (trigName=="HLT_Mu3_Track0_Jpsi") hlt_mu3trk0 = hltflag;
678 if (trigName=="HLT_Mu0_TkMu0_Jpsi") hlt_mu0trkmu0 = hltflag;
679 if (trigName=="HLT_Mu3_TkMu0_Jpsi") hlt_mu3trkmu0 = hltflag;
680 if (trigName=="HLT_L1MuOpen") hlt_L1muOpen = hltflag;
681 if (trigName=="HLT_L1DoubleMuOpen") hlt_L12muOpen = hltflag;
682 if (trigName=="HLT_L1DoubleMuOpen_Tight") hlt_L12muOpenTight = hltflag;
683 if (trigName=="HLT_L2DoubleMu0") hlt_2mu0L2 = hltflag;
684 if (trigName=="HLT_Mu0_L1MuOpen") hlt_mu0_L1muOpen = hltflag;
685 if (trigName=="HLT_Mu3_L1MuOpen") hlt_mu3_L1muOpen = hltflag;
686 if (trigName=="HLT_Mu5_L1MuOpen") hlt_mu5_L1muOpen = hltflag;
687 if (trigName=="HLT_Mu5_Track0_Jpsi") hlt_mu5trk0 = hltflag;
688 if (trigName=="HLT_Mu0_TkMu0_Jpsi_NoCharge") hlt_mu0trkmu0_nocharge = hltflag;
689
690 if (trigName=="HLT_Mu0_TkMu0_OST_Jpsi") hlt_mu0trkmu0_OST = hltflag;
691 if (trigName=="HLT_DoubleMu0_Quarkonium_v1") hlt_2mu0_Qv1 = hltflag;
692 if (trigName=="HLT_DoubleMu0_Quarkonium_LS_v1") hlt_2mu0_QLSv1 = hltflag;
693 if (trigName=="HLT_Mu0_TkMu0_OST_Jpsi_Tight_v1") hlt_mu0trkmu0_OST_tightV1 = hltflag;
694 if (trigName=="HLT_Mu0_TkMu0_OST_Jpsi_Tight_v2") hlt_mu0trkmu0_OST_tightV2 = hltflag;
695 if (trigName=="HLT_Mu0_TkMu0_OST_Jpsi_Tight_v3") hlt_mu0trkmu0_OST_tightV3 = hltflag;
696
697 }
698
699 } else {
700 for (int itrig=0; itrig< ntrigs; itrig++) {
701 TString trigName = triggerNames_.triggerName(itrig);
702 int hltflag = (*hltresults)[itrig].accept();
703 if (trigName=="HLT_L1DoubleMuOpen") hlt_L12muOpen = hltflag;
704 if (trigName=="HLT_Mu0_Track0_Jpsi") hlt_mu0trk0 = hltflag;
705 if (trigName=="HLT_DoubleMu0") hlt_2mu0 = hltflag;
706 }
707
708 edm::ESHandle<L1GtTriggerMenu> l1GtMenu;
709 iSetup.get<L1GtTriggerMenuRcd>().get(l1GtMenu);
710 const L1GtTriggerMenu* m_l1GtMenu = l1GtMenu.product();
711 (const_cast<L1GtTriggerMenu*> (m_l1GtMenu))->buildGtConditionMap(); //...ugly
712
713 const AlgorithmMap& algorithmAliasMap = l1GtMenu->gtAlgorithmAliasMap();
714
715 edm::InputTag m_l1GtReadoutRecordTag("gtDigis");
716 edm::Handle<L1GlobalTriggerReadoutRecord> gtReadoutRecord;
717 iEvent.getByLabel(m_l1GtReadoutRecordTag, gtReadoutRecord);
718 const L1GlobalTriggerReadoutRecord* gtReadoutRecordPtr = gtReadoutRecord.product();
719 const std::vector<bool>& gtDecisionWord = gtReadoutRecordPtr->decisionWord();
720
721 bool pass_L1SingleMuOpen = gtDecisionWord[(algorithmAliasMap.find("L1_SingleMuOpen")->second).algoBitNumber()];
722 bool pass_L1SingleMu0 = gtDecisionWord[(algorithmAliasMap.find("L1_SingleMu0")->second).algoBitNumber()];
723 bool pass_L1SingleMu3 = gtDecisionWord[(algorithmAliasMap.find("L1_SingleMu3")->second).algoBitNumber()];
724
725 hlt_mu0trkmu0 = (pass_L1SingleMuOpen||pass_L1SingleMu0||pass_L1SingleMu3) && pass_HLT_Mu0_TkMu0_Jpsi(iEvent);
726 hlt_mu0trkmu0_OST = hlt_mu0trkmu0;
727 hlt_mu0trkmu0_OST_tightV1 = (pass_L1SingleMuOpen||pass_L1SingleMu0||pass_L1SingleMu3) && pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v1(iEvent);
728 hlt_mu0trkmu0_OST_tightV2 = pass_L1SingleMu0 && hlt_mu0trkmu0_OST_tightV1;
729 hlt_mu0trkmu0_OST_tightV3 = pass_L1SingleMu0 && hlt_mu0trkmu0_OST_tightV1;
730 hlt_2mu0_Qv1 = pass_L1SingleMuOpen && pass_HLT_DoubleMu0_Quarkonium_v1(iEvent);
731 }
732
733
734 //loop on Lambdas and muons
735 //Lambdas
736 if(theLambdas->size()>0 && theMuons->size()>=2){
737 for(reco::VertexCompositeCandidateCollection::const_iterator iLambda=theLambdas->begin(); iLambda!=theLambdas->end();++iLambda){
738
739 const Candidate * dau0Cand = iLambda->daughter(0);
740 const Candidate * dau1Cand = iLambda->daughter(1);
741 //check for opposite charge of daughters
742 if(dau0Cand->charge()==dau1Cand->charge()){
743 //cout<<"Lambda0 daughters have same sign. Continue."<<endl;
744 continue;
745 }
746 CompositeCandidate LambdaCand;
747 LambdaCand.addDaughter(*dau0Cand);
748 LambdaCand.addDaughter(*dau1Cand);
749
750 AddFourMomenta add;
751 add.set(LambdaCand);
752
753 const reco::RecoChargedCandidate* daughter0 = static_cast<const reco::RecoChargedCandidate*>(iLambda->daughter(0));
754 const reco::RecoChargedCandidate* daughter1 = static_cast<const reco::RecoChargedCandidate*>(iLambda->daughter(1));
755 vector<RecoChargedCandidate> LambdaDaughters;
756 vector<TrackRef> LambdaDauTracks;
757
758 //check momentum and sort: pion first, proton second
759 if( daughter0->momentum().mag2() < daughter1->momentum().mag2() ){
760 LambdaDaughters.push_back(*(dynamic_cast<const reco::RecoChargedCandidate *> (daughter0)) );
761 LambdaDaughters.push_back(*(dynamic_cast<const reco::RecoChargedCandidate *> (daughter1)) );
762 }else{
763 LambdaDaughters.push_back(*(dynamic_cast<const reco::RecoChargedCandidate *> (daughter1)) );
764 LambdaDaughters.push_back(*(dynamic_cast<const reco::RecoChargedCandidate *> (daughter0)) );
765
766 }
767 for(unsigned int j=0; j<LambdaDaughters.size(); ++j){
768 LambdaDauTracks.push_back(LambdaDaughters[j].track() );
769 }
770
771 //loop on all Kshors in the event to check for common tracks with Lambda
772 if(excludeKshortTracks_){
773 for(reco::VertexCompositeCandidateCollection::const_iterator iKshort=theKshorts->begin(); iKshort!=theKshorts->end(); ++iKshort){
774 for(unsigned int idau=0; idau<LambdaDauTracks.size();++idau){
775 if(iKshort->daughter(0)->charge()==LambdaDauTracks[idau]->charge() &&
776 iKshort->daughter(0)->momentum()==LambdaDauTracks[idau]->momentum() ){
777 //cout<<"Lambda0 dau track matched with Ks dau track. Continue."<<endl;
778 continue;
779 }
780 if(iKshort->daughter(1)->charge()==LambdaDauTracks[idau]->charge() &&
781 iKshort->daughter(1)->momentum()==LambdaDauTracks[idau]->momentum() ){
782 //cout<<"Lambda0 dau track matched with Ks dau track. Continue."<<endl;
783 continue;
784 }
785 }
786 }//end of loop on Kshorts
787 }// end of if(excludeKshortTracks)
788
789
790 //now loop on muons
791 //first the positive muon
792 for(reco::MuonCollection::const_iterator iMuonP=theMuons->begin(); iMuonP!=theMuons->end(); ++iMuonP){
793 if(iMuonP->charge() == 1){
794 const Candidate &muPCand = *iMuonP;
795 TrackRef trackP = iMuonP->track();
796 if(trackP.isNull()){continue;}
797 bool match = false;
798 for(unsigned int j=0;j<LambdaDauTracks.size();++j){
799 if(trackP->charge()==LambdaDauTracks[j]->charge() &&
800 trackP->momentum() == LambdaDauTracks[j]->momentum() ){
801 //cout<<"Match found between muP and Lambda track with P = "<<trackP->momentum()<<" and "<<LambdaDauTracks[j]->momentum()<<endl;
802 match = true;
803 }
804 if(match) break;
805 }
806 if(match){// track is already used in making the Lambda
807 //cout<<"Match found between muP and Lambda track"<<endl;
808 continue;
809 }
810 //now the negative muon
811 for(reco::MuonCollection::const_iterator iMuonM=theMuons->begin(); iMuonM!=theMuons->end(); ++iMuonM){
812 if(iMuonM->charge() == -1){
813 const Candidate &muMCand = *iMuonM;
814 TrackRef trackM = iMuonM->track();
815 if(trackM.isNull()){continue;}
816 bool match = false;
817 for(unsigned int j=0;j<LambdaDauTracks.size();++j){
818 if(trackM->charge()==LambdaDauTracks[j]->charge() &&
819 trackM->momentum() == LambdaDauTracks[j]->momentum() ){
820 //cout<<"Match found between muM and Lambda track with P = "<<trackM->momentum()<<" and "<<LambdaDauTracks[j]->momentum()<<endl;
821 match = true;
822 }
823 if(match) break;
824 }
825 if(match){// track is already used in making the Lambda
826 //cout<<"Match found between muM and Lambda track"<<endl;
827 continue;
828 }
829 //four good tracks found
830
831 //KalmanVertexFitter to fit the four particles
832 edm::ESHandle<TransientTrackBuilder> theB;
833 iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB);
834
835 vector<TransientTrack> MuMuTT;
836 MuMuTT.push_back((*theB).build(&trackP));
837 MuMuTT.push_back((*theB).build(&trackM));
838
839 vector<TransientTrack> PiProtonTT;
840 PiProtonTT.push_back((*theB).build(&LambdaDauTracks[0]));
841 PiProtonTT.push_back((*theB).build(&LambdaDauTracks[1]));
842
843
844 KalmanVertexFitter kvf;
845 TransientVertex tvJpsi = kvf.vertex(MuMuTT);
846 TransientVertex tvLambda = kvf.vertex(PiProtonTT);
847 if(tvJpsi.isValid() && tvLambda.isValid()){
848 nLambdaB++;
849 mupChi2->push_back(trackP->chi2());
850 mumChi2->push_back(trackM->chi2());
851 mupNDOF->push_back(trackP->ndof());
852 mumNDOF->push_back(trackM->ndof());
853 mupD0->push_back(trackP->d0());
854 mumD0->push_back(trackM->d0());
855 mupNvalidHits->push_back(trackP->numberOfValidHits());
856 mumNvalidHits->push_back(trackM->numberOfValidHits());
857 mupIsGlobal->push_back(iMuonP->isGlobalMuon());
858 mumIsGlobal->push_back(iMuonM->isGlobalMuon());
859 mupIsTracker->push_back(iMuonP->isTrackerMuon());
860 mumIsTracker->push_back(iMuonM->isTrackerMuon());
861
862
863 CompositeCandidate JpsiCand;
864 JpsiCand.addDaughter(muPCand);
865 JpsiCand.addDaughter(muMCand);
866
867 AddFourMomenta add;
868 add.set(JpsiCand);
869 //create a composite candidate
870 CompositeCandidate LambdaB;
871 LambdaB.addDaughter(*dau0Cand);
872 LambdaB.addDaughter(*dau1Cand);
873 LambdaB.addDaughter(muPCand);
874 LambdaB.addDaughter(muMCand);
875
876 //AddFourMomenta add;
877 add.set(LambdaB);
878
879 //fill candidates information in vectors
880 LambdaBMass->push_back(LambdaB.mass());
881 LambdaBPx->push_back(LambdaB.px());
882 LambdaBPy->push_back(LambdaB.py());
883 LambdaBPz->push_back(LambdaB.pz());
884 JpsiPx->push_back(muPCand.px() + muMCand.px());
885 JpsiPy->push_back(muPCand.py() + muMCand.py());
886 JpsiPz->push_back(muPCand.pz() + muMCand.pz());
887 LambdaPx->push_back(dau0Cand->px() + dau1Cand->px());
888 LambdaPy->push_back(dau0Cand->py() + dau1Cand->py());
889 LambdaPz->push_back(dau0Cand->pz() + dau1Cand->pz());
890 JpsiMass->push_back(JpsiCand.mass());
891 LambdaMass->push_back(LambdaCand.mass());
892 mupPx->push_back(muPCand.px());
893 mupPy->push_back(muPCand.py());
894 mupPz->push_back(muPCand.pz());
895 mumPx->push_back(muMCand.px());
896 mumPy->push_back(muMCand.py());
897 mumPz->push_back(muMCand.pz());
898 for(uint k=0; k<iLambda->numberOfDaughters(); ++k){
899 if(iLambda->daughter(k)->charge()>0){
900
901 trackPPx->push_back(iLambda->daughter(k)->px());
902 trackPPy->push_back(iLambda->daughter(k)->py());
903 trackPPz->push_back(iLambda->daughter(k)->pz());
904 }
905 if(iLambda->daughter(k)->charge()<0){
906 trackMPx->push_back(iLambda->daughter(k)->px());
907 trackMPy->push_back(iLambda->daughter(k)->py());
908 trackMPz->push_back(iLambda->daughter(k)->pz());
909 }
910 }
911
912 ///////////////////////////////////////////////////////////
913 //primary vertex selection
914
915 //select best PV from the min opening angle for LambdaB
916 float minOpeningAngle = -999;
917 reco::VertexCollection::const_iterator bestVtx = thePrimaryVertex->begin();
918 for (reco::VertexCollection::const_iterator vtx = thePrimaryVertex->begin();
919 vtx != thePrimaryVertex->end(); ++vtx){
920 //opening for LambdaB, before kin fit
921 float openingAngleLambdaB_tmp = ( ( (tvJpsi.position().x()-vtx->x())*LambdaB.px()
922 + (tvJpsi.position().y()-vtx->y())*LambdaB.py()
923 + (tvJpsi.position().z()-vtx->z())*LambdaB.pz() ) /
924 sqrt( ( (tvJpsi.position().x()-vtx->x())*(tvJpsi.position().x()-vtx->x()) +
925 (tvJpsi.position().y()-vtx->y())*(tvJpsi.position().y()-vtx->y()) +
926 (tvJpsi.position().z()-vtx->z())*(tvJpsi.position().z()-vtx->z()) ) *
927 ( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() ) )
928 );
929 if(openingAngleLambdaB_tmp > minOpeningAngle){
930 minOpeningAngle = openingAngleLambdaB_tmp;
931 bestVtx = vtx;
932 }
933
934 }
935
936 //fill primary vertex info
937 //primary vertex
938 priVtxX->push_back(bestVtx->x());
939 priVtxY->push_back(bestVtx->y());
940 priVtxZ->push_back(bestVtx->z());
941 priVtxXE->push_back(bestVtx->xError());
942 priVtxYE->push_back(bestVtx->yError());
943 priVtxZE->push_back(bestVtx->zError());
944 priVtxCL->push_back(ChiSquaredProbability((double)(bestVtx->chi2()),(double)(bestVtx->ndof())) );
945
946 //select best primary vertex with BS constraint from min opening angle for LambdaB
947 minOpeningAngle = -999;
948 reco::VertexCollection::const_iterator bestVtxBS = thePrimaryVertexBS->begin();
949 for (reco::VertexCollection::const_iterator vtxBS = thePrimaryVertexBS->begin();
950 vtxBS != thePrimaryVertexBS->end(); ++vtxBS){
951 //opening for LambdaB, before kin fit
952 float openingAngleLambdaB_tmp = ( ( (tvJpsi.position().x()-vtxBS->x())*LambdaB.px()
953 + (tvJpsi.position().y()-vtxBS->y())*LambdaB.py()
954 + (tvJpsi.position().z()-vtxBS->z())*LambdaB.pz() ) /
955 sqrt( ( (tvJpsi.position().x()-vtxBS->x())*(tvJpsi.position().x()-vtxBS->x()) +
956 (tvJpsi.position().y()-vtxBS->y())*(tvJpsi.position().y()-vtxBS->y()) +
957 (tvJpsi.position().z()-vtxBS->z())*(tvJpsi.position().z()-vtxBS->z()) ) *
958 ( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() ) )
959 );
960 if(openingAngleLambdaB_tmp > minOpeningAngle){
961 minOpeningAngle = openingAngleLambdaB_tmp;
962 bestVtxBS = vtxBS;
963 }
964 }
965
966 //fill primary vertex with BS constraint info
967 priVtxBSx->push_back(bestVtxBS->x());
968 priVtxBSy->push_back(bestVtxBS->y());
969 priVtxBSz->push_back(bestVtxBS->z());
970 priVtxBSxE->push_back(bestVtxBS->xError());
971 priVtxBSyE->push_back(bestVtxBS->yError());
972 priVtxBSzE->push_back(bestVtxBS->zError());
973 priVtxBSCL->push_back(ChiSquaredProbability((double)(bestVtxBS->chi2()),(double)(bestVtxBS->ndof())) );
974
975
976
977 //end of primary vertex selection
978 //////////////////////////////////////////////
979
980
981 TVector3 LambdaVtx;
982 TVector3 JpsiVtx;
983 VertexDistanceXY vdistXY;
984 VertexDistance3D vdist3D;
985
986 LambdaVtx.SetXYZ(iLambda->vx(), iLambda->vy(), iLambda->vz());
987 JpsiVtx.SetXYZ(tvJpsi.position().x(), tvJpsi.position().y(), tvJpsi.position().z());
988
989 Measurement1D JpsiLambdaVtxDistXYMeasure = vdistXY.distance(tvJpsi, tvLambda);
990 Measurement1D JpsiLambdaVtxDist3DMeasure = vdist3D.distance(tvJpsi, tvLambda);
991
992 JpsiLambdaVtxDistXY->push_back(JpsiLambdaVtxDistXYMeasure.value());
993 JpsiLambdaVtxDist3D->push_back(JpsiLambdaVtxDist3DMeasure.value());
994
995 Vertex vertexJpsi = tvJpsi;
996 JpsiVtxProb->push_back(TMath::Prob(vertexJpsi.chi2(),(int)vertexJpsi.ndof()));
997 JpsiVtxCL->push_back(ChiSquaredProbability( (double)(vertexJpsi.chi2()),(double)(vertexJpsi.ndof()) ));
998
999
1000
1001 JpsiVtxX->push_back(tvJpsi.position().x());
1002 JpsiVtxY->push_back(tvJpsi.position().y());
1003 JpsiVtxZ->push_back(tvJpsi.position().z());
1004 LambdaVtxX->push_back(iLambda->vx());
1005 LambdaVtxY->push_back(iLambda->vy());
1006 LambdaVtxZ->push_back(iLambda->vz());
1007
1008 Vertex vertexLambda = tvLambda;
1009 LambdaVtxProb->push_back(TMath::Prob(vertexLambda.chi2(),(int)vertexLambda.ndof()));
1010 LambdaVtxCL->push_back(ChiSquaredProbability( (double)(vertexLambda.chi2()),(double)(vertexLambda.ndof()) ) );
1011
1012 //cos of angle between momentum and decay length
1013 //opening for LambdaB, before kin fit
1014 openingAngleLambdaB->push_back( ( (tvJpsi.position().x()-bestVtx->x())*LambdaB.px()
1015 + (tvJpsi.position().y()-bestVtx->y())*LambdaB.py()
1016 + (tvJpsi.position().z()-bestVtx->z())*LambdaB.pz() ) /
1017 sqrt( ( (tvJpsi.position().x()-bestVtx->x())*(tvJpsi.position().x()-bestVtx->x()) +
1018 (tvJpsi.position().y()-bestVtx->y())*(tvJpsi.position().y()-bestVtx->y()) +
1019 (tvJpsi.position().z()-bestVtx->z())*(tvJpsi.position().z()-bestVtx->z()) ) *
1020 ( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() ) )
1021 );
1022
1023 openingAngleLambda->push_back( ( (iLambda->vx()-tvJpsi.position().x())*(dau0Cand->px() + dau1Cand->px()) +
1024 (iLambda->vy()-tvJpsi.position().y())*(dau0Cand->py() + dau1Cand->py()) +
1025 (iLambda->vz()-tvJpsi.position().z())*(dau0Cand->pz() + dau1Cand->pz()) ) /
1026 sqrt( (
1027 (iLambda->vx()-tvJpsi.position().x())*(iLambda->vx()-tvJpsi.position().x()) +
1028 (iLambda->vy()-tvJpsi.position().y())*(iLambda->vy()-tvJpsi.position().y()) +
1029 (iLambda->vz()-tvJpsi.position().z())*(iLambda->vz()-tvJpsi.position().z()) ) *
1030 (
1031 (dau0Cand->px() + dau1Cand->px())*(dau0Cand->px() + dau1Cand->px()) +
1032 (dau0Cand->py() + dau1Cand->py())*(dau0Cand->py() + dau1Cand->py()) +
1033 (dau0Cand->pz() + dau1Cand->pz())*(dau0Cand->pz() + dau1Cand->pz())
1034 )
1035 )
1036 );
1037
1038 //decay lendth L and proper decay length ctau=L*m/p with beam spot and PV
1039 LxyPV->push_back( ( (tvJpsi.position().x()-bestVtx->x())*LambdaB.px() + (tvJpsi.position().y()-bestVtx->y())*LambdaB.py() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py()) );
1040 LxyBS->push_back( ( (tvJpsi.position().x()-theBeamSpot->position().x())*LambdaB.px() + (tvJpsi.position().y()-theBeamSpot->position().y())*LambdaB.py() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() ) );
1041 L3DPV->push_back( ( (tvJpsi.position().x()-bestVtx->x())*LambdaB.px() + (tvJpsi.position().y()-bestVtx->y())*LambdaB.py() + (tvJpsi.position().z()-bestVtx->z())*LambdaB.pz()) /sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() ) );
1042 L3DBS->push_back( ( (tvJpsi.position().x()-theBeamSpot->position().x())*LambdaB.px() + (tvJpsi.position().y()-theBeamSpot->position().y())*LambdaB.py() +(tvJpsi.position().z()-theBeamSpot->position().z())*LambdaB.pz() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() ) );
1043
1044 float LxyPVtmp = ( (tvJpsi.position().x()-bestVtx->x())*LambdaB.px() + (tvJpsi.position().y()-bestVtx->y())*LambdaB.py() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py());
1045 float LxyBStmp = ( (tvJpsi.position().x()-theBeamSpot->position().x())*LambdaB.px() + (tvJpsi.position().y()-theBeamSpot->position().y())*LambdaB.py() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() );
1046 float L3DPVtmp = ( (tvJpsi.position().x()-bestVtx->x())*LambdaB.px() + (tvJpsi.position().y()-bestVtx->y())*LambdaB.py() + (tvJpsi.position().z()-bestVtx->z())*LambdaB.pz()) /sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() );
1047 float L3DBStmp = ( (tvJpsi.position().x()-theBeamSpot->position().x())*LambdaB.px() + (tvJpsi.position().y()-theBeamSpot->position().y())*LambdaB.py() +(tvJpsi.position().z()-theBeamSpot->position().z())*LambdaB.pz() )/sqrt( LambdaB.px()*LambdaB.px() + LambdaB.py()*LambdaB.py() + LambdaB.pz()*LambdaB.pz() );
1048
1049 ctauLambdaBxyPV->push_back(LxyPVtmp*LambdaB.mass()/LambdaB.pt());
1050 ctauLambdaBxyBS->push_back(LxyBStmp*LambdaB.mass()/LambdaB.pt());
1051 ctauLambdaB3DPV->push_back(L3DPVtmp*LambdaB.mass()/LambdaB.p());
1052 ctauLambdaB3DBS->push_back(L3DBStmp*LambdaB.mass()/LambdaB.p());
1053
1054 float LambdaLxytmp = ( (iLambda->vx()-tvJpsi.position().x())*(dau0Cand->px()+dau1Cand->px()) +
1055 (iLambda->vy()-tvJpsi.position().y())*(dau0Cand->py()+dau1Cand->py()) ) /
1056 sqrt( (dau0Cand->px()+dau1Cand->px())*(dau0Cand->px()+dau1Cand->px()) + (dau0Cand->py()+dau1Cand->py())*(dau0Cand->py()+dau1Cand->py()) );
1057 ctauLambdaXY->push_back(LambdaLxytmp*LambdaCand.mass()/LambdaCand.pt() );
1058
1059
1060 ///////////////////////////////////////////////
1061 // Kinematic fit
1062
1063 KinematicParticleFactoryFromTransientTrack pFactory;
1064 KinematicParticleVertexFitter Fitter;
1065 KinematicParticleFitter constFitter;
1066 KinematicConstrainedVertexFitter constVertexFitter;
1067
1068 const ParticleMass muonMass(0.1056583);
1069 float muonSigma = 1E-10;
1070
1071 const ParticleMass protonMass(0.938272013);
1072 float protonSigma = 1E-10;
1073
1074 const ParticleMass pionMass(0.139570);
1075 float pionSigma = 0.000016;
1076
1077 vector<RefCountedKinematicParticle> allParticlesMu;
1078 allParticlesMu.push_back(pFactory.particle (MuMuTT[0], muonMass, float(0), float(0), muonSigma));
1079 allParticlesMu.push_back(pFactory.particle (MuMuTT[1], muonMass, float(0), float(0), muonSigma));
1080 RefCountedKinematicTree JpsiTree = Fitter.fit(allParticlesMu);
1081 if(JpsiTree->isEmpty()) {push_anything_back(); continue;}
1082
1083 double nominalJpsiMass = 3.096916;
1084 float jpsiMsigma = 0.00004;
1085 KinematicConstraint *jpsi_const = new MassKinematicConstraint( nominalJpsiMass, jpsiMsigma);
1086
1087 JpsiTree = constFitter.fit(jpsi_const,JpsiTree);
1088
1089 if(JpsiTree->isEmpty()) {push_anything_back();continue;}
1090
1091 JpsiTree->movePointerToTheTop();
1092 RefCountedKinematicParticle Jpsi_branch = JpsiTree->currentParticle();
1093
1094 if (!Jpsi_branch->currentState().isValid()) {push_anything_back();continue;}
1095
1096 vector<RefCountedKinematicParticle> allLambdaDaughters;
1097 allLambdaDaughters.push_back(pFactory.particle (PiProtonTT[0], pionMass, float(0), float(0), pionSigma));
1098 allLambdaDaughters.push_back(pFactory.particle (PiProtonTT[1], protonMass, float(0), float(0), protonSigma));
1099 RefCountedKinematicTree LambdaTree = Fitter.fit(allLambdaDaughters);
1100 if(LambdaTree->isEmpty()) {push_anything_back();continue;}
1101
1102 double nominalLambdaMass = 1.115683;
1103 float lambdaMsigma = 0.000006;
1104 KinematicConstraint *lambda_const = new MassKinematicConstraint( nominalLambdaMass, lambdaMsigma);
1105
1106 LambdaTree = constFitter.fit(lambda_const,LambdaTree);
1107
1108 if(LambdaTree->isEmpty()) {push_anything_back();continue;}
1109
1110 LambdaTree->movePointerToTheTop();
1111 RefCountedKinematicParticle Lambda_branch = LambdaTree->currentParticle();
1112 if (!Lambda_branch->currentState().isValid()) {push_anything_back();continue;}
1113
1114 LambdaTree->movePointerToTheFirstChild();
1115 RefCountedKinematicParticle Pion_branch = LambdaTree->currentParticle();
1116 if (!Pion_branch->currentState().isValid()) {push_anything_back();continue;}
1117
1118 LambdaTree->movePointerToTheNextChild();
1119 RefCountedKinematicParticle Proton_branch = LambdaTree->currentParticle();
1120 if (!Proton_branch->currentState().isValid()) {push_anything_back();continue;}
1121
1122
1123 vector<RefCountedKinematicParticle> allLambdaBDaughters;
1124 allLambdaBDaughters.push_back(pFactory.particle (MuMuTT[0], muonMass, float(0), float(0), muonSigma));
1125 allLambdaBDaughters.push_back(pFactory.particle (MuMuTT[1], muonMass, float(0), float(0), muonSigma));
1126 allLambdaBDaughters.push_back(Lambda_branch);
1127
1128 MultiTrackKinematicConstraint *jpsi_mtc = new TwoTrackMassKinematicConstraint(nominalJpsiMass);
1129 RefCountedKinematicTree LambdaBTree = constVertexFitter.fit(allLambdaBDaughters,jpsi_mtc);
1130 if(LambdaBTree->isEmpty()) {push_anything_back();continue;}
1131
1132
1133 LambdaBTree->movePointerToTheTop();
1134 RefCountedKinematicParticle LambdaB_branch = LambdaBTree->currentParticle();
1135 if (!LambdaB_branch->currentState().isValid()) {push_anything_back();continue;}
1136 RefCountedKinematicVertex LambdaB_DecayVtx = LambdaBTree->currentDecayVertex();
1137
1138 AlgebraicVector7 LambdaB_par = LambdaB_branch->currentState().kinematicParameters().vector();
1139 AlgebraicMatrix77 LambdaB_err = LambdaB_branch->currentState().kinematicParametersError().matrix();
1140 GlobalVector LambdaBvector(LambdaB_par[3], LambdaB_par[4], LambdaB_par[5]); // the fitted momentum vector of the LambdaB
1141
1142
1143 LambdaBTree->movePointerToTheFirstChild();
1144 RefCountedKinematicParticle MuP_branch = LambdaBTree->currentParticle();
1145 if (!MuP_branch->currentState().isValid()) {push_anything_back();continue;}
1146
1147 LambdaBTree->movePointerToTheNextChild();
1148 RefCountedKinematicParticle MuM_branch = LambdaBTree->currentParticle();
1149 if (!MuM_branch->currentState().isValid()) {push_anything_back();continue;}
1150
1151 LambdaBTree->movePointerToTheNextChild();
1152 Lambda_branch = LambdaBTree->currentParticle();
1153 if (!Lambda_branch->currentState().isValid()) {push_anything_back();continue;}
1154
1155 //fill candidates information in vectors
1156 LambdaBMass_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().mass());
1157 LambdaBPx_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().momentum().x());
1158 LambdaBPy_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().momentum().y());
1159 LambdaBPz_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().momentum().z());
1160 LambdaBX_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().position().x());
1161 LambdaBY_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().position().y());
1162 LambdaBZ_afterKinFit->push_back(LambdaB_branch->currentState().kinematicParameters().position().z());
1163 LambdaBEta_afterKinFit->push_back(LambdaBvector.eta());
1164 LambdaBPhi_afterKinFit->push_back(LambdaBvector.phi());
1165 LambdaBPt_afterKinFit->push_back(LambdaBvector.perp());
1166 LambdaBNDOF_afterKinFit->push_back((int)LambdaB_branch->degreesOfFreedom());
1167 LambdaBchiSquared_afterKinFit->push_back(LambdaB_branch->chiSquared());
1168 LambdaBVtxProb_afterKinFit->push_back(TMath::Prob(LambdaB_branch->chiSquared(), (int)LambdaB_branch->degreesOfFreedom()) );
1169 LambdaBVtxX_afterKinFit->push_back(LambdaB_DecayVtx->position().x());
1170 LambdaBVtxY_afterKinFit->push_back(LambdaB_DecayVtx->position().y());
1171 LambdaBVtxZ_afterKinFit->push_back(LambdaB_DecayVtx->position().z());
1172
1173 mumPx_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[3]);
1174 mumPy_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[4]);
1175 mumPz_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[5]);
1176
1177 mupPx_afterKinFit->push_back((MuP_branch->currentState().kinematicParameters().vector())[3]);
1178 mupPy_afterKinFit->push_back((MuP_branch->currentState().kinematicParameters().vector())[4]);
1179 mupPz_afterKinFit->push_back((MuP_branch->currentState().kinematicParameters().vector())[5]);
1180
1181 JpsiPx_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[3]+(MuP_branch->currentState().kinematicParameters().vector())[3]);
1182 JpsiPy_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[4]+(MuP_branch->currentState().kinematicParameters().vector())[4]);
1183 JpsiPz_afterKinFit->push_back((MuM_branch->currentState().kinematicParameters().vector())[5]+(MuP_branch->currentState().kinematicParameters().vector())[5]);
1184
1185 LambdaPx_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[3]);
1186 LambdaPy_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[4]);
1187 LambdaPz_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[5]);
1188
1189 JpsiVtxX_afterKinFit->push_back((LambdaB_branch->currentState().kinematicParameters().vector())[0]);
1190 JpsiVtxY_afterKinFit->push_back((LambdaB_branch->currentState().kinematicParameters().vector())[1]);
1191 JpsiVtxZ_afterKinFit->push_back((LambdaB_branch->currentState().kinematicParameters().vector())[2]);
1192
1193 LambdaVtxX_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[0]);
1194 LambdaVtxY_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[1]);
1195 LambdaVtxZ_afterKinFit->push_back((Lambda_branch->currentState().kinematicParameters().vector())[2]);
1196
1197 LambdaVtxProb_afterKinFit->push_back(TMath::Prob(Lambda_branch->chiSquared(),(int)Lambda_branch->degreesOfFreedom()));
1198 JpsiVtxProb_afterKinFit->push_back(TMath::Prob(Jpsi_branch->chiSquared(),(int)Jpsi_branch->degreesOfFreedom()));
1199
1200 float LambdaBVtxX_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[0];
1201 float LambdaBVtxY_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[1];
1202 float LambdaBVtxZ_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[2];
1203
1204 float LambdaBPx_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[3];
1205 float LambdaBPy_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[4];
1206 float LambdaBPz_fit_tmp = (LambdaB_branch->currentState().kinematicParameters().vector())[5];
1207
1208 float PVx_tmp = bestVtx->x();
1209 float PVy_tmp = bestVtx->y();
1210 float PVz_tmp = bestVtx->z();
1211
1212 float LambdaBMass_fit_tmp = LambdaB_branch->currentState().kinematicParameters().mass();
1213 float LambdaMass_fit_tmp = Lambda_branch->currentState().kinematicParameters().mass();
1214
1215
1216 float oAngleLambdaB_fit_tmp =
1217 ( (LambdaBVtxX_fit_tmp - PVx_tmp)*LambdaBPx_fit_tmp
1218 + (LambdaBVtxY_fit_tmp - PVy_tmp)*LambdaBPy_fit_tmp
1219 + (LambdaBVtxZ_fit_tmp - PVz_tmp)*LambdaBPz_fit_tmp ) /
1220 sqrt( ((LambdaBVtxX_fit_tmp - PVx_tmp)*(LambdaBVtxX_fit_tmp - PVx_tmp) +
1221 (LambdaBVtxY_fit_tmp - PVy_tmp)*(LambdaBVtxY_fit_tmp - PVy_tmp) +
1222 (LambdaBVtxZ_fit_tmp - PVz_tmp)*(LambdaBVtxZ_fit_tmp - PVz_tmp)) * (
1223 (LambdaBPx_fit_tmp * LambdaBPx_fit_tmp) +
1224 (LambdaBPy_fit_tmp * LambdaBPy_fit_tmp) +
1225 (LambdaBPz_fit_tmp * LambdaBPz_fit_tmp) ) );
1226
1227
1228 float LambdaVtxX_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[0];
1229 float LambdaVtxY_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[1];
1230 float LambdaVtxZ_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[2];
1231
1232 float LambdaPx_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[3];
1233 float LambdaPy_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[4];
1234 float LambdaPz_fit_tmp = (Lambda_branch->currentState().kinematicParameters().vector())[5];
1235
1236
1237
1238 // float oAngleLambda_fit_tmp =
1239 // ( (LambdaVtxX_fit_tmp - LambdaBVtxX_fit_tmp)*LambdaPx_fit_tmp
1240 // + (LambdaVtxY_fit_tmp - LambdaBVtxY_fit_tmp)*LambdaPy_fit_tmp
1241 // + (LambdaVtxZ_fit_tmp - LambdaBVtxZ_fit_tmp)*LambdaPz_fit_tmp )/
1242 // sqrt( ((LambdaVtxX_fit_tmp - LambdaBVtxX_fit_tmp)*(LambdaVtxX_fit_tmp - LambdaBVtxX_fit_tmp) +
1243 // (LambdaVtxY_fit_tmp - LambdaBVtxY_fit_tmp)*(LambdaVtxY_fit_tmp - LambdaBVtxY_fit_tmp) +
1244 // (LambdaVtxZ_fit_tmp - LambdaBVtxZ_fit_tmp)*(LambdaVtxZ_fit_tmp - LambdaBVtxZ_fit_tmp)) * (
1245 // (LambdaPx_fit_tmp*LambdaPx_fit_tmp) +
1246 // (LambdaPy_fit_tmp*LambdaPy_fit_tmp) +
1247 // (LambdaPz_fit_tmp*LambdaPz_fit_tmp) ) );
1248
1249
1250
1251
1252 openingAngleLambdaB_afterKinFit->push_back(oAngleLambdaB_fit_tmp);
1253 //openingAngleLambda_afterKinFit->push_back(oAngleLambda_fit_tmp);
1254
1255
1256
1257 //cout<<"two muons have valid vtx"<<endl;
1258 tmpNvtx++;
1259
1260
1261 float L3DPV_fit_tmp = ( (LambdaBVtxX_fit_tmp-PVx_tmp)*LambdaBPx_fit_tmp + (LambdaBVtxY_fit_tmp-PVy_tmp)*LambdaBPy_fit_tmp + (LambdaBVtxZ_fit_tmp-PVz_tmp)*LambdaBPz_fit_tmp )/
1262 sqrt(LambdaBPx_fit_tmp*LambdaBPx_fit_tmp + LambdaBPy_fit_tmp*LambdaBPy_fit_tmp + LambdaBPz_fit_tmp*LambdaBPz_fit_tmp);
1263 ctauLambdaB3DPV_afterKinFit->push_back(L3DPV_fit_tmp*LambdaBMass_fit_tmp/sqrt(LambdaBPx_fit_tmp*LambdaBPx_fit_tmp+LambdaBPy_fit_tmp*LambdaBPy_fit_tmp+LambdaBPz_fit_tmp*LambdaBPz_fit_tmp) );
1264
1265 float LxyPV_fit_tmp = ( (LambdaBVtxX_fit_tmp-PVx_tmp)*LambdaBPx_fit_tmp + (LambdaBVtxY_fit_tmp-PVy_tmp)*LambdaBPy_fit_tmp )/
1266 sqrt(LambdaBPx_fit_tmp*LambdaBPx_fit_tmp + LambdaBPy_fit_tmp*LambdaBPy_fit_tmp);
1267 ctauLambdaBxyPV_afterKinFit->push_back(LxyPV_fit_tmp*LambdaBMass_fit_tmp/sqrt(LambdaBPx_fit_tmp*LambdaBPx_fit_tmp+LambdaBPy_fit_tmp*LambdaBPy_fit_tmp) );
1268
1269 /*
1270 float LambdaLxy_fit_tmp = ( (LambdaVtxX_fit_tmp-LambdaBVtxX_fit_tmp)*LambdaPx_fit_tmp + (LambdaVtxY_fit_tmp-LambdaBVtxY_fit_tmp)*LambdaPy_fit_tmp )/
1271 sqrt(LambdaPx_fit_tmp*LambdaPx_fit_tmp + LambdaPy_fit_tmp*LambdaPy_fit_tmp);
1272
1273 ctauLambdaXY_afterKinFit->push_back(LambdaLxy_fit_tmp*LambdaMass_fit_tmp/sqrt(LambdaPx_fit_tmp*LambdaPx_fit_tmp+LambdaPy_fit_tmp*LambdaPy_fit_tmp) );
1274
1275 */
1276
1277
1278 }//end of if MuMu vertex is valid
1279 }//end of if(iMuonP->charge() == -1)
1280 }//end of loop of iMuonP
1281 }//end of if(iMuonP->charge() == 1)
1282 }//end of loop of iMuonP
1283 }//end of loop on Lambdas
1284 }//end of loop on Lambdas and muons combined (if ...)
1285
1286 // Beam Spot
1287 beamSpotX = theBeamSpot->position().x();
1288 beamSpotY = theBeamSpot->position().y();
1289 beamSpotZ = theBeamSpot->position().z();
1290
1291
1292
1293 ///////////////////////////////
1294 /// MC part
1295 if(nLambdaB > 0 && MCtruth_){
1296
1297 // Beam Spot
1298 beamSpotX = theBeamSpot->position().x();
1299 beamSpotY = theBeamSpot->position().y();
1300 beamSpotZ = theBeamSpot->position().z();
1301
1302 Handle<GenParticleCollection>theGenParticles;
1303 iEvent.getByLabel("genParticlePlusGEANT", theGenParticles);
1304
1305 genLambdaJpsi = -1;
1306 for(size_t iGenParticle=0; iGenParticle<theGenParticles->size();++iGenParticle){
1307 const Candidate & genLambdab = (*theGenParticles)[iGenParticle];
1308 if(abs( genLambdab.pdgId() )==5122){//if Lambda^0_b is generated
1309 //cout<<"found gen Lambda_b";
1310 int iJpsi(-1), iLambda(-1);
1311 bool wrong = false;
1312 for(uint i=0; i<genLambdab.numberOfDaughters();++i){
1313 //check Lambda_b for Jpsi and Lambda daughters
1314 const Candidate * genDau = genLambdab.daughter(i);
1315 //cout<<" ="<<genDau->pdgId();
1316 int imu1(-1), imu2(-1), ipi(-1), ip(-1);
1317 for(uint j=0; j<genDau->numberOfDaughters();++j){
1318 //cout<<" =="<<genDau->daughter(j)->pdgId();
1319 if(genDau->daughter(j)->pdgId()==13 && genDau->pdgId()==443)imu1=j;
1320 if(genDau->daughter(j)->pdgId()==-13 && genDau->pdgId()==443)imu2=j;
1321 if(abs(genDau->daughter(j)->pdgId())==211 && abs(genDau->pdgId())==3122 )ipi=j;
1322 if(abs(genDau->daughter(j)->pdgId())==2212 && abs(genDau->pdgId())==3122)ip=j;
1323
1324 if(genDau->pdgId()==443 && abs(genDau->daughter(j)->pdgId())!=13 && genDau->daughter(j)->pdgId()!=22)wrong=true;
1325 if(abs(genDau->pdgId())==3122 && abs(genDau->daughter(j)->pdgId())!=211 && abs(genDau->daughter(j)->pdgId())!=2212 && genDau->daughter(j)->pdgId()!=22)wrong = true;
1326 }//end of loop on daughters of daughters
1327 if(genDau->pdgId()!=443 && abs(genDau->pdgId())!=3122 && genDau->pdgId()!=22)wrong = true;
1328 if(imu1!=-1 && imu2!=-1 && !wrong)iJpsi = i;
1329 if(ipi!=-1 && ip!=-1 && !wrong)iLambda = i;
1330 }//end of loop on Lambda^0_b MC daughters
1331 //cout<<endl;
1332 if(iJpsi!=-1 && iLambda!=-1 && !wrong){
1333 //cout<<"setting genLambdaJpsi = 1"<<endl;
1334 genLambdaJpsi = 1;
1335
1336 //fill MC truth info
1337 trueLambdaBMass = genLambdab.p4().M();
1338
1339 //write out info from daughters
1340 const Candidate * genJpsi = genLambdab.daughter(iJpsi);
1341 const Candidate * genLambda = genLambdab.daughter(iLambda);
1342
1343 truePriVtxX = genLambdab.vx();
1344 truePriVtxY = genLambdab.vy();
1345 truePriVtxZ = genLambdab.vz();
1346
1347 trueLambdaBPx = genLambdab.px();
1348 trueLambdaBPy = genLambdab.py();
1349 trueLambdaBPz = genLambdab.pz();
1350
1351 trueLambdaBDecayVtxX = genJpsi->vx();
1352 trueLambdaBDecayVtxY = genJpsi->vy();
1353 trueLambdaBDecayVtxZ = genJpsi->vz();
1354
1355 trueJpsiPx = genJpsi->px();
1356 trueJpsiPy = genJpsi->py();
1357 trueJpsiPz = genJpsi->pz();
1358
1359 for(uint j=0; j<genJpsi->numberOfDaughters(); ++j){
1360 if(genJpsi->daughter(j)->pdgId()==13){//13 is a mu-
1361 trueMumPx = genJpsi->daughter(j)->px();
1362 trueMumPy = genJpsi->daughter(j)->py();
1363 trueMumPz = genJpsi->daughter(j)->pz();
1364
1365 trueJpsiDecayVtxX = genJpsi->daughter(j)->vx();
1366 trueJpsiDecayVtxY = genJpsi->daughter(j)->vy();
1367 trueJpsiDecayVtxZ = genJpsi->daughter(j)->vz();
1368 }
1369 if(genJpsi->daughter(j)->pdgId()==-13){//-13 is a mu+
1370 trueMupPx = genJpsi->daughter(j)->px();
1371 trueMupPy = genJpsi->daughter(j)->py();
1372 trueMupPz = genJpsi->daughter(j)->pz();
1373 }
1374 }
1375 trueLambdaPx = genLambda->px();
1376 trueLambdaPy = genLambda->py();
1377 trueLambdaPz = genLambda->pz();
1378
1379 //cout<<"starting daughters loop with genLambda id = "<<genLambda->pdgId()<<" and with "<<genLambda->numberOfDaughters()<<" daughters"<<endl;
1380 for(uint j=0; j<genLambda->numberOfDaughters(); ++j){
1381 //cout<<"genLambda daughter has id = "<<genLambda->daughter(j)->pdgId()<<" with ndau "<<genLambda->daughter(j)->numberOfDaughters()<<endl;
1382 //cout<<"no granddaughters, filling genLambda info"<<endl;
1383 if(genLambda->daughter(j)->charge()>0){
1384 trueTrackPPx = genLambda->daughter(j)->px();
1385 trueTrackPPy = genLambda->daughter(j)->py();
1386 trueTrackPPz = genLambda->daughter(j)->pz();
1387 trueLambdaDecayVtxX = genLambda->daughter(j)->vx();
1388 trueLambdaDecayVtxY = genLambda->daughter(j)->vy();
1389 trueLambdaDecayVtxZ = genLambda->daughter(j)->vz();
1390 }
1391 if(genLambda->daughter(j)->charge()<0){
1392 trueTrackMPx = genLambda->daughter(j)->px();
1393 trueTrackMPy = genLambda->daughter(j)->py();
1394 trueTrackMPz = genLambda->daughter(j)->pz();
1395 }
1396 }
1397
1398 //true decay lengths
1399 trueLxyPV = ( (trueLambdaBDecayVtxX-truePriVtxX)*trueLambdaBPx +
1400 (trueLambdaBDecayVtxY-truePriVtxY)*trueLambdaBPy ) /
1401 sqrt(trueLambdaBPx*trueLambdaBPx + trueLambdaBPy*trueLambdaBPy);
1402
1403 trueL3DPV = ( (trueLambdaBDecayVtxX-truePriVtxX)*trueLambdaBPx +
1404 (trueLambdaBDecayVtxY-truePriVtxY)*trueLambdaBPy +
1405 (trueLambdaBDecayVtxZ-truePriVtxZ)*trueLambdaBPz) /
1406 sqrt(trueLambdaBPx*trueLambdaBPx + trueLambdaBPy*trueLambdaBPy + trueLambdaBPz*trueLambdaBPz);
1407
1408 truectauLambdaBxyPV = trueLxyPV*trueLambdaBMass/sqrt(trueLambdaBPx*trueLambdaBPx + trueLambdaBPy*trueLambdaBPy );
1409 truectauLambdaB3DPV = trueL3DPV*trueLambdaBMass/sqrt(trueLambdaBPx*trueLambdaBPx + trueLambdaBPy*trueLambdaBPy + trueLambdaBPz*trueLambdaBPz);
1410
1411
1412 //determine MC truth
1413
1414 //calculate true eta and phi for all tracks
1415 float trueMupPhi = atan(trueMupPy/trueMupPx);
1416 if(trueMupPx<0 && trueMupPy<0) trueMupPhi -=pi;
1417 if(trueMupPx<0 && trueMupPy>0) trueMupPhi +=pi;
1418 float trueMupP = sqrt(trueMupPx*trueMupPx + trueMupPy*trueMupPy + trueMupPz*trueMupPz);
1419 float trueMupEta = 0.5*(log( (trueMupP + trueMupPz)/(trueMupP - trueMupPz)) );
1420
1421 float trueMumPhi = atan(trueMumPy/trueMumPx);
1422 if(trueMumPx<0 && trueMumPy<0) trueMumPhi -=pi;
1423 if(trueMumPx<0 && trueMumPy>0) trueMumPhi +=pi;
1424 float trueMumP = sqrt( trueMumPx*trueMumPx + trueMumPy*trueMumPy + trueMumPz*trueMumPz);
1425 float trueMumEta = 0.5*(log( (trueMumP + trueMumPz)/(trueMumP - trueMumPz)) );
1426
1427 float trueTrackPPhi = atan(trueTrackPPy/trueTrackPPx);
1428 if ( trueTrackPPx < 0 && trueTrackPPy < 0 ) trueTrackPPhi -= pi;
1429 if ( trueTrackPPx < 0 && trueTrackPPy > 0 ) trueTrackPPhi += pi;
1430 float trueTrackPP = sqrt( trueTrackPPx*trueTrackPPx + trueTrackPPy*trueTrackPPy + trueTrackPPz*trueTrackPPz );
1431 float trueTrackPEta = 0.5*log( (trueTrackPP + trueTrackPPz)/(trueTrackPP - trueTrackPPz) );
1432
1433 float trueTrackMPhi = atan(trueTrackMPy/trueTrackMPx);
1434 if ( trueTrackMPx < 0 && trueTrackMPy < 0 ) trueTrackMPhi -= pi;
1435 if ( trueTrackMPx < 0 && trueTrackMPy > 0 ) trueTrackMPhi += pi;
1436 float trueTrackMP = sqrt( trueTrackMPx*trueTrackMPx + trueTrackMPy*trueTrackMPy + trueTrackMPz*trueTrackMPz );
1437 float trueTrackMEta = 0.5*log( (trueTrackMP + trueTrackMPz)/(trueTrackMP - trueTrackMPz) );
1438
1439 float RcutMu = 100;
1440 float RcutTrack = 100;
1441 float RcutVtx = 100;
1442
1443 truthMatch.clear(); truthLambda.clear(); truthJpsi.clear();
1444
1445 //loop to check all LambdaB candidates found
1446 for (uint i=0; i<mupPx->size(); ++i){
1447
1448 bool istrueMup = false;
1449 bool istrueMum = false;
1450 bool istrueTrackP = false;
1451 bool istrueTrackM = false;
1452 bool istrueLambda = false;
1453 bool istrueJpsi = false;
1454 bool istrueLambdaB = false;
1455
1456 //calculate eta and phi for all tracks in B candidate
1457 float mupPhi = atan(mupPy->at(i)/mupPx->at(i));
1458 if( mupPx->at(i)<0 && mupPy->at(i)<0 ) mupPhi -=pi;
1459 if( mupPx->at(i)<0 && mupPy->at(i)>0 ) mupPhi +=pi;
1460 float mupP = sqrt( mupPx->at(i)*mupPx->at(i) + mupPy->at(i)*mupPy->at(i) + mupPz->at(i)*mupPz->at(i) );
1461 float mupEta = 0.5*log( (mupP + mupPz->at(i))/(mupP - mupPz->at(i)) );
1462
1463 float mumPhi = atan(mumPy->at(i)/mumPx->at(i));
1464 if ( mumPx->at(i) < 0 && mumPy->at(i) < 0 ) mumPhi -= pi;
1465 if ( mumPx->at(i) < 0 && mumPy->at(i) > 0 ) mumPhi += pi;
1466 float mumP = sqrt( mumPx->at(i)*mumPx->at(i) + mumPy->at(i)*mumPy->at(i) + mumPz->at(i)*mumPz->at(i) );
1467 float mumEta = 0.5*log( (mumP + mumPz->at(i))/(mumP - mumPz->at(i)) );
1468
1469 float trackPPhi = atan(trackPPy->at(i)/trackPPx->at(i));
1470 if ( trackPPx->at(i) < 0 && trackPPy->at(i) < 0 ) trackPPhi -= pi;
1471 if ( trackPPx->at(i) < 0 && trackPPy->at(i) > 0 ) trackPPhi += pi;
1472 float trackPP = sqrt( trackPPx->at(i)*trackPPx->at(i) + trackPPy->at(i)*trackPPy->at(i) + trackPPz->at(i)*trackPPz->at(i) );
1473 float trackPEta = 0.5*log( (trackPP + trackPPz->at(i))/(trackPP - trackPPz->at(i)) );
1474
1475 float trackMPhi = atan(trackMPy->at(i)/trackMPx->at(i));
1476 if ( trackMPx->at(i) < 0 && trackMPy->at(i) < 0 ) trackMPhi -= pi;
1477 if ( trackMPx->at(i) < 0 && trackMPy->at(i) > 0 ) trackMPhi += pi;
1478 float trackMP = sqrt( trackMPx->at(i)*trackMPx->at(i) + trackMPy->at(i)*trackMPy->at(i) + trackMPz->at(i)*trackMPz->at(i) );
1479 float trackMEta = 0.5*log( (trackMP + trackMPz->at(i))/(trackMP - trackMPz->at(i)) );
1480
1481
1482
1483 float deltaRmup = sqrt( (mupEta-trueMupEta)*(mupEta-trueMupEta) +
1484 (mupPhi-trueMupPhi)*(mupPhi-trueMupPhi) );
1485 if(deltaRmup<RcutMu)istrueMup = true;
1486 float deltaRmum = sqrt( (mumEta-trueMumEta)*(mumEta-trueMumEta) +
1487 (mumPhi-trueMumPhi)*(mumPhi-trueMumPhi) );
1488 if(deltaRmum<RcutMu)istrueMum = true;
1489
1490 float deltaRtrackp =sqrt( (trackPEta-trueTrackPEta)*(trackPEta-trueTrackPEta) +
1491 (trackPPhi-trueTrackPPhi)*(trackPPhi-trueTrackPPhi) );
1492 if(deltaRtrackp<RcutTrack)istrueTrackP = true;
1493
1494 float deltaRtrackm =sqrt( (trackMEta-trueTrackMEta)*(trackMEta-trueTrackMEta) +
1495 (trackMPhi-trueTrackMPhi)*(trackMPhi-trueTrackMPhi) );
1496 if(deltaRtrackm<RcutTrack)istrueTrackM = true;
1497
1498 float deltaRLambdaVtx = sqrt( (trueLambdaDecayVtxX-LambdaVtxX->at(i))*(trueLambdaDecayVtxX-LambdaVtxX->at(i)) +
1499 (trueLambdaDecayVtxY-LambdaVtxY->at(i))*(trueLambdaDecayVtxY-LambdaVtxY->at(i)) +
1500 (trueLambdaDecayVtxZ-LambdaVtxZ->at(i))*(trueLambdaDecayVtxZ-LambdaVtxZ->at(i)) );
1501
1502 dRmup->push_back(deltaRmup);
1503 dRmum->push_back(deltaRmum);
1504 dRtrackp->push_back(deltaRtrackp);
1505 dRtrackm->push_back(deltaRtrackm);
1506 dRLambdaVtx->push_back(deltaRLambdaVtx);
1507
1508 if(istrueMup && istrueMum)istrueJpsi = true;
1509 if(istrueTrackP && istrueTrackM && (deltaRLambdaVtx<RcutVtx) )istrueLambda = true;
1510 if(istrueJpsi && istrueLambda)istrueLambdaB = true;
1511
1512 if(istrueJpsi){
1513 truthJpsi.push_back(1);
1514 }else{
1515 truthJpsi.push_back(-1);
1516 }
1517 if(istrueLambda){
1518 truthLambda.push_back(1);
1519 }else{
1520 truthLambda.push_back(-1);
1521 }
1522 if(istrueLambdaB){
1523 truthMatch.push_back(1);
1524 }else{
1525 truthMatch.push_back(-1);
1526 }
1527
1528
1529 //cout<<"istrueMup = "<<istrueMup<<", istrueMum = "<<istrueMum<<", istrueTrackP = "<<istrueTrackP<<", istrueTrackM = "<<istrueTrackM<<endl;
1530 //cout<<"istrueJpsi = "<<istrueJpsi<<", istrueLambda = "<<istrueLambda<<", istrueLambdaB = "<<istrueLambdaB<<endl;
1531 //cout<<"trueLambdaBMass = "<<trueLambdaBMass<<endl;
1532
1533 }//end of loop to check the reco Lambda B with MC
1534 }//end of if the gen Lambdab is correct
1535 }//end of if Lambda^0_b
1536 }//end of loop on genParticles
1537 ///
1538 ///////////////////////////////
1539
1540 //fill truthMatch with zeros if we didn't fidn a truth match
1541 if ( truthMatch.size()==0 ) { //if no signal event has been found, fill with zeros
1542 for (uint i = 0; i<mupPx->size(); i++) {
1543 truthJpsi.push_back(0);
1544 truthLambda.push_back(0);
1545 truthMatch.push_back(0);
1546 }
1547 }
1548
1549
1550 }//end of if nLambdaB candidates > 0
1551
1552 //fill tree and clear vectors
1553 if(nLambdaB > 0){
1554 tree_->Fill();
1555 }
1556 nLambdaB = 0;
1557 run = 0;
1558 lumiBlock = 0;
1559 event = 0;
1560
1561 priVtxX->clear();
1562 priVtxY->clear();
1563 priVtxZ->clear();
1564 priVtxXE->clear();
1565 priVtxYE->clear();
1566 priVtxZE->clear();
1567 priVtxCL->clear();
1568 priVtxBSx->clear();
1569 priVtxBSy->clear();
1570 priVtxBSz->clear();
1571 priVtxBSxE->clear();
1572 priVtxBSyE->clear();
1573 priVtxBSzE->clear();
1574 priVtxBSCL->clear()
1575 ;
1576 beamSpotX = 0;
1577 beamSpotY = 0;
1578 beamSpotZ = 0;
1579 hlt_mu3 = 0;
1580 hlt_mu5 = 0;
1581 hlt_mu7 = 0;
1582 hlt_mu9 = 0;
1583 hlt_2mu0 = 0;
1584 hlt_2mu3 = 0;
1585 hlt_2mu0L2 = 0;
1586 hlt_2mu3JPsi = 0;
1587 hlt_BJPsiMuMu = 0;
1588 hlt_mu0trk0 = 0;
1589 hlt_mu3trk0 = 0;
1590 hlt_mu0trkmu0 = 0;
1591 hlt_mu3trkmu0 = 0;
1592 hlt_L1muOpen = 0;
1593 hlt_L12muOpen = 0;
1594 hlt_L12muOpenTight = 0;
1595 hlt_mu0_L1muOpen= 0;
1596 hlt_mu3_L1muOpen= 0;
1597 hlt_mu5_L1muOpen= 0;
1598 hlt_mu5trk0 = 0;
1599 hlt_mu0trkmu0_nocharge = 0;
1600 hlt_mu0trkmu0_OST = 0;
1601 hlt_2mu0_Qv1 = 0;
1602 hlt_2mu0_QLSv1 = 0;
1603 hlt_mu0trkmu0_OST_tightV1 = 0;
1604 hlt_mu0trkmu0_OST_tightV2 = 0;
1605 hlt_mu0trkmu0_OST_tightV3 = 0;
1606 JpsiVtxX->clear();
1607 JpsiVtxY->clear();
1608 JpsiVtxZ->clear();
1609 LambdaVtxX->clear();
1610 LambdaVtxY->clear();
1611 LambdaVtxZ->clear();
1612 JpsiVtxProb->clear();
1613 JpsiVtxCL->clear();
1614 LambdaVtxProb->clear();
1615 LambdaVtxCL->clear();
1616 LambdaBMass->clear();
1617 LambdaBPx->clear();
1618 LambdaBPy->clear();
1619 LambdaBPz->clear();
1620 JpsiPx->clear();
1621 JpsiPy->clear();
1622 JpsiPz->clear();
1623 LambdaPx->clear();
1624 LambdaPy->clear();
1625 LambdaPz->clear();
1626 JpsiMass->clear();
1627 LambdaMass->clear();
1628 mupPx->clear();
1629 mupPy->clear();
1630 mupPz->clear();
1631 mumPx->clear();
1632 mumPy->clear();
1633 mumPz->clear();
1634 trackPPx->clear();
1635 trackPPy->clear();
1636 trackPPz->clear();
1637 trackMPx->clear();
1638 trackMPy->clear();
1639 trackMPz->clear();
1640 JpsiLambdaVtxDistXY->clear();
1641 JpsiLambdaVtxDist3D->clear();
1642 openingAngleLambdaB->clear();
1643 openingAngleLambda->clear();
1644 LxyPV->clear();
1645 LxyBS->clear();
1646 L3DPV->clear();
1647 L3DBS->clear();
1648 ctauLambdaBxyPV->clear();
1649 ctauLambdaBxyBS->clear();
1650 ctauLambdaB3DPV->clear();
1651 ctauLambdaB3DBS->clear();
1652 dRmup->clear();
1653 dRmum->clear();
1654 dRtrackp->clear();
1655 dRtrackm->clear();
1656 dRLambdaVtx->clear();
1657 mupChi2->clear();
1658 mumChi2->clear();
1659 mupNDOF->clear();
1660 mumNDOF->clear();
1661 mupD0->clear();
1662 mumD0->clear();
1663 mupNvalidHits->clear();
1664 mumNvalidHits->clear();
1665 mupIsGlobal->clear();
1666 mumIsGlobal->clear();
1667 mupIsTracker->clear();
1668 mumIsTracker->clear();
1669 ctauLambdaXY->clear();
1670
1671 JpsiVtxX_afterKinFit->clear();
1672 JpsiVtxY_afterKinFit->clear();
1673 JpsiVtxZ_afterKinFit->clear();
1674 LambdaVtxX_afterKinFit->clear();
1675 LambdaVtxY_afterKinFit->clear();
1676 LambdaVtxZ_afterKinFit->clear();
1677 JpsiVtxProb_afterKinFit->clear();
1678 JpsiVtxCL_afterKinFit->clear();
1679 LambdaVtxProb_afterKinFit->clear();
1680 LambdaVtxCL_afterKinFit->clear();
1681 LambdaBMass_afterKinFit->clear();
1682 LambdaBPx_afterKinFit->clear();
1683 LambdaBPy_afterKinFit->clear();
1684 LambdaBPz_afterKinFit->clear();
1685 JpsiPx_afterKinFit->clear();
1686 JpsiPy_afterKinFit->clear();
1687 JpsiPz_afterKinFit->clear();
1688 LambdaPx_afterKinFit->clear();
1689 LambdaPy_afterKinFit->clear();
1690 LambdaPz_afterKinFit->clear();
1691 JpsiMass_afterKinFit->clear();
1692 LambdaMass_afterKinFit->clear();
1693 mupPx_afterKinFit->clear();
1694 mupPy_afterKinFit->clear();
1695 mupPz_afterKinFit->clear();
1696 mumPx_afterKinFit->clear();
1697 mumPy_afterKinFit->clear();
1698 mumPz_afterKinFit->clear();
1699 trackPPx_afterKinFit->clear();
1700 trackPPy_afterKinFit->clear();
1701 trackPPz_afterKinFit->clear();
1702 trackMPx_afterKinFit->clear();
1703 trackMPy_afterKinFit->clear();
1704 trackMPz_afterKinFit->clear();
1705 JpsiLambdaVtxDistXY_afterKinFit->clear();
1706 JpsiLambdaVtxDist3D_afterKinFit->clear();
1707 openingAngleLambdaB_afterKinFit->clear();
1708 openingAngleLambda_afterKinFit->clear();
1709 LxyPV_afterKinFit->clear();
1710 LxyBS_afterKinFit->clear();
1711 L3DPV_afterKinFit->clear();
1712 L3DBS_afterKinFit->clear();
1713 ctauLambdaBxyPV_afterKinFit->clear();
1714 ctauLambdaBxyBS_afterKinFit->clear();
1715 ctauLambdaB3DPV_afterKinFit->clear();
1716 ctauLambdaB3DBS_afterKinFit->clear();
1717 dRmup_afterKinFit->clear();
1718 dRmum_afterKinFit->clear();
1719 dRtrackp_afterKinFit->clear();
1720 dRtrackm_afterKinFit->clear();
1721 dRLambdaVtx_afterKinFit->clear();
1722 mupChi2_afterKinFit->clear();
1723 mumChi2_afterKinFit->clear();
1724 mupNDOF_afterKinFit->clear();
1725 mumNDOF_afterKinFit->clear();
1726 mupD0_afterKinFit->clear();
1727 mumD0_afterKinFit->clear();
1728 mupNvalidHits_afterKinFit->clear();
1729 mumNvalidHits_afterKinFit->clear();
1730 mupIsGlobal_afterKinFit->clear();
1731 mumIsGlobal_afterKinFit->clear();
1732 mupIsTracker_afterKinFit->clear();
1733 mumIsTracker_afterKinFit->clear();
1734 ctauLambdaXY_afterKinFit->clear();
1735 LambdaBX_afterKinFit->clear();
1736 LambdaBY_afterKinFit->clear();
1737 LambdaBZ_afterKinFit->clear();
1738 LambdaBVtxX_afterKinFit->clear();
1739 LambdaBVtxY_afterKinFit->clear();
1740 LambdaBVtxZ_afterKinFit->clear();
1741 LambdaBNDOF_afterKinFit->clear();
1742 LambdaBchiSquared_afterKinFit->clear();
1743 LambdaBVtxProb_afterKinFit->clear();
1744 LambdaBPt_afterKinFit->clear();
1745 LambdaBEta_afterKinFit->clear();
1746 LambdaBPhi_afterKinFit->clear();
1747
1748
1749 genLambdaJpsi = 0;
1750
1751
1752
1753 //MC truth
1754 truePriVtxX = 0;
1755 truePriVtxY = 0;
1756 truePriVtxZ = 0;
1757 trueLambdaBPx = 0;
1758 trueLambdaBPy = 0;
1759 trueLambdaBPz = 0;
1760 trueLambdaBDecayVtxX = 0;
1761 trueLambdaBDecayVtxY = 0;
1762 trueLambdaBDecayVtxZ = 0;
1763 trueJpsiPx = 0;
1764 trueJpsiPy = 0;
1765 trueJpsiPz = 0;
1766 trueJpsiDecayVtxX = 0;
1767 trueJpsiDecayVtxY = 0;
1768 trueJpsiDecayVtxZ = 0;
1769 trueMumPx = 0;
1770 trueMumPy = 0;
1771 trueMumPz = 0;
1772 trueMupPx = 0;
1773 trueMupPy = 0;
1774 trueMupPz = 0;
1775 trueLambdaPx = 0;
1776 trueLambdaPy = 0;
1777 trueLambdaPz = 0;
1778 trueLambdaDecayVtxX = 0;
1779 trueLambdaDecayVtxY = 0;
1780 trueLambdaDecayVtxZ = 0;
1781 trueTrackPPx = 0;
1782 trueTrackPPy = 0;
1783 trueTrackPPz = 0;
1784 trueTrackMPx = 0;
1785 trueTrackMPy = 0;
1786 trueTrackMPz = 0;
1787 trueLxyPV = 0;
1788 trueLxyBS = 0;
1789 trueL3DPV = 0;
1790 trueL3DBS = 0;
1791 truectauLambdaBxyPV = 0;
1792 truectauLambdaB3DPV = 0;
1793 trueLambdaBMass = 0;
1794
1795 truthMatch.clear();
1796 truthLambda.clear();
1797 truthJpsi.clear();
1798
1799
1800
1801
1802
1803
1804
1805 #ifdef THIS_IS_AN_EVENT_EXAMPLE
1806 Handle<ExampleData> pIn;
1807 iEvent.getByLabel("example",pIn);
1808 #endif
1809
1810 #ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
1811 ESHandle<SetupData> pSetup;
1812 iSetup.get<SetupRecord>().get(pSetup);
1813 #endif
1814 }
1815
1816
1817 // ------------ method called once each job just before starting event loop ------------
1818 void
1819 LambdaAnalyzer::beginJob()
1820 {
1821 tmpNvtx = 0;
1822 cout<<"begin job"<<endl;
1823 edm::Service<TFileService> fs;
1824 tree_ = fs->make<TTree>("ntpl","LambdaB ntuple");
1825
1826 //branches for reco particles
1827 tree_->Branch("nLambdaB",&nLambdaB,"nBLambda/i");
1828 tree_->Branch("run",&run,"run/f");
1829 tree_->Branch("lumiBlock",&lumiBlock,"lumiBlock/f");
1830 tree_->Branch("event",&event,"event/f");
1831
1832 tree_->Branch("priVtxX",&priVtxX);
1833 tree_->Branch("priVtxY",&priVtxY);
1834 tree_->Branch("priVtxZ",&priVtxZ);
1835 tree_->Branch("priVtxXE",&priVtxXE);
1836 tree_->Branch("priVtxYE",&priVtxYE);
1837 tree_->Branch("priVtxZE",&priVtxZE);
1838 tree_->Branch("priVtxCL",&priVtxCL);
1839 tree_->Branch("priVtxBSx",&priVtxBSx);
1840 tree_->Branch("priVtxBSy",&priVtxBSy);
1841 tree_->Branch("priVtxBSz",&priVtxBSz);
1842 tree_->Branch("priVtxBSxE",&priVtxBSxE);
1843 tree_->Branch("priVtxBSyE",&priVtxBSyE);
1844 tree_->Branch("priVtxBSzE",&priVtxBSzE);
1845 tree_->Branch("priVtxBSCL",&priVtxBSCL);
1846 tree_->Branch("beamSpotX",&beamSpotX, "beamSpotX/f");
1847 tree_->Branch("beamSpotY",&beamSpotY, "beamSpotY/f");
1848 tree_->Branch("beamSpotZ",&beamSpotZ, "beamSpotZ/f");
1849 tree_->Branch("hlt_mu3",&hlt_mu3,"hlt_mu3/i");
1850 tree_->Branch("hlt_mu5",&hlt_mu5,"hlt_mu5/i");
1851 tree_->Branch("hlt_mu7",&hlt_mu7,"hlt_mu7/i");
1852 tree_->Branch("hlt_mu9",&hlt_mu9,"hlt_mu9/i");
1853 tree_->Branch("hlt_2mu0",&hlt_2mu0,"hlt_2mu0/i");
1854 tree_->Branch("hlt_2mu3",&hlt_2mu3,"hlt_2mu3/i");
1855 tree_->Branch("hlt_2mu0L2",&hlt_2mu0L2,"hlt_2mu0L2/i");
1856 tree_->Branch("hlt_2mu3JPsi",&hlt_2mu3JPsi,"hlt_2mu3JPsi/i");
1857 tree_->Branch("hlt_BJPsiMuMu",&hlt_BJPsiMuMu,"hlt_BJPsiMuMu/i");
1858 tree_->Branch("hlt_mu0trk0",&hlt_mu0trk0,"hlt_mu0trk0/i");
1859 tree_->Branch("hlt_mu3trk0",&hlt_mu3trk0,"hlt_mu3trk0/i");
1860 tree_->Branch("hlt_mu0trkmu0",&hlt_mu0trkmu0,"hlt_mu0trkmu0/i");
1861 tree_->Branch("hlt_mu3trkmu0",&hlt_mu3trkmu0,"hlt_mu3trkmu0/i");
1862 tree_->Branch("hlt_L1muOpen",&hlt_L1muOpen,"hlt_L1muOpen/i");
1863 tree_->Branch("hlt_L12muOpen",&hlt_L12muOpen,"hlt_L12muOpen/i");
1864 tree_->Branch("hlt_L12muOpenTight",&hlt_L12muOpenTight,"hlt_L12muOpenTight/i");
1865 tree_->Branch("hlt_mu0_L1muOpen",&hlt_mu0_L1muOpen,"hlt_mu0_L1muOpen/i");
1866 tree_->Branch("hlt_mu3_L1muOpen",&hlt_mu3_L1muOpen,"hlt_mu3_L1muOpen/i");
1867 tree_->Branch("hlt_mu5_L1muOpen",&hlt_mu5_L1muOpen,"hlt_mu5_L1muOpen/i");
1868 tree_->Branch("hlt_mu5trk0",&hlt_mu5trk0,"hlt_mu5trk0/i");
1869 tree_->Branch("hlt_mu0trkmu0_nocharge",&hlt_mu0trkmu0_nocharge,"hlt_mu0trkmu0_nocharge/i");
1870 tree_->Branch("hlt_mu0trkmu0_OST",&hlt_mu0trkmu0_OST,"hlt_mu0trkmu0_OST/i");
1871 tree_->Branch("hlt_2mu0_Qv1",&hlt_2mu0_Qv1,"hlt_2mu0_Qv1/i");
1872 tree_->Branch("hlt_2mu0_QLSv1",&hlt_2mu0_QLSv1,"hlt_2mu0_QLSv1/i");
1873 tree_->Branch("hlt_mu0trkmu0_OST_tightV1",&hlt_mu0trkmu0_OST_tightV1,"hlt_mu0trkmu0_OST_tightV1/i");
1874 tree_->Branch("hlt_mu0trkmu0_OST_tightV2",&hlt_mu0trkmu0_OST_tightV2,"hlt_mu0trkmu0_OST_tightV2/i");
1875 tree_->Branch("hlt_mu0trkmu0_OST_tightV3",&hlt_mu0trkmu0_OST_tightV3,"hlt_mu0trkmu0_OST_tightV3/i");
1876 tree_->Branch("JpsiVtxX",&JpsiVtxX);
1877 tree_->Branch("JpsiVtxY",&JpsiVtxY);
1878 tree_->Branch("JpsiVtxZ",&JpsiVtxZ);
1879 tree_->Branch("LambdaVtxX",&LambdaVtxX);
1880 tree_->Branch("LambdaVtxY",&LambdaVtxY);
1881 tree_->Branch("LambdaVtxZ",&LambdaVtxZ);
1882 tree_->Branch("JpsiVtxProb",&JpsiVtxProb);
1883 tree_->Branch("JpsiVtxCL",&JpsiVtxCL);
1884 tree_->Branch("LambdaVtxProb",&LambdaVtxProb);
1885 tree_->Branch("LambdaVtxCL",&LambdaVtxCL);
1886 tree_->Branch("LambdaBMass",&LambdaBMass);
1887 tree_->Branch("LambdaBPx",&LambdaBPx);
1888 tree_->Branch("LambdaBPy",&LambdaBPy);
1889 tree_->Branch("LambdaBPz",&LambdaBPz);
1890 tree_->Branch("JpsiPx",&JpsiPx);
1891 tree_->Branch("JpsiPy",&JpsiPy);
1892 tree_->Branch("JpsiPz",&JpsiPz);
1893 tree_->Branch("LambdaPx",&LambdaPx);
1894 tree_->Branch("LambdaPy",&LambdaPy);
1895 tree_->Branch("LambdaPz",&LambdaPz);
1896 tree_->Branch("JpsiMass",&JpsiMass);
1897 tree_->Branch("LambdaMass",&LambdaMass);
1898 tree_->Branch("mupPx",&mupPx);
1899 tree_->Branch("mupPy",&mupPy);
1900 tree_->Branch("mupPz",&mupPz);
1901 tree_->Branch("mumPx",&mumPx);
1902 tree_->Branch("mumPy",&mumPy);
1903 tree_->Branch("mumPz",&mumPz);
1904 tree_->Branch("trackPPx",&trackPPx);
1905 tree_->Branch("trackPPy",&trackPPy);
1906 tree_->Branch("trackPPz",&trackPPz);
1907 tree_->Branch("trackMPx",&trackMPx);
1908 tree_->Branch("trackMPy",&trackMPy);
1909 tree_->Branch("trackMPz",&trackMPz);
1910 tree_->Branch("JpsiLambdaVtxDistXY",&JpsiLambdaVtxDistXY);
1911 tree_->Branch("JpsiLambdaVtxDist3D",&JpsiLambdaVtxDist3D);
1912 tree_->Branch("openingAngleLambdaB",&openingAngleLambdaB);
1913 tree_->Branch("openingAngleLambda",&openingAngleLambda);
1914 tree_->Branch("LxyPV",&LxyPV);
1915 tree_->Branch("LxyBS",&LxyBS);
1916 tree_->Branch("L3DPV",&L3DPV);
1917 tree_->Branch("L3DBS",&L3DBS);
1918 tree_->Branch("ctauLambdaBxyPV",&ctauLambdaBxyPV);
1919 tree_->Branch("ctauLambdaBxyBS",&ctauLambdaBxyBS);
1920 tree_->Branch("ctauLambdaB3DPV",&ctauLambdaB3DPV);
1921 tree_->Branch("ctauLambdaB3DBS",&ctauLambdaB3DBS);
1922 tree_->Branch("dRmup",&dRmup);
1923 tree_->Branch("dRmum",&dRmum);
1924 tree_->Branch("dRtrackp",&dRtrackp);
1925 tree_->Branch("dRtrackm",&dRtrackm);
1926 tree_->Branch("dRLambdaVtx",&dRLambdaVtx);
1927 tree_->Branch("mupChi2",&mupChi2);
1928 tree_->Branch("mumChi2",&mumChi2);
1929 tree_->Branch("mupNDOF",&mupNDOF);
1930 tree_->Branch("mumNDOF",&mumNDOF);
1931 tree_->Branch("mupD0",&mupD0);
1932 tree_->Branch("mumD0",&mumD0);
1933 tree_->Branch("mupNvalidHits",&mupNvalidHits);
1934 tree_->Branch("mumNvalidHits",&mumNvalidHits);
1935 tree_->Branch("mupIsGlobal",&mupIsGlobal);
1936 tree_->Branch("mumIsGlobal",&mumIsGlobal);
1937 tree_->Branch("mupIsTracker",&mupIsTracker);
1938 tree_->Branch("mumIsTracker",&mumIsTracker);
1939 tree_->Branch("ctauLambdaXY",&ctauLambdaXY);
1940
1941 //after kinematic fit
1942 tree_->Branch("JpsiVtxX_afterKinFit",&JpsiVtxX_afterKinFit);
1943 tree_->Branch("JpsiVtxY_afterKinFit",&JpsiVtxY_afterKinFit);
1944 tree_->Branch("JpsiVtxZ_afterKinFit",&JpsiVtxZ_afterKinFit);
1945 tree_->Branch("LambdaVtxX_afterKinFit",&LambdaVtxX_afterKinFit);
1946 tree_->Branch("LambdaVtxY_afterKinFit",&LambdaVtxY_afterKinFit);
1947 tree_->Branch("LambdaVtxZ_afterKinFit",&LambdaVtxZ_afterKinFit);
1948 tree_->Branch("JpsiVtxProb_afterKinFit",&JpsiVtxProb_afterKinFit);
1949 tree_->Branch("JpsiVtxCL_afterKinFit",&JpsiVtxCL_afterKinFit);
1950 tree_->Branch("LambdaVtxProb_afterKinFit",&LambdaVtxProb_afterKinFit);
1951 tree_->Branch("LambdaVtxCL_afterKinFit",&LambdaVtxCL_afterKinFit);
1952 tree_->Branch("LambdaBMass_afterKinFit",&LambdaBMass_afterKinFit);
1953 tree_->Branch("LambdaBPx_afterKinFit",&LambdaBPx_afterKinFit);
1954 tree_->Branch("LambdaBPy_afterKinFit",&LambdaBPy_afterKinFit);
1955 tree_->Branch("LambdaBPz_afterKinFit",&LambdaBPz_afterKinFit);
1956 tree_->Branch("JpsiPx_afterKinFit",&JpsiPx_afterKinFit);
1957 tree_->Branch("JpsiPy_afterKinFit",&JpsiPy_afterKinFit);
1958 tree_->Branch("JpsiPz_afterKinFit",&JpsiPz_afterKinFit);
1959 tree_->Branch("LambdaPx_afterKinFit",&LambdaPx_afterKinFit);
1960 tree_->Branch("LambdaPy_afterKinFit",&LambdaPy_afterKinFit);
1961 tree_->Branch("LambdaPz_afterKinFit",&LambdaPz_afterKinFit);
1962 tree_->Branch("JpsiMass_afterKinFit",&JpsiMass_afterKinFit);
1963 tree_->Branch("LambdaMass_afterKinFit",&LambdaMass_afterKinFit);
1964 tree_->Branch("mupPx_afterKinFit",&mupPx_afterKinFit);
1965 tree_->Branch("mupPy_afterKinFit",&mupPy_afterKinFit);
1966 tree_->Branch("mupPz_afterKinFit",&mupPz_afterKinFit);
1967 tree_->Branch("mumPx_afterKinFit",&mumPx_afterKinFit);
1968 tree_->Branch("mumPy_afterKinFit",&mumPy_afterKinFit);
1969 tree_->Branch("mumPz_afterKinFit",&mumPz_afterKinFit);
1970 tree_->Branch("trackPPx_afterKinFit",&trackPPx_afterKinFit);
1971 tree_->Branch("trackPPy_afterKinFit",&trackPPy_afterKinFit);
1972 tree_->Branch("trackPPz_afterKinFit",&trackPPz_afterKinFit);
1973 tree_->Branch("trackMPx_afterKinFit",&trackMPx_afterKinFit);
1974 tree_->Branch("trackMPy_afterKinFit",&trackMPy_afterKinFit);
1975 tree_->Branch("trackMPz_afterKinFit",&trackMPz_afterKinFit);
1976 tree_->Branch("JpsiLambdaVtxDistXY_afterKinFit",&JpsiLambdaVtxDistXY_afterKinFit);
1977 tree_->Branch("JpsiLambdaVtxDist3D_afterKinFit",&JpsiLambdaVtxDist3D_afterKinFit);
1978 tree_->Branch("openingAngleLambdaB_afterKinFit",&openingAngleLambdaB_afterKinFit);
1979 tree_->Branch("openingAngleLambda_afterKinFit",&openingAngleLambda_afterKinFit);
1980 tree_->Branch("LxyPV_afterKinFit",&LxyPV_afterKinFit);
1981 tree_->Branch("LxyBS_afterKinFit",&LxyBS_afterKinFit);
1982 tree_->Branch("L3DPV_afterKinFit",&L3DPV_afterKinFit);
1983 tree_->Branch("L3DBS_afterKinFit",&L3DBS_afterKinFit);
1984 tree_->Branch("ctauLambdaBxyPV_afterKinFit",&ctauLambdaBxyPV_afterKinFit);
1985 tree_->Branch("ctauLambdaBxyBS_afterKinFit",&ctauLambdaBxyBS_afterKinFit);
1986 tree_->Branch("ctauLambdaB3DPV_afterKinFit",&ctauLambdaB3DPV_afterKinFit);
1987 tree_->Branch("ctauLambdaB3DBS_afterKinFit",&ctauLambdaB3DBS_afterKinFit);
1988 tree_->Branch("dRmup_afterKinFit",&dRmup_afterKinFit);
1989 tree_->Branch("dRmum_afterKinFit",&dRmum_afterKinFit);
1990 tree_->Branch("dRtrackp_afterKinFit",&dRtrackp_afterKinFit);
1991 tree_->Branch("dRtrackm_afterKinFit",&dRtrackm_afterKinFit);
1992 tree_->Branch("dRLambdaVtx_afterKinFit",&dRLambdaVtx_afterKinFit);
1993 tree_->Branch("mupChi2_afterKinFit",&mupChi2_afterKinFit);
1994 tree_->Branch("mumChi2_afterKinFit",&mumChi2_afterKinFit);
1995 tree_->Branch("mupNDOF_afterKinFit",&mupNDOF_afterKinFit);
1996 tree_->Branch("mumNDOF_afterKinFit",&mumNDOF_afterKinFit);
1997 tree_->Branch("mupD0_afterKinFit",&mupD0_afterKinFit);
1998 tree_->Branch("mumD0_afterKinFit",&mumD0_afterKinFit);
1999 tree_->Branch("mupNvalidHits_afterKinFit",&mupNvalidHits_afterKinFit);
2000 tree_->Branch("mumNvalidHits_afterKinFit",&mumNvalidHits_afterKinFit);
2001 tree_->Branch("mupIsGlobal_afterKinFit",&mupIsGlobal_afterKinFit);
2002 tree_->Branch("mumIsGlobal_afterKinFit",&mumIsGlobal_afterKinFit);
2003 tree_->Branch("mupIsTracker_afterKinFit",&mupIsTracker_afterKinFit);
2004 tree_->Branch("mumIsTracker_afterKinFit",&mumIsTracker_afterKinFit);
2005 tree_->Branch("ctauLambdaXY_afterKinFit",&ctauLambdaXY_afterKinFit);
2006 tree_->Branch("LambdaBX_afterKinFit",&LambdaBX_afterKinFit);
2007 tree_->Branch("LambdaBY_afterKinFit",&LambdaBY_afterKinFit);
2008 tree_->Branch("LambdaBZ_afterKinFit",&LambdaBZ_afterKinFit);
2009 tree_->Branch("LambdaBVtxX_afterKinFit",&LambdaBVtxX_afterKinFit);
2010 tree_->Branch("LambdaBVtxY_afterKinFit",&LambdaBVtxY_afterKinFit);
2011 tree_->Branch("LambdaBVtxZ_afterKinFit",&LambdaBVtxZ_afterKinFit);
2012 tree_->Branch("LambdaBNDOF_afterKinFit",&LambdaBNDOF_afterKinFit);
2013 tree_->Branch("LambdaBchiSquared_afterKinFit",&LambdaBchiSquared_afterKinFit);
2014 tree_->Branch("LambdaBVtxProb_afterKinFit",&LambdaBVtxProb_afterKinFit);
2015 tree_->Branch("LambdaBPt_afterKinFit",&LambdaBPt_afterKinFit);
2016 tree_->Branch("LambdaBEta_afterKinFit",&LambdaBEta_afterKinFit);
2017 tree_->Branch("LambdaBPhi_afterKinFit",&LambdaBPhi_afterKinFit);
2018
2019 tree_->Branch("genLambdaJpsi", &genLambdaJpsi, "genLambdaJpsi/I");
2020
2021 // branches for MC truth
2022 tree_->Branch("truePriVtxX", &truePriVtxX, "truePriVtxX/f");
2023 tree_->Branch("truePriVtxY", &truePriVtxY, "truePriVtxY/f");
2024 tree_->Branch("truePriVtxZ", &truePriVtxZ, "truePriVtxZ/f");
2025 tree_->Branch("trueLambdaBPx",&trueLambdaBPx, "trueLambdaBPx/f");
2026 tree_->Branch("trueLambdaBPy",&trueLambdaBPy, "trueLambdaBPy/f");
2027 tree_->Branch("trueLambdaBPz",&trueLambdaBPz, "trueLambdaBPz/f");
2028 tree_->Branch("trueLambdaBDecayVtxX",&trueLambdaBDecayVtxX, "trueLambdaBDecayVtxX/f");
2029 tree_->Branch("trueLambdaBDecayVtxY",&trueLambdaBDecayVtxY, "trueLambdaBDecayVtxY/f");
2030 tree_->Branch("trueLambdaBDecayVtxZ",&trueLambdaBDecayVtxZ, "trueLambdaBDecayVtxZ/f");
2031 tree_->Branch("trueJpsiPx",&trueJpsiPx, "trueJpsiPx/f");
2032 tree_->Branch("trueJpsiPy",&trueJpsiPy, "trueJpsiPy/f");
2033 tree_->Branch("trueJpsiPz",&trueJpsiPz, "trueJpsiPz/f");
2034 tree_->Branch("trueJpsiDecayVtxX",&trueJpsiDecayVtxX, "trueJpsiDecayVtxX/f");
2035 tree_->Branch("trueJpsiDecayVtxY",&trueJpsiDecayVtxY, "trueJpsiDecayVtxY/f");
2036 tree_->Branch("trueJpsiDecayVtxZ",&trueJpsiDecayVtxZ, "trueJpsiDecayVtxZ/f");
2037 tree_->Branch("trueMumPx",&trueMumPx, "trueMumPx/f");
2038 tree_->Branch("trueMumPy",&trueMumPy, "trueMumPy/f");
2039 tree_->Branch("trueMumPz",&trueMumPz, "trueMumPz/f");
2040 tree_->Branch("trueMupPx",&trueMupPx, "trueMupPx/f");
2041 tree_->Branch("trueMupPy",&trueMupPy, "trueMupPy/f");
2042 tree_->Branch("trueMupPz",&trueMupPz, "trueMupPz/f");
2043 tree_->Branch("trueLambdaPx",&trueLambdaPx, "trueLambdaPx/f");
2044 tree_->Branch("trueLambdaPy",&trueLambdaPy, "trueLambdaPy/f");
2045 tree_->Branch("trueLambdaPz",&trueLambdaPz, "trueLambdaPz/f");
2046 tree_->Branch("trueLambdaDecayVtxX",&trueLambdaDecayVtxX, "trueLambdaDecayVtxX/f");
2047 tree_->Branch("trueLambdaDecayVtxY",&trueLambdaDecayVtxY, "trueLambdaDecayVtxY/f");
2048 tree_->Branch("trueLambdaDecayVtxZ",&trueLambdaDecayVtxZ, "trueLambdaDecayVtxZ/f");
2049 tree_->Branch("trueTrackPPx",&trueTrackPPx, "trueTrackPPx/f");
2050 tree_->Branch("trueTrackPPy",&trueTrackPPy, "trueTrackPPy/f");
2051 tree_->Branch("trueTrackPPz",&trueTrackPPz, "trueTrackPPz/f");
2052 tree_->Branch("trueTrackMPx",&trueTrackMPx, "trueTrackMPx/f");
2053 tree_->Branch("trueTrackMPy",&trueTrackMPy, "trueTrackMPy/f");
2054 tree_->Branch("trueTrackMPz",&trueTrackMPz, "trueTrackMPz/f");
2055 tree_->Branch("trueLxyPV",&trueLxyPV, "trueLxyPV/f");
2056 tree_->Branch("trueLxyBS",&trueLxyBS, "trueLxyBS/f");
2057 tree_->Branch("trueL3DPV",&trueL3DPV, "trueL3DPV/f");
2058 tree_->Branch("trueL3DBS",&trueL3DBS, "trueL3DBS/f");
2059 tree_->Branch("truectauLambdaBxyPV",&truectauLambdaBxyPV, "truectauLambdaBxyPV/f");
2060 tree_->Branch("truectauLambdaB3DPV",&truectauLambdaB3DPV, "truectauLambdaB3DPV/f");
2061 tree_->Branch("trueLambdaBMass",&trueLambdaBMass, "trueLambdaBMass/f");
2062
2063 tree_->Branch("truthMatch",&truthMatch);
2064 tree_->Branch("truthLambda",&truthLambda);
2065 tree_->Branch("truthJpsi",&truthJpsi);
2066
2067
2068
2069 }
2070
2071 // ------------ method called once each job just after ending the event loop ------------
2072 void
2073 LambdaAnalyzer::endJob() {
2074 cout<<"number of valid vertices: "<<tmpNvtx<<endl;
2075 // tree_->GetDirectory()->cd();
2076 // tree_->Write();
2077 }
2078
2079
2080
2081 bool LambdaAnalyzer::pass_HLT_Mu0_TkMu0_Jpsi(const edm::Event& iEvent){
2082
2083 edm::Handle<trigger::TriggerEvent> hltevent;
2084 iEvent.getByLabel(triggerEventTag_, hltevent);
2085 const trigger::TriggerObjectCollection&
2086 hltobjects(hltevent->getObjects());
2087
2088 trigger::size_type index(0);
2089 edm::InputTag it_track("hltMu0TrackJpsiTrackMassFiltered::REDIGI36X");
2090
2091 index=hltevent->filterIndex(it_track);
2092 if (index<hltevent->sizeFilters()) {
2093 const trigger::Keys& KEYS (hltevent->filterKeys(index));
2094 const trigger::size_type nTO(KEYS.size());
2095 for (trigger::size_type i=0; i!=nTO; ++i) {
2096 const trigger::TriggerObject& TO1( hltobjects.at(KEYS[i]) );
2097 double pt1 = TO1.pt(); double eta1 = TO1.eta(); double phi1 = TO1.phi();
2098 i++;
2099 const trigger::TriggerObject& TO2( hltobjects.at(KEYS[i]) );
2100 double pt2 = TO2.pt(); double eta2 = TO2.eta(); double phi2 = TO2.phi();
2101 double massJpsi = sqrt(2*pt1*pt2*(TMath::CosH(eta1)*TMath::CosH(eta2)-TMath::SinH(eta1)*TMath::SinH(eta2)-TMath::Sin(phi1)*TMath::Sin(phi2)-TMath::Cos(phi1)*TMath::Cos(phi2)));
2102 if (massJpsi>2.5 && massJpsi<4.1) return true;
2103 }
2104 }
2105
2106 return false;
2107 }
2108
2109 bool LambdaAnalyzer::pass_HLT_Mu0_TkMu0_OST_Jpsi(const edm::Event& iEvent) {
2110 return pass_HLT_Mu0_TkMu0_Jpsi(iEvent);
2111 }
2112
2113 bool LambdaAnalyzer::pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v1(const edm::Event& iEvent)
2114 {
2115
2116 edm::Handle<reco::MuonCollection> theMuons;
2117 iEvent.getByLabel( edm::InputTag("muons"), theMuons);
2118
2119 edm::Handle<trigger::TriggerEvent> hltevent;
2120 iEvent.getByLabel(triggerEventTag_, hltevent);
2121 const trigger::TriggerObjectCollection&
2122 hltobjects(hltevent->getObjects());
2123
2124 trigger::size_type index(0);
2125 edm::InputTag it_track("hltMu0TrackJpsiTrackMassFiltered::REDIGI36X");
2126
2127 index=hltevent->filterIndex(it_track);
2128 if (index<hltevent->sizeFilters()) {
2129 const trigger::Keys& KEYS (hltevent->filterKeys(index));
2130 const trigger::size_type nTO(KEYS.size());
2131 for (trigger::size_type i=0; i!=nTO; ++i) {
2132 const trigger::TriggerObject& TO1( hltobjects.at(KEYS[i]) );
2133 double pt1 = TO1.pt(); double eta1 = TO1.eta(); double phi1 = TO1.phi();
2134 i++;
2135 const trigger::TriggerObject& TO2( hltobjects.at(KEYS[i]) );
2136 double pt2 = TO2.pt(); double eta2 = TO2.eta(); double phi2 = TO2.phi();
2137 double minDR1 = 999.;
2138 int q1 = 0;
2139 for(reco::MuonCollection::const_iterator iMuon=theMuons->begin(); iMuon!=theMuons->end(); ++iMuon){
2140 double dR = deltaR(eta1,phi1,iMuon->eta(),iMuon->phi());
2141 if(dR<minDR1){
2142 minDR1 = dR;
2143 q1 = iMuon->charge();
2144 }
2145 }
2146 double minDR2 = 999.;
2147 int q2 = 0;
2148 for(reco::MuonCollection::const_iterator iMuon=theMuons->begin(); iMuon!=theMuons->end(); ++iMuon){
2149 double dR = deltaR(eta2,phi2,iMuon->eta(),iMuon->phi());
2150 if(dR<minDR2){
2151 minDR2 = dR;
2152 q2 = iMuon->charge();
2153 }
2154 }
2155 double massJpsi = sqrt(2*pt1*pt2*(TMath::CosH(eta1)*TMath::CosH(eta2)-TMath::SinH(eta1)*TMath::SinH(eta2)-TMath::Sin(phi1)*TMath::Sin(phi2)-TMath::Cos(phi1)*TMath::Cos(phi2)));
2156 bool cowboy=false;
2157 if(minDR1<minDR2) cowboy = q1*deltaPhi(phi1,phi2) > 0;
2158 else cowboy = q2*deltaPhi(phi2,phi1) > 0;
2159 if (massJpsi>2.5 && massJpsi<4.1 && !cowboy) return true;
2160 }
2161 }
2162
2163
2164 return false;
2165 }
2166
2167 bool LambdaAnalyzer::pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v2(const
2168 edm::Event& iEvent)
2169 {
2170 return pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v1(iEvent);
2171 }
2172
2173 bool LambdaAnalyzer::pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v3(const
2174 edm::Event& iEvent)
2175 {
2176 return pass_HLT_Mu0_TkMu0_OST_Jpsi_Tight_v2(iEvent);
2177 }
2178
2179
2180
2181 bool LambdaAnalyzer::pass_HLT_DoubleMu0_Quarkonium_v1(const edm::Event& iEvent)
2182 {
2183
2184 edm::Handle<trigger::TriggerEvent> hltevent;
2185 iEvent.getByLabel(triggerEventTag_, hltevent);
2186 const trigger::TriggerObjectCollection& hltobjects(hltevent->getObjects());
2187
2188 trigger::size_type index(0);
2189 edm::InputTag it_track("hltDiMuonL3PreFiltered0::REDIGI36X");
2190
2191 index=hltevent->filterIndex(it_track);
2192 if (index<hltevent->sizeFilters()) {
2193 const trigger::Keys& KEYS (hltevent->filterKeys(index));
2194 const trigger::size_type nTO(KEYS.size());
2195 for (trigger::size_type i=0; i!=nTO; ++i) {
2196 const trigger::TriggerObject& TO1( hltobjects.at(KEYS[i]) );
2197 double pt1 = TO1.pt(); double eta1 = TO1.eta(); double phi1 = TO1.phi();
2198 for (trigger::size_type j=i+1; j!=nTO; ++j) {
2199 const trigger::TriggerObject& TO2( hltobjects.at(KEYS[j]) );
2200 double pt2 = TO2.pt(); double eta2 = TO2.eta(); double phi2 = TO2.phi();
2201 double massJpsi = sqrt(2*pt1*pt2*(TMath::CosH(eta1)*TMath::CosH(eta2)-TMath::SinH(eta1)*TMath::SinH(eta2)-TMath::Sin(phi1)*TMath::Sin(phi2)-TMath::Cos(phi1)*TMath::Cos(phi2)));
2202 if (massJpsi>1.5 && massJpsi<14.5) return true;
2203 }
2204 }
2205 }
2206
2207 return false;
2208 }
2209
2210
2211 void LambdaAnalyzer::push_anything_back()
2212 {
2213
2214 LambdaBMass_afterKinFit->push_back(-99999);
2215 LambdaBPx_afterKinFit->push_back(-99999);
2216 LambdaBPy_afterKinFit->push_back(-99999);
2217 LambdaBPz_afterKinFit->push_back(-99999);
2218 LambdaBX_afterKinFit->push_back(-99999);
2219 LambdaBY_afterKinFit->push_back(-99999);
2220 LambdaBZ_afterKinFit->push_back(-99999);
2221 LambdaBEta_afterKinFit->push_back(-99999);
2222 LambdaBPhi_afterKinFit->push_back(-99999);
2223 LambdaBPt_afterKinFit->push_back(-99999);
2224 LambdaBNDOF_afterKinFit->push_back(-99999);
2225 LambdaBchiSquared_afterKinFit->push_back(-99999);
2226 LambdaBVtxProb_afterKinFit->push_back(-99999);
2227 LambdaBVtxX_afterKinFit->push_back(-99999);
2228 LambdaBVtxY_afterKinFit->push_back(-99999);
2229 LambdaBVtxZ_afterKinFit->push_back(-99999);
2230
2231 mumPx_afterKinFit->push_back(-99999);
2232 mumPy_afterKinFit->push_back(-99999);
2233 mumPz_afterKinFit->push_back(-99999);
2234
2235 mupPx_afterKinFit->push_back(-99999);
2236 mupPy_afterKinFit->push_back(-99999);
2237 mupPz_afterKinFit->push_back(-99999);
2238
2239 JpsiPx_afterKinFit->push_back(-99999);
2240 JpsiPy_afterKinFit->push_back(-99999);
2241 JpsiPz_afterKinFit->push_back(-99999);
2242
2243 LambdaPx_afterKinFit->push_back(-99999);
2244 LambdaPy_afterKinFit->push_back(-99999);
2245 LambdaPz_afterKinFit->push_back(-99999);
2246
2247 JpsiVtxX_afterKinFit->push_back(-99999);
2248 JpsiVtxY_afterKinFit->push_back(-99999);
2249 JpsiVtxZ_afterKinFit->push_back(-99999);
2250
2251 LambdaVtxX_afterKinFit->push_back(-99999);
2252 LambdaVtxY_afterKinFit->push_back(-99999);
2253 LambdaVtxZ_afterKinFit->push_back(-99999);
2254
2255 LambdaVtxProb_afterKinFit->push_back(-99999);
2256 JpsiVtxProb_afterKinFit->push_back(-99999);
2257
2258 openingAngleLambdaB_afterKinFit->push_back(-99999);
2259 ctauLambdaB3DPV_afterKinFit->push_back(-99999);
2260 ctauLambdaBxyPV_afterKinFit->push_back(-99999);
2261 ctauLambdaXY_afterKinFit->push_back(-99999);
2262
2263 return;
2264
2265 }
2266
2267
2268 //define this as a plug-in
2269 DEFINE_FWK_MODULE(LambdaAnalyzer);