CbmTrdHitProducerQa.cxx

/* Copyright (C) 2005-2019 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Matus Kalisky, Florian Uhlig, Denis Bertini [committer] */

// -----------------------------------------------------------------------
// -----                     CbmTrdHitProducerQa                     -----
// -----               Created 13/12/05  by M. Kalisky               -----
// -----------------------------------------------------------------------

#include "CbmTrdHitProducerQa.h"

#include "CbmDigiManager.h"
#include "CbmMCDataArray.h"
#include "CbmMCDataManager.h"
#include "CbmMCTrack.h"
#include "CbmMatch.h"
#include "CbmQaCanvas.h"
#include "CbmTrdDigi.h"
#include "CbmTrdHit.h"
#include "CbmTrdPoint.h"

#include "FairBaseParSet.h"
#include "FairDetector.h"
#include "FairRootManager.h"
#include "FairRunAna.h"
#include "FairRuntimeDb.h"

#include "TClonesArray.h"
#include "TH1F.h"
#include "TMath.h"
#include "TObjArray.h"

#include <iostream>
using std::cout;
using std::endl;

// ---- Default constructor -------------------------------------------------

CbmTrdHitProducerQa::CbmTrdHitProducerQa() : CbmTrdHitProducerQa("TrdHitProducerQa", "") {}
// --------------------------------------------------------------------------

// ---- Standard constructor ------------------------------------------------
CbmTrdHitProducerQa::CbmTrdHitProducerQa(const char* name, const char*)
  : FairTask(name)
  , fOutFolder("TrdHitProducerQA", "TrdHitProducerQA")
{
}
// --------------------------------------------------------------------------

// ---- Destructor ---------------------------------------------------------
CbmTrdHitProducerQa::~CbmTrdHitProducerQa() {}
// --------------------------------------------------------------------------

// ---- Initialisation ------------------------------------------------------
InitStatus CbmTrdHitProducerQa::Init()
{
  fOutFolder.Clear();
  histFolder = fOutFolder.AddFolder("hist", "Histogramms");

  const int nStations = fNoTrdStations * fNoTrdPerStation;

  for (int i = 0; i < nStations; i++) {
    fvhHitPullX.push_back(new TH1F(Form("L%iHitPullX", i), "", 500, -50, 50));
    fvhHitPullY.push_back(new TH1F(Form("L%iHitPullY", i), "", 500, -50, 50));
    fvhHitPullT.push_back(new TH1F(Form("L%iHitPullT", i), "", 500, -50, 50));
    fvhHitPullX[i]->SetCanExtend(TH1::kAllAxes);
    fvhHitPullY[i]->SetCanExtend(TH1::kAllAxes);
    fvhHitPullT[i]->SetCanExtend(TH1::kAllAxes);
    histFolder->Add(fvhHitPullX[i]);
    histFolder->Add(fvhHitPullY[i]);
    histFolder->Add(fvhHitPullT[i]);

    fvhHitResX.push_back(new TH1F(Form("L%iHitResX", i), "", 500, -50, 50));
    fvhHitResY.push_back(new TH1F(Form("L%iHitResY", i), "", 500, -50, 50));
    fvhHitResT.push_back(new TH1F(Form("L%iHitResT", i), "", 500, -50, 50));
    fvhHitResX[i]->SetCanExtend(TH1::kAllAxes);
    fvhHitResY[i]->SetCanExtend(TH1::kAllAxes);
    fvhHitResT[i]->SetCanExtend(TH1::kAllAxes);
    histFolder->Add(fvhHitResX[i]);
    histFolder->Add(fvhHitResY[i]);
    histFolder->Add(fvhHitResT[i]);

    fvhedEcut.push_back(new TH1F(Form("L%iedEcut", i), Form("dEdx of e- for layer %i, mom. cut", i), 600, 0., 60.));
    fvhedEall.push_back(new TH1F(Form("L%iedEall", i), Form("dEdx of e- for layer %i", i), 600, 0., 60.));

    fvhpidEcut.push_back(new TH1F(Form("L%ipidEcut", i), Form("dEdx of pi- for layer %i, mom. cut", i), 600, 0., 60.));
    fvhpidEall.push_back(new TH1F(Form("L%ipidEall", i), Form("dEdx of pi- for layer %i", i), 600, 0., 60.));

    histFolder->Add(fvhedEcut[i]);
    histFolder->Add(fvhedEall[i]);
    histFolder->Add(fvhpidEcut[i]);
    histFolder->Add(fvhpidEall[i]);

    fvdECanvas.push_back(new CbmQaCanvas(Form("cL%iEnergyLoss", i), Form("Energy Loss Layer %i", i), 2 * 400, 2 * 400));
    fvdECanvas[i]->Divide2D(4);
    fOutFolder.Add(fvdECanvas[i]);

    fvPullCanvas.push_back(
      new CbmQaCanvas(Form("cL%iPull", i), Form("Pull Distribution Layer %i", i), 3 * 400, 2 * 400));
    fvPullCanvas[i]->Divide2D(6);
    fOutFolder.Add(fvPullCanvas[i]);
  }

  // Get pointer to the ROOT I/O manager
  FairRootManager* rootMgr = FairRootManager::Instance();
  if (nullptr == rootMgr) {
    cout << "-E- CbmTrdHitProducerQa::Init : "
         << "ROOT manager is not instantiated !" << endl;
    return kFATAL;
  }

  // Get a pointer to the previous already existing data level
  fDigiMan = CbmDigiManager::Instance();
  if (nullptr == fDigiMan) {
    cout << "-W- CbmTrdHitProducerQa::Init : "
         << "no digi manager found !" << endl;
    return kERROR;
  }
  fDigiMan->Init();

  CbmMCDataManager* mcManager = (CbmMCDataManager*) FairRootManager::Instance()->GetObject("MCDataManager");

  // Get pointer to TRD point array
  fTrdPoints = mcManager->InitBranch("TrdPoint");
  if (nullptr == fTrdPoints) {
    cout << "-W- CbmTrdHitProducerQa::Init : "
         << "no TRD point array !" << endl;
    return kERROR;
  }

  // Get pointer to TRD hit array
  fTrdHitCollection = (TClonesArray*) rootMgr->GetObject("TrdHit");
  if (nullptr == fTrdHitCollection) {
    cout << "-W- CbmTrdHitProducerQa::Init : "
         << "no TRD hit array !" << endl;
    return kERROR;
  }

  // Get MCTrack array
  fMCTrackArray = mcManager->InitBranch("MCTrack");
  if (nullptr == fMCTrackArray) {
    cout << "-E- CbmTrdHitProducerQa::Init : No MCTrack array!" << endl;
    return kFATAL;
  }

  // Get pointer to TRD digi array match
  if (!fDigiMan->IsMatchPresent(ECbmModuleId::kTrd)) {
    cout << GetName() << ": no TRD match branch in digi manager." << endl;
    return kERROR;
  }

  return kSUCCESS;
}

// --------------------------------------------------------------------------


// ---- Task execution ------------------------------------------------------
void CbmTrdHitProducerQa::Exec(Option_t*)
{
  // Loop over TRD hits
  for (int iHit = 0; iHit < fTrdHitCollection->GetEntriesFast(); iHit++) {
    const CbmTrdHit* trdHit = (CbmTrdHit*) fTrdHitCollection->At(iHit);
    if (nullptr == trdHit) continue;

    const CbmMatch* trdDigiMatch = fDigiMan->GetMatch(ECbmModuleId::kTrd, trdHit->GetRefId());
    if (nullptr == trdDigiMatch) continue;

    const CbmTrdPoint* trdPoint =
      dynamic_cast<CbmTrdPoint*>(fTrdPoints->Get(trdDigiMatch->GetMatchedLink()));  // file, event, object
    if (nullptr == trdPoint) continue;

    const int planeId = trdHit->GetPlaneId();

    if (planeId >= fNoTrdStations * fNoTrdPerStation) {
      cout << GetName() << ": Warning, TRD plane out of bounds, skipping hit." << endl;
      continue;
    }

    //get particle ID of track corresponding to point
    const int fileId        = trdDigiMatch->GetMatchedLink().GetFile();
    const int event         = trdDigiMatch->GetMatchedLink().GetEntry();
    const int index         = trdPoint->GetTrackID();
    const CbmMCTrack* track = dynamic_cast<CbmMCTrack*>(fMCTrackArray->Get(fileId, event, index));
    const int partID        = track->GetPdgCode();

    const float momentum = TMath::Sqrt((trdPoint->GetPx() * trdPoint->GetPx()) + (trdPoint->GetPy() * trdPoint->GetPy())
                                       + (trdPoint->GetPz() * trdPoint->GetPz()));

    //electrons
    if ((TMath::Abs(partID) == 11) && (momentum > fMomCutLower) && (momentum < fMomCutUpper)) {
      fvhedEcut[planeId]->Fill((trdHit->GetELoss()) * 1000000);
    }
    if (TMath::Abs(partID) == 11) { fvhedEall[planeId]->Fill((trdHit->GetELoss()) * 1000000); }

    //pions
    if ((TMath::Abs(partID) == 211) && (momentum > fMomCutLower) && (momentum < fMomCutUpper)) {
      fvhpidEcut[planeId]->Fill((trdHit->GetELoss()) * 1000000);
    }
    if (TMath::Abs(partID) == 211) { fvhpidEall[planeId]->Fill((trdHit->GetELoss()) * 1000000); }

    // compute the Hit pulls for X and Y coordinate
    const float hitPosX   = trdHit->GetX();
    const float pointPosX = (trdPoint->GetX() + trdPoint->GetXOut()) / 2.;
    const float hitResX   = hitPosX - pointPosX;
    const float hitPullX  = hitResX / trdHit->GetDx();

    const float hitPosY   = trdHit->GetY();
    const float pointPosY = (trdPoint->GetY() + trdPoint->GetYOut()) / 2.;
    const float hitResY   = hitPosY - pointPosY;
    const float hitPullY  = hitResY / trdHit->GetDy();

    const double hitPosT   = trdHit->GetTime();
    const double pointPosT = trdPoint->GetTime();
    const double hitResT   = hitPosT - pointPosT;
    const double hitPullT  = hitResT / trdHit->GetTimeError();

    // fill histograms
    fvhHitPullX[planeId]->Fill(hitPullX);
    fvhHitPullY[planeId]->Fill(hitPullY);
    fvhHitPullT[planeId]->Fill(hitPullT);

    fvhHitResX[planeId]->Fill(hitResX);
    fvhHitResY[planeId]->Fill(hitResY);
    fvhHitResT[planeId]->Fill(hitResT);
  }
}
// --------------------------------------------------------------------------

// ---- Finish --------------------------------------------------------------
void CbmTrdHitProducerQa::Finish() { WriteHistograms(); }
// --------------------------------------------------------------------------

// ---- Write test histograms ------------------------------------------------
void CbmTrdHitProducerQa::WriteHistograms()
{
  for (size_t i = 0; i < fvdECanvas.size(); i++) {
    fvdECanvas[i]->cd(1);
    fvhedEcut[i]->DrawCopy("", "");

    fvdECanvas[i]->cd(2);
    fvhedEall[i]->DrawCopy("", "");

    fvdECanvas[i]->cd(3);
    fvhpidEcut[i]->DrawCopy("", "");

    fvdECanvas[i]->cd(4);
    fvhpidEall[i]->DrawCopy("", "");
  }

  for (size_t i = 0; i < fvPullCanvas.size(); i++) {
    fvPullCanvas[i]->cd(1);
    fvhHitPullX[i]->DrawCopy("", "");

    fvPullCanvas[i]->cd(2);
    fvhHitPullY[i]->DrawCopy("", "");

    fvPullCanvas[i]->cd(3);
    fvhHitPullT[i]->DrawCopy("", "");

    fvPullCanvas[i]->cd(4);
    fvhHitResX[i]->DrawCopy("", "");

    fvPullCanvas[i]->cd(5);
    fvhHitResY[i]->DrawCopy("", "");

    fvPullCanvas[i]->cd(6);
    fvhHitResT[i]->DrawCopy("", "");
  }

  FairSink* sink = FairRootManager::Instance()->GetSink();
  sink->WriteObject(&fOutFolder, nullptr);
}
// --------------------------------------------------------------------------

ClassImp(CbmTrdHitProducerQa);