GEMAmplification.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#include "TPCSimulation/GEMAmplification.h"
#include <TStopwatch.h>
#include "MathUtils/CachingTF1.h"
#include <TFile.h>
#include "TPCBase/CDBInterface.h"
#include <fstream>
#include "Framework/Logger.h"
#include <filesystem>
 
using namespace o2::tpc;
using namespace o2::math_utils;
using boost::format;
 
GEMAmplification::GEMAmplification()
  : mRandomGaus(),
    mRandomFlat(RandomRing<>::RandomType::Flat),
    mGain()
{
  updateParameters();
 
  TStopwatch watch;
  watch.Start();
  const float sigmaOverMu = mGasParam->SigmaOverMu;
  const float kappa = 1 / (sigmaOverMu * sigmaOverMu);
  boost::format polya("1/(TMath::Gamma(%1%)*%2%) * TMath::Power(x/%3%, %4%) * TMath::Exp(-x/%5%)");
 
  const char* polyaFileName = "tpc_polya.root";
  TFile* outfile;
  auto cacheexists = std::filesystem::exists(polyaFileName);
  if (cacheexists) {
    LOG(info) << "TPC: GEM setup from existing cache";
    outfile = TFile::Open(polyaFileName);
  } else {
    LOG(info) << "TPC: GEM setup - initialization from scratch";
    outfile = TFile::Open(polyaFileName, "CREATE");
  }
 
  for (int i = 0; i < 4; ++i) {
    float s = mGEMParam->AbsoluteGain[i] / kappa;
    polya % kappa % s % s % (kappa - 1) % s;
    std::string name = polya.str();
    o2::math_utils::CachingTF1* polyaDistribution = nullptr;
    if (!cacheexists) {
      polyaDistribution = new o2::math_utils::CachingTF1("polya", name.c_str(), 0, 10.f * mGEMParam->AbsoluteGain[i]);
      polyaDistribution->SetNpx(100000);
    } else {
      polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get(TString::Format("func%d", i).Data());
      // FIXME: verify that distribution corresponds to the parameters used here
    }
    mGain[i].initialize(*polyaDistribution);
 
    if (!cacheexists) {
      outfile->WriteTObject(polyaDistribution, TString::Format("func%d", i).Data());
    }
    delete polyaDistribution;
  }
 
  const float gainStack = mGEMParam->TotalGainStack;
  const float kappaStack = mGEMParam->KappaStack;
  const float effStack = mGEMParam->EfficiencyStack;
  // Correct for electron losses occurring when changing kappa and the efficiency
  const float sStack = gainStack / (kappaStack * effStack);
  polya % kappaStack % sStack % sStack % (kappaStack - 1) % sStack;
  std::string name = polya.str();
  o2::math_utils::CachingTF1* polyaDistribution = nullptr;
  if (!cacheexists) {
    polyaDistribution = new o2::math_utils::CachingTF1("polya", name.c_str(), 0, 25.f * gainStack);
    polyaDistribution->SetNpx(50000);
  } else {
    polyaDistribution = (o2::math_utils::CachingTF1*)outfile->Get("polyaStack");
  }
  mGainFullStack.initialize(*polyaDistribution);
 
  if (!cacheexists) {
    outfile->WriteTObject(polyaDistribution, "polyaStack");
  }
  delete polyaDistribution;
 
  if (outfile) {
    outfile->Close();
  }
  watch.Stop();
  LOG(info) << "TPC: GEM setup (polya) took " << watch.CpuTime();
}
 
void GEMAmplification::updateParameters()
{
  auto& cdb = CDBInterface::instance();
  mGEMParam = &(ParameterGEM::Instance());
  mGasParam = &(ParameterGas::Instance());
  mGainMap = &(cdb.getGainMap());
}
 
int GEMAmplification::getStackAmplification(int nElectrons)
{
  const int nElectronsGEM1 = getSingleGEMAmplification(nElectrons, 0);
  const int nElectronsGEM2 = getSingleGEMAmplification(nElectronsGEM1, 1);
  const int nElectronsGEM3 = getSingleGEMAmplification(nElectronsGEM2, 2);
  const int nElectronsGEM4 = getSingleGEMAmplification(nElectronsGEM3, 3);
  return nElectronsGEM4;
}
 
int GEMAmplification::getEffectiveStackAmplification(int nElectrons)
{
  int nElectronsGEM = 0;
  for (int i = 0; i < nElectrons; ++i) {
    if (mRandomFlat.getNextValue() > mGEMParam->EfficiencyStack) {
      continue;
    }
    nElectronsGEM += mGainFullStack.getNextValue();
  }
  return nElectronsGEM;
}
 
int GEMAmplification::getSingleGEMAmplification(int nElectrons, int GEM)
{
  int collectionGEM = getElectronLosses(nElectrons, mGEMParam->CollectionEfficiency[GEM]);
  int amplificationGEM = getGEMMultiplication(collectionGEM, GEM);
  int extractionGEM = getElectronLosses(amplificationGEM, mGEMParam->ExtractionEfficiency[GEM]);
  return extractionGEM;
}
 
int GEMAmplification::getGEMMultiplication(int nElectrons, int GEM)
{
  if (nElectrons < 1) {
    return 0;
  } else if (nElectrons > 500) {
    return ((mRandomGaus.getNextValue() * std::sqrt(static_cast<float>(nElectrons)) *
             mGasParam->SigmaOverMu) +
            nElectrons) *
           mGEMParam->AbsoluteGain[GEM];
  } else {
    int electronsOut = 0;
    for (int i = 0; i < nElectrons; ++i) {
      electronsOut += mGain[GEM].getNextValue();
    }
    return electronsOut;
  }
}
 
int GEMAmplification::getElectronLosses(int nElectrons, float probability)
{
  float electronsFloat = static_cast<float>(nElectrons);
  if (nElectrons < 1 || probability < 0.00001) {
    return 0;
  } else if (probability > 0.99999) {
    return nElectrons;
  } else if (electronsFloat * probability >= 5.f && electronsFloat * (1.f - probability) >= 5.f) {
    return (mRandomGaus.getNextValue() * std::sqrt(electronsFloat * probability * (1 - probability))) +
           electronsFloat * probability + 0.5;
  } else {
    int nElectronsOut = 0;
    for (int i = 0; i < nElectrons; ++i) {
      if (mRandomFlat.getNextValue() < probability) {
        ++nElectronsOut;
      }
    }
    return nElectronsOut;
  }
}