ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/codeHeavyBaryons/HeavyBaryonsAnalysis/Lambda/LambdaAnalyzer/src/SigmaAnalyzer.cc
Revision: 1.1
Committed: Mon Dec 13 15:58:16 2010 UTC (14 years, 5 months ago) by mirena
Content type: text/plain
Branch: MAIN
Branch point for: LambdaBAnalyzer
Log Message:
Initial revision

File Contents

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