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

File Contents

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