chargeCalculator.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef MDT_DIGITIZATION_CHARGECALCULATOR_H
#define MDT_DIGITIZATION_CHARGECALCULATOR_H
/*-----------------------------------------------
 
Created 20-8-2011 by Oleg.Bulekov@cern.ch
Function chargeCalculator returns the value of electric charge for Qball particle.
The information about charge coded in the pdgid of Qball particle is used.
It is supposed that the first digit of decimal number of charge corresponds
to the third digit of decimal number of the  pdgid.
(e.g. Qball pdgid=10000200 corresponds to Q=2 )
 
-----------------------------------------------*/
 
#include <iostream>
#include <sstream>
#include <vector>
 
#include "CLHEP/Units/PhysicalConstants.h"
#include "GeneratorObjects/HepMcParticleLink.h"
#include "MuonSimData/MuonSimData.h"
#include "MuonSimData/MuonSimDataCollection.h"
#include "MuonSimEvent/MDTSimHitCollection.h"
// SB
#include "AtlasHepMC/GenParticle.h"
//
double chargeCalculator(const EventContext& ctx, const MDTSimHit& hit, unsigned short eventId = 0) {
  const HepMcParticleLink trkParticle = HepMcParticleLink::getRedirectedLink(hit.particleLink(), eventId, ctx); // This link should now correctly resolve to the TruthEvent McEventCollection in the main StoreGateSvc.
  HepMC::ConstGenParticlePtr genParticle = trkParticle.cptr();
  double qcharge = 1.;
  if (genParticle) {
    const int particleEncoding = genParticle->pdg_id();
    if (((int)(std::abs(particleEncoding) / 10000000) == 1) && ((int)(std::abs(particleEncoding) / 100000) == 100)) { // TODO use a function from TruthUtils/AtlasPID.h
      qcharge = ((std::abs(particleEncoding) / 100000.0) - 100.0) * 1000.0;
      if (particleEncoding < 0.0) qcharge = -qcharge;
    }
    else if (((int)(std::abs(particleEncoding) / 10000000) == 2) && ((int)(std::abs(particleEncoding) / 100000) == 200)) { // TODO use a function from TruthUtils/AtlasPID.h
      qcharge = (double)((std::abs(particleEncoding) / 1000) % 100) / (double)((std::abs(particleEncoding) / 10) % 100);
      if (particleEncoding < 0.0) qcharge = -qcharge;
    }
  }
  return qcharge;
}
 
#endif