ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/codeHeavyBaryons/HeavyBaryonsAnalysis/LambdaBAnalyzer/src/LambdaBAnalyzer.cc
Revision: 1.1.1.1 (vendor branch)
Committed: Fri Dec 10 16:36:11 2010 UTC (14 years, 5 months ago) by mirena
Content type: text/plain
Branch: LambdaBAnalyzer, codeHeavyBaryons, HeavyBaryonsAnalysis
Changes since 1.1: +1 -2 lines
Log Message:
*** empty log message ***

File Contents

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