CbmFastSim.cxx

Go to the documentation of this file.

/* Copyright (C) 2016 Justus-Liebig-Universitaet Giessen, Giessen
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Julian Book [committer] */
 
//
//  FastSim
//
 
//C++ class headers
#include <fstream>
#include <string>
 
//Fair class headers
#include "FairRootManager.h"
#include <Logger.h>
//#include "FairRunAna.h"
#include "FairRun.h"
#include "FairRuntimeDb.h"
 
//CBM class headers
#include "CbmMCTrack.h"
#include "CbmStack.h"
 
#include "PairAnalysisHelper.h"
 
//ROOT class headers
#include "CbmFastSim.h"
 
#include <TClonesArray.h>
#include <TDatabasePDG.h>
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <THnBase.h>
#include <THnSparse.h>
#include <TMath.h>
#include <TMatrixD.h>
#include <TPDGCode.h>
#include <TParticle.h>
#include <TParticlePDG.h>
#include <TProfile3D.h>
#include <TRandom3.h>
#include <TVector3.h>
#include <TVectorD.h>
#include <TVirtualMC.h>
 
 
using std::cout;
using std::endl;
using std::ifstream;
using std::string;
 
 
// -----   Default constructor   -------------------------------------------
CbmFastSim::CbmFastSim(bool persist)
  : FairTask("Cbm Fast Simulation")
  , fRand(new TRandom3(0))
  , fdbPdg(TDatabasePDG::Instance())
{
  for (Int_t idx = 0; idx < fgkNParts; idx++) {
    fEFF[idx]      = NULL;
    fEFFdenom[idx] = NULL;
    for (Int_t j = 0; j < 10; j++)
      fEFFproj[idx][j] = NULL;
 
    fRFp[idx] = NULL;
  }
}
// -------------------------------------------------------------------------
 
// -----   Destructor   ----------------------------------------------------
CbmFastSim::~CbmFastSim()
{
  FairRootManager* fManager = FairRootManager::Instance();
  fManager->Write();
 
  if (fFastTracks) {
    fFastTracks->Delete();
    delete fFastTracks;
  }
  if (fRand) delete fRand;
}
// -------------------------------------------------------------------------
 
 
void CbmFastSim::Register()
{
  fFastTracks = new TClonesArray("TParticle");
  FairRootManager::Instance()->Register("FastTrack", "FastSim", fFastTracks, kTRUE);
}
 
 
// -----   Public method Init   --------------------------------------------
void CbmFastSim::InitTask() { Init(); }
 
// -----   Public method Init   --------------------------------------------
InitStatus CbmFastSim::Init()
{
 
  for (Int_t idx = 0; idx < fgkNParts; idx++) {  
 
    if (fEFF[idx] && fEFFdenom[idx]) {
      Info("Init", "Prepare efficiency map and projections for index %d and method %d", idx, fMethod);
      //      fEFF[idx]->Print("a");
 
      // first get projections and calculate integrated effiencies
      // NOTE: if statistic are bad these are used
      for (Int_t j = 0; j < fEFF[idx]->GetNdimensions(); j++) {
 
        switch (fMethod) {
          case kIgnoreFluct:
            continue;
            break;
          case kLastDim:
            if (j < (fEFF[idx]->GetNdimensions() - 1)) continue;
          case kAverage:
          case kFactorize: {
            fEFFproj[idx][j] = fEFF[idx]->Projection(j, "E");
            fEFFproj[idx][j]->SetName(Form("NOMidx%dj%d", idx, j));
 
            TH1* den = fEFFdenom[idx]->Projection(j, "E");
            den->SetName(Form("DENidx%dj%d", idx, j));
            fEFFproj[idx][j]->Divide(den);
 
            // factoring using last dimension!=1 histograms
            if (fMethod == kFactorize && j > 0)
              fEFFproj[idx][j]->Scale(1. / fEFFproj[idx][j]->GetBinContent(fEFFproj[idx][j]->GetMaximumBin()));
 
            if (den) delete den;
          } break;
          case kInterpolate: {
            if (j > 0) continue;
            fEFFproj[idx][j] = fEFF[idx]->Projection(j, 1, 2, "E");
            fEFFproj[idx][j]->SetName(Form("NOMidx%dj%d12", idx, j));
 
            TH1* den = fEFFdenom[idx]->Projection(j, 1, 2, "E");
            den->SetName(Form("DENidx%dj%d12", idx, j));
            fEFFproj[idx][j]->Divide(den);
 
            if (den) delete den;
          } break;
          default: break;
        }
      }
 
      // calculate main effiency map
      fEFF[idx]->Divide(fEFFdenom[idx]);
 
    }  // endif: nominator
 
 
    if (fRFp[idx]) {
      Info("Init", "Momentum smearing map found for index %d", idx);
      PairAnalysisHelper::CumulateSlicesX(fRFp[idx], kFALSE, kTRUE);
    }
 
  }  //end: particle species
 
  // Create and register output array
  Register();
 
  return kSUCCESS;
}
 
void CbmFastSim::SetParContainers()
{
 
  // Get run and runtime database
  FairRun* run = FairRun::Instance();
  if (!run) Fatal("SetParContainers", "No analysis run");
 
  FairRuntimeDb* db = run->GetRuntimeDb();
  if (!db) Fatal("SetParContainers", "No runtime database");
}
// -------------------------------------------------------------------------
 
void CbmFastSim::SetSeed(unsigned int seed)
{
  fRand->SetSeed(seed);
}
 
// -----   Public method Finish   --------------------------------------------
void CbmFastSim::Finish()
{
}
// -------------------------------------------------------------------------
 
// -----   Public method Exec   --------------------------------------------
void CbmFastSim::Exec(Option_t* opt)
{
 
  CbmStack* fStack = (CbmStack*) gMC->GetStack();
  Int_t nTracks    = fStack->GetNtrack();
  //  Info("Exec","number of tracks **** %d",Tracks);
 
  if (fFastTracks->GetEntriesFast() != 0) fFastTracks->Clear("C");
 
  TClonesArray& chrgCandidates = *fFastTracks;  // Charged Candidates
 
  Int_t iKeep = 0;
  Int_t iConv = 0;
 
  Double_t* vals = new Double_t[4];  //fEFF->GetNdimensions()];
  Int_t* coord   = new Int_t[4];
 
  for (Int_t iTrack = 0; iTrack < nTracks; iTrack++) {
 
    TParticle* t = fStack->GetParticle(iTrack);
    int pdg      = t->GetPdgCode();
    //    t->Print();
 
    if (!fdbPdg->GetParticle(pdg)) continue;
    if (TMath::Abs(fdbPdg->GetParticle(pdg)->Charge()) < 3) continue;
 
    vals[0] = t->Pt();
    vals[1] = t->Theta();
    vals[2] = TVector2::Phi_0_2pi(t->Phi());
    vals[3] = t->Vz();
 
    if (TMath::Abs(pdg) == 11 && fStack->GetParticle(t->GetMother(0))
        && fStack->GetParticle(t->GetMother(0))->GetPdgCode() == 22) {
      if (!PassEfficiencyFilter(22, vals, coord)) continue;
      iConv++;
    }
    else {
      if (!PassEfficiencyFilter(pdg, vals, coord)) continue;
    }
 
    TLorentzVector p4(t->Px(), t->Py(), t->Pz(), t->Energy());
 
    //    TVector3 stvtx(t->Vx(),t->Vy(),t->Vz());
 
    if (!Smearing(&p4, pdg)) { /*ERROR*/
      continue;
    }
 
    new ((*fFastTracks)[iKeep]) TParticle(pdg, t->GetStatusCode(), iTrack, 0, 0, 0, p4.Px(), p4.Py(), p4.Pz(),
                                          p4.Energy(), t->Vx(), t->Vy(), t->Vz(), 0);
 
    iKeep++;
 
  }  //trackloop
 
  // cleanup
  delete[] vals;
  delete[] coord;
 
  Info("Exec", "stack size: \t %d \n", fStack->GetNtrack());
  Info("Exec", "kept tracks: \t %d \n", iKeep);
  Info("Exec", "kept converison tracks: \t %d \n", iConv);
}
 
Bool_t CbmFastSim::PassEfficiencyFilter(Int_t pdg, Double_t* vals, Int_t* coord)
{
  //
  // check if partilce passed the reconstruction cuts
  //
  Int_t idx = -1;
  switch (TMath::Abs(pdg)) {
    case kGamma: idx = CbmFastSim::kGam; break;
    case kElectron: idx = CbmFastSim::kEle; break;
    case kMuonMinus: idx = CbmFastSim::kMuo; break;
    case kPiPlus: idx = CbmFastSim::kPio; break;
    case kKPlus: idx = CbmFastSim::kKao; break;
    case kProton: idx = CbmFastSim::kPro; break;
    default: /*Warning("PassEfficiencyFilter","Pdg code %d not supported, particle will not be processed",pdg);*/
      return kFALSE;
  }
 
  if (!fEFF[idx]) {
    Warning("PassEfficiencyFilter", "No filter for pdg code %d", pdg);
    return kFALSE;
  }
 
  //  Info("PassEfficiencyFilter","Efficiency for pt:%.2e theta:%.2e phi:%.2e", vals[0],vals[1],vals[2]);
 
  Long_t bin = fEFF[idx]->GetBin(vals, kFALSE);
 
  Double_t eff = 0.0;
  Double_t err = 0.0;
  if (bin > 0) {
    eff = fEFF[idx]->GetBinContent(bin, coord);  // write bin coordinates into coord
 
    if (fMethod != kIgnoreFluct) {
 
      err = fEFF[idx]->GetBinError(bin);
      //    if (eff>0. && err/eff>=0.60 ) return kFALSE;
      if (eff > 0. && err / eff >= 0.10) {
        Double_t effm = GetEfficiency(idx, vals, coord);  // vals was NULL
        //      Info("PassEfficiencyFilter","statistical error %f too small (eff:%.2e, new:%.2e)",(eff>0.?err/eff:0.),eff,effm);
        eff = effm;
      }
    }
  }
  //  else {
  //    eff= GetEfficiency(idx, vals, NULL);
  //    Info("PassEfficiencyFilter","No entry in map, caculate via projections: %.2e (mean eff.)",eff);
  //  }
 
  Double_t rndm = fRand->Rndm();
 
  if (rndm > eff) return kFALSE;
  else
    return kTRUE;
}
 
Bool_t CbmFastSim::Smearing(TLorentzVector* p4, Int_t pdg)
{
 
  Int_t idx = -1;
  switch (TMath::Abs(pdg)) {
    case kElectron: idx = CbmFastSim::kEle; break;
    case kMuonMinus: idx = CbmFastSim::kMuo; break;
    case kPiPlus: idx = CbmFastSim::kPio; break;
    case kKPlus: idx = CbmFastSim::kKao; break;
    case kProton: idx = CbmFastSim::kPro; break;
    default: return kFALSE;
  }
 
  Double_t tP     = p4->P();
  Double_t tTheta = p4->Theta();
  Double_t tPhi   = TVector2::Phi_0_2pi(p4->Phi());
  Double_t mass   = fdbPdg->GetParticle(pdg)->Mass();
 
  Double_t sP     = AnalyzeMap(fRFp[idx], tP);
  Double_t sTheta = tTheta;
  Double_t sPhi   = tPhi;
 
  Double_t sPx = sP * TMath::Sin(sTheta) * TMath::Cos(sPhi);
  Double_t sPy = sP * TMath::Sin(sTheta) * TMath::Sin(sPhi);
  Double_t sPz = sP * TMath::Cos(sTheta);
 
  Double_t sPt = TMath::Sqrt(sPx * sPx + sPy * sPy);
  Double_t eta = -TMath::Log(TMath::Tan(sTheta / 2));
 
 
  //  p4->Print();
  p4->SetPtEtaPhiM(sPt, eta, sPhi, mass);
  // p4->Print();
  // printf("\n\n");
 
  return kTRUE;
}
 
Double_t CbmFastSim::AnalyzeMap(TH2F* map, Double_t refValue)
{
 
  if (!map) return refValue;
 
  TAxis* xRec = map->GetXaxis();
  TAxis* yGen = map->GetYaxis();
 
  Int_t binGen = yGen->FindBin(refValue);
 
  if (binGen < 1 || binGen > yGen->GetNbins()) { return (fRand->Gaus(refValue, refValue * 0.02)); }
 
  Float_t* arrayh = map->GetArray();
  Int_t offset    = (yGen->GetNbins() + 2) * (binGen) + 1;
 
  Float_t r1   = fRand->Rndm();
  Int_t binRec = TMath::BinarySearch(xRec->GetNbins() + 2, arrayh + offset - 1, r1);
 
  Double_t smearedValue = xRec->GetBinLowEdge(binRec + 1);
  if (r1 > arrayh[offset - 1 + binRec])
    smearedValue += xRec->GetBinWidth(binRec + 1) * (r1 - arrayh[offset - 1 + binRec])
                    / (arrayh[offset - 1 + binRec + 1] - arrayh[offset - 1 + binRec]);
 
  // printf("binGen: %d, binRec: %d , offset: %d , \t r1: %f \t in: %f --> out: %f \n",
  //     binGen,binRec,offset, r1,refValue,smearedValue);
 
  return smearedValue;
}
 
 
void CbmFastSim::SetLookupEfficiency(THnBase* nom, THnBase* denom, EParticleType part)
{
 
  if (!nom || !denom) Fatal("SetLookupEfficiency", "Either nominator or denominator is NULL ");
  fEFF[part]      = nom;
  fEFFdenom[part] = denom;
  /*
  // first get projections and calculate integrated effiencies
  // NOTE: if statistic are bad these are used
  for (Int_t j = 0; j < nom->GetNdimensions(); ++j) {
 
    switch(fMethod) {
    case kIgnoreFluct:
      continue;
      break;
    case kLastDim:
      if( j<(nom->GetNdimensions()-1) ) continue;
      break;
    case kAverage:
    case kFactorize: {
      fEFFproj[part][j] = nom->Projection(j,"E");
      fEFFproj[part][j]->SetName(Form("NOMidx%dj%d",part,j));
 
      TH1 *den=denom->Projection(j,"E");
      den->SetName(Form("DENidx%dj%d",part,j));
      fEFFproj[part][j]->Divide( den );
 
      // factoring using last dimension!=1 histograms
      if(fMethod==kFactorize && j>0)    fEFFproj[part][j]->Scale( 1./fEFFproj[part][j]->GetBinContent( fEFFproj[part][j]->GetMaximumBin() ) );
 
      if(den) delete den;
    }
      break;
    case kInterpolate: {
      if(j>0 && j<3)  continue;
      fEFFproj[part][j] = nom->Projection(j,1,2,"E");
      fEFFproj[part][j]->SetName(Form("NOMidx%dj%d12",part,j));
 
      TH1 *den=denom->Projection(j,1,2,"E");
      den->SetName(Form("DENidx%dj%d12",part,j));
      fEFFproj[part][j]->Divide( den );
 
      if(den) delete den;
    }
      break;
    default:
      break;
    }
 
 
  }
 
  // calculate main effiency map
  nom->Divide(denom);
  */
}
 
Double_t CbmFastSim::GetEfficiency(Int_t idx, Double_t* vals, Int_t* coord)
{
 
  if (fMethod == kIgnoreFluct) return 0.;
 
 
  Int_t ndim       = 0;
  Double_t meanEff = (fMethod == kFactorize ? 1. : 0.);
  Double_t eff     = 0.;
  //  for (Int_t j = 0; j < 10; ++j) {
  for (Int_t j = 10; j >= 0; j--) {
 
    if (!fEFFproj[idx][j]) continue;
 
    if (coord) eff = fEFFproj[idx][j]->GetBinContent(coord[j]);
    else
      eff = fEFFproj[idx][j]->GetBinContent(fEFFproj[idx][j]->FindBin(vals[j]));
    //    Info("GetEfficiency","int eff for j:%d and idx:%d and value: %.2e is %.2e!",j,idx,vals[j],eff);
 
    switch (fMethod) {
      case kIgnoreFluct: return 0.;
      case kInterpolate: {
        if (coord) {
          for (Int_t c = 0; c < 3; c++) {
            if (coord[c] == 0) coord[c] += 2;
            if (coord[c] == 1) coord[c] += 1;
            if (coord[c] == fEFF[idx]->GetAxis(c)->GetNbins()) coord[c] -= 2;
            if (coord[c] == fEFF[idx]->GetAxis(c)->GetNbins() - 1) coord[c] -= 1;
          }
          eff = fEFFproj[idx][j]->Interpolate(fEFFproj[idx][j]->GetXaxis()->GetBinCenter(coord[0]),
                                              fEFFproj[idx][j]->GetYaxis()->GetBinCenter(coord[1]),
                                              fEFFproj[idx][j]->GetZaxis()->GetBinCenter(coord[2]));
        }
        else
          eff = fEFFproj[idx][j]->Interpolate(vals[0], vals[1], vals[2]);
        return eff;
      } break;
      case kAverage:
      case kLastDim:
        if (eff < 1.e-10) return 0.;  
        meanEff += eff;
        ndim++;
        break;
      case kFactorize:
        if (eff < 1.e-10) return 0.;  
        meanEff *= eff;
        ndim = 1;  //      ndim++;
        break;
    }
  }
 
  // protect against division by 0
  if (!ndim) {
    Error("GetEfficiency", "number of dimensions for idx:%d is zero!", idx);
    return 0.;
  }
 
  // return mean efficiency
  return (meanEff / ndim);
}
 
Bool_t CbmFastSim::IsSelected(TObject* sel, Int_t opt)
{
  //
  // check if partilce passed the cuts
  //
 
  Int_t pdg     = opt;
  TVector3* mom = dynamic_cast<TVector3*>(sel);
 
  // Info("IsSelected","Compare pdg code %d to e(%d),mu(%d),pi(%d),K(%d),p(%d)",
  //      TMath::Abs(pdg),kElectron,kMuonMinus,kPiPlus,kKPlus,kProton);
 
  Int_t idx = -1;
  switch (TMath::Abs(pdg)) {
    case kGamma: idx = CbmFastSim::kGam; break;
    case kElectron: idx = CbmFastSim::kEle; break;
    case kMuonMinus: idx = CbmFastSim::kMuo; break;
    case kPiPlus: idx = CbmFastSim::kPio; break;
    case kKPlus: idx = CbmFastSim::kKao; break;
    case kProton: idx = CbmFastSim::kPro; break;
    default: /*Warning("IsSelected","Pdg code %d not supported, particle will not be processed",pdg);*/ return kFALSE;
  }
 
  // Double_t rndm2 = fRand->Rndm();
  // if(rndm2<0.1) return kTRUE;
  // else         return kFALSE;
 
  if (!fEFF[idx]) {
    Warning("IsSelected", "No filter for pdg code %d", pdg);
    return kTRUE;
  }
 
  Double_t* vals = new Double_t[3];
 
  vals[0] = mom->Pt();
  vals[1] = mom->Theta();
  vals[2] = TVector2::Phi_0_2pi(mom->Phi());
 
  Long_t bin = fEFF[idx]->GetBin(vals, kFALSE);
 
  delete vals;
 
  Double_t eff = 0.0;
  if (bin > 0) eff = fEFF[idx]->GetBinContent(bin);
  Double_t rndm = fRand->Rndm();
 
 
  if (rndm > eff) return kFALSE;
  else
    return kTRUE;
}
 
ClassImp(CbmFastSim)