DaVinci Tutorial 2

The purpose of these exercises is to allow you to write a complete though simple selection algorithms for a typical decay: Bs → J/ψφ. Let's do the J/ψ→μμ first. We start for the algorithm created in DaVinciTutorial1.

Add some global variables

  • We want to be able to cut on the J/psi mass and the vertex chi2 at least. And change the cut values by options.
  double m_jPsiMassWin ; ///< Mass window 
  double m_jPsiChi2 ; ///< Max J/psi chi^2

  • Let's also store a few things we will need in the execute()
LHCb::ParticleID m_jPsiID ; ///< J/psi ID 
double m_jPsiMass ; ///< J/psi mass 
  • It is mantatory to initialise everything in the constructor (in the .cpp file).
TutorialAlgorithm::TutorialAlgorithm( const std::string& name, ISvcLocator* pSvcLocator ) : DaVinciTupleAlgorithm( name, pSvcLocator ) 
  , m_jPsiID(0) 
  , m_jPsiMass(0.) 
  declareProperty("MassWindow", m_jPsiMassWin = 10.*Gaudi::Units::GeV);   
  declareProperty("MaxChi2", m_jPsiChi2 = 1000.); 

Note that while the compiler issues a warning when you use an uninitialized local variable in a method, it doesn't for uninitialized class members. People have run into horrible troubles impossible to debug because of such mistakes.

  • Of course, we want to get some real numbers in there. This happens in initialize()
debug() << "==> Initialize" << endmsg; 
const LHCb::ParticleProperty* Jpsi = ppSvc()->find( "J/psi(1S)" );
Now get the mass and the PID of the J/psi from the ParticleProperty and store these values in the private data members, for example:
m_jPsiID = LHCb::ParticleID(Jpsi->pdgID()); 
To find out the method you need to get the mass look at the doxygen documentation.

The list of official particle names and their properties is stored in a database. To dump the current list you can run in a new window:

lb-run DaVinci/latest dump_particle_properties | tee ParticleTable.txt
This will create the ParticleTable.txt file in the current directory, which lists all particles and their properties.

Make the J/ψ

  • First separate positive and negative muons. This is best done using the DaVinci::filter functions. Add a new private method StatusCode makeJpsi(const LHCb::Particle::ConstVector&) ; and implement it:
#include "Kernel/ParticleFilters.h"
#include <boost/lambda/lambda.hpp>
#include <boost/lambda/bind.hpp>

using namespace boost::lambda;  
LHCb::Particle::ConstVector MuPlus, MuMinus;
StatusCode sc = StatusCode::SUCCESS;

size_t nMuons = DaVinci::filter(muons, bind(&LHCb::Particle::charge,_1)<0, MuMinus);
if (nMuons>0) nMuons += DaVinci::filter(muons, bind(&LHCb::Particle::charge,_1)>0, MuPlus);
Assuming that the muons you got from the TES are in a vector called muons.

  • Loop over the muons, you already know how. You need two nested for loop, one for positive muons, one for negative muons.

  • Cut on the dimuon mass before you make the vertex, as the latter takes much more CPU. Inside the loop:
Gaudi::LorentzVector twoMu = (*imp)->momentum() + (*imm)->momentum();
if ( fabs(twoMu.M() - m_jPsiMass) > m_jPsiMassWin ) continue; // mass cuts
You probably also want to plot the mass before the cut, and maybe to print it out in debug() or verbose() level. imp and imm are the iterator on MuPlus and MuMinus

  • Now call the vertex fitter. It will return you a vertex and the J/ψ Particle.
LHCb::Vertex MuMuVertex; 
LHCb::Particle Jpsi( m_jPsiID );
StatusCode scFit = vertexFitter()->fit( *(*imp), *(*imm), MuMuVertex, Jpsi );
if ( !scFit ) {
  Warning("Fit error").ignore(); 
Presently the fitter returns a failure if the fit failed. This is not a good reason to stop the execution. So catch this error and handle it properly. The ignore() is just there because Warning() returns a StatusCode and all StatusCodes should be checked or explicitly ignored. Else you'll get a list of unchecked StatusCodes at the end of the logfile.

  • Cut on the chi2.
if ( MuMuVertex.chi2() > m_jPsiChi2 ) continue; // chi2 cut
But maybe you should fill a histogram with the chi2 first.

  • Decare the the J/ψ as to be kept.
setFilterPassed( true ); // Mandatory. Set to true if event is accepted
Warning: It does not save it yet.

Since you have a J/ψ you can set FilterPassed to true, which tells the sequencer the algorithm found what it was looking for.

Now you have only to add this function into the execute method:

sc = makeJpsi(muons);
  if (!sc) return sc;
remember to delete the setFilterPassed(true) statement, because it is evaluate in the makeJpsi function.

Need some counters?

Counters are booked on-demand like plots. Just add at an appropriate place
and you'll get a count of the number of candidates.

Run it

  • You need to set the appropriate options
# J/psi->mumu selection
jpsi2mumu = TutorialAlgorithm("Jpsi2MuMu");
jpsi2mumu.InputLocations = [ "Phys/StdAllLooseMuons" ]
from GaudiKernel.SystemOfUnits import MeV
jpsi2mumu.MassWindow = 30*MeV


Don't panic!

You'll see some of those:

LHCb::ParticleID m_jPsiID ; ///< J/psi ID 
Jpsi2MuMu                                                   WARNING TutorialAlgorithm:: Fit error StatusCode=FAILURE
Jpsi2MuMu                                                   WARNING TutorialAlgorithm:: Fit error StatusCode=FAILURE
Jpsi2MuMu                                                   WARNING TutorialAlgorithm:: Fit error StatusCode=FAILURE

That just means the vertex fit failed, which can happen.


The solution is given in solutions/DaVinci2.

Some plots

Try to write the macro that does the plots of the slides.


  • If you'd like to know how the same algorithm can select J/ψ→μμ and φ(1020)→KK go to DaVinciTutorial3.
  • If you'd like to use standard code to do the same thing go to DaVinciTutorial4.

-- PatrickKoppenburg - 01 Oct 2007 -- PatrickKoppenburg - 13 Jun 2008 -- PatrickKoppenburg - 05 Jan 2009 -- PatrickSKoppenburg - 16-Oct-2012 -- PatrickSKoppenburg - 30-Sep-2013

Topic attachments
I Attachment History Action Size Date Who Comment
PNGpng DM.png r1 manage 16.0 K 2012-10-16 - 15:36 PatrickSKoppenburg  
PNGpng DiMuM.png r1 manage 16.0 K 2012-10-16 - 15:36 PatrickSKoppenburg  
Edit | Attach | Watch | Print version | History: r29 < r28 < r27 < r26 < r25 | Backlinks | Raw View | WYSIWYG | More topic actions
Topic revision: r29 - 2017-12-19 - EduardoRodrigues
    • Cern Search Icon Cern Search
    • TWiki Search Icon TWiki Search
    • Google Search Icon Google Search

    LHCb All webs login

This site is powered by the TWiki collaboration platform Powered by PerlCopyright & 2008-2019 by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback