PndFsmMdtPid.cxx

Go to the documentation of this file.

//--------------------------------------------------------------------------
// File and Version Information:
//      $Id: PndFsmMdtPid.cxx $
//
// Description:
//      Class PndFsmMdtPid
//
//  Implementation of the MDT for the FastSim
//
//  This software was developed for the PANDA collaboration.  If you
//  use all or part of it, please give an appropriate acknowledgement.
//
// Author List:
//      Ralf Kliemt                     Original Author
//
// Copyright Information:
//      Copyright (C) 2014              GSI
//
//------------------------------------------------------------------------

//-----------------------
// This Class's Header --
//-----------------------
#include "PndFsmMdtPid.h"

//-------------
// C Headers --
//-------------

//---------------
// C++ Headers --
//---------------
#include <math.h>
#include <iostream>

using std::cout;
using std::endl;
using std::ostream;
using std::string;

//-------------------------------
// Collaborating Class Headers --
//-------------------------------

#include "ArgList.h"
#include "PndFsmResponse.h"
#include "PndFsmTrack.h"
#include "TH3F.h"
#include "TH1D.h"
#include "TFile.h"

//-----------------------------------------------------------------------
// Local Macros, Typedefs, Structures, Unions and Forward Declarations --
//-----------------------------------------------------------------------

//----------------
// Constructors --
//----------------

PndFsmMdtPid::PndFsmMdtPid()
{
  initParameters();

  _thtMin=_thtMin*M_PI/180.0;
  _thtMax=_thtMax*M_PI/180.0;
  readParameters();
  //print(std::cout);
}

PndFsmMdtPid::PndFsmMdtPid(ArgList &par)
{
  initParameters();
  //set default parameter values and parses a parameter list
  //i.e. std::list<std::string> of the form
  //"a=1" "b=2" "c=3"
  parseParameterList(par);

  _thtMin=_thtMin*M_PI/180.0;
  _thtMax=_thtMax*M_PI/180.0;
  readParameters();
  //print(std::cout);
}

//--------------
// Destructor --
//--------------

PndFsmMdtPid::~PndFsmMdtPid()
{
}

//--------------
// Operations --
//--------------

PndFsmResponse*
PndFsmMdtPid::respond(PndFsmTrack *t)
{
  PndFsmResponse *result=new PndFsmResponse();

  result->setDetector(this);
  bool wasDetected=detected(t);
  result->setDetected(wasDetected);

  if (wasDetected)
  {
          if (_useFlat)  // just use very simple likelihoods for muons and pions
          {
                        result->setLHElectron(0.);
                        result->setLHMuon( 1-_misId );
                        result->setLHPion( _misId );
                        result->setLHKaon( 0. );
                        result->setLHProton( 0. );
          }
          else   // use sophisticated methode with histogram lookup likelihoods
          {
                        int    type   = t->pdt();
                        double p      = t->p4().Vect().Mag();
                        double tht    = t->p4().Theta()*57.296;
                        double charge = t->charge();

                        // convert pdg code in type: 0=e+-, 1=mu+-, 2=pi+-, 3=K+-, 4=p, 5=p-bar
                        int idx;
                        if (abs(type)==11) idx=0;
                        else if (abs(type)==13) idx=1;
                        else if (abs(type)==211) idx=2;
                        else if (abs(type)==321) idx=3;
                        else if (type==2212) idx=4;
                        else idx=5;

                        // histogram max p = 8.0
                        if (p>7.99) p=7.99;

                        // get the bin corresponding to p
                        int currbinT = _mdtPidPdf[0]->GetYaxis()->FindBin(tht);
                        int currbinP = _mdtPidPdf[0]->GetZaxis()->FindBin(p);

                        // get the slice containing the pdf emcecal(p) for true particle type
                        TH1D *hpdf = _mdtPidPdf[idx]->ProjectionX("_tmppdf",currbinT, currbinT, currbinP, currbinP);

                        // get a random emcecal value from the true distribution
                        // if distribution is empty, choose 0
                        double xsig = 0.;
                        if (hpdf->Integral()>0) xsig = hpdf->GetRandom();

                        // store the value in the response object
                        result->setMuoIron(xsig);

                        // find the bin corresponding to the value
                        int xsigbin = hpdf->FindBin(xsig);

                        // get the probability values for the different particle types for this emcecal value; sum needed for normalization
                        double P[6], Psum=0.;

                        // sum up e ... K
                        for (int k=0;k<6;++k)
                        {
                                P[k] = _mdtPidPdf[k]->GetBinContent(xsigbin,currbinT,currbinP);
                                if (k<4) Psum += P[k];
                        }
                        // add P_p or P_pbar depending on charge of particle
                        if (charge>0) Psum+=P[4];
                        else Psum+=P[5];

                        if (Psum<=0.) Psum=1.;

                        result->setLHElectron(P[0]/Psum);
                        result->setLHMuon(P[1]/Psum);
                        result->setLHPion(P[2]/Psum);
                        result->setLHKaon(P[3]/Psum);
                        if (charge>0) result->setLHProton(P[4]/Psum);
                        else result->setLHProton(P[5]/Psum);

          }
  }

  return result;
}

bool
PndFsmMdtPid::detected(PndFsmTrack *t) const
{
    double theta = t->p4().Theta();
    double p=t->p4().Vect().Mag();
    double charge=t->charge();
    int type=abs(t->pdt());

        if (_useFlat) return ( charge!=0.0 && theta>=_thtMin && theta<=_thtMax && p>_pmin && ((type==13 && _rand->Rndm()<=_efficiency) || (type==211 && _rand->Rndm()<=_misId)) );

        return ( charge!=0.0 && theta>=_thtMin && theta<=_thtMax && p>_pmin );
}

void
PndFsmMdtPid::print(ostream &o)
{
  o <<"Parameters for detector <"<<detName()<<">"<<endl;
  o  <<"  _thtMin = "<<_thtMin<<endl;
  o  <<"  _thtMax = "<<_thtMax<<endl;
  o  <<"  _pmin = "<<_pmin<<endl;
  o  <<"  _misId = "<<_misId<<endl;
  o  <<"  _efficiency = "<<_efficiency<<endl;
}

void
PndFsmMdtPid::initParameters()
{
  _detName         = "MdtPid";
  _thtMin          = 5.;
  _thtMax          = 160.;
  _pmin            = 1.0;
  _misId                   = 0.01;
  _efficiency      = 1.0;
  _parFileName     = "$VMCWORKDIR/fastsim/MdtPidPdf3D.root";
  _useFlat                 = false;
}

bool
PndFsmMdtPid::setParameter(std::string &name, double value)
{
  // *****************
  // include here all parameters which should be settable via tcl
  // *****************

  bool knownName=true;

  if (name == "thtMin")
    _thtMin=value;
  else
  if (name == "thtMax")
    _thtMax=value;
  else
  if (name == "pmin")
        _pmin=value;
  else
  if (name == "misId")
        _misId=value;
  else
  if (name == "efficiency")
    _efficiency=value;
  else
    knownName=false;

  return knownName;
}

bool PndFsmMdtPid::readParameters()
{

  TFile *f=new TFile(_parFileName.c_str());

  for (int i=0;i<6;i++)
  {
    _mdtPidPdf[i]=0;
  }

  if (f->IsZombie())
  {
    cout <<" -W-  (PndFsmMdtPid::readParameters) - file "<<_parFileName.c_str()
         <<" doesn't exist. Using simple likelihood."<<endl;
        _useFlat = true;
  }
  else
  {
    _mdtPidPdf[0]=(TH3F*)f->Get("hpdf_e");
    _mdtPidPdf[1]=(TH3F*)f->Get("hpdf_mu");
    _mdtPidPdf[2]=(TH3F*)f->Get("hpdf_pi");
    _mdtPidPdf[3]=(TH3F*)f->Get("hpdf_k");
    _mdtPidPdf[4]=(TH3F*)f->Get("hpdf_p");
    _mdtPidPdf[5]=(TH3F*)f->Get("hpdf_pb");

        cout << _mdtPidPdf[0]<<endl;

    for (int i=0;i<6;i++) _mdtPidPdf[i]->SetDirectory(0);

    f->Close();
  }
  delete f;

  return true;
}