ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/codeHeavyBaryons/HeavyBaryonsAnalysis/Lambda/LambdaAnalyzer/src/SigmaAnalyzer.cc
Revision: 1.1.1.1 (vendor branch)
Committed: Mon Dec 13 15:58:16 2010 UTC (14 years, 5 months ago) by mirena
Content type: text/plain
Branch: LambdaBAnalyzer
CVS Tags: START
Changes since 1.1: +0 -0 lines
Log Message:

File Contents

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